Specify Cost Function for Nonlinear MPC

While traditional linear MPC controllers optimize control actions to minimize a quadratic cost function, nonlinear MPC controllers support generic custom cost functions. For example, you can specify your cost function as a combination of linear or nonlinear functions of the system states and inputs. To improve computational efficiency, you can also specify an analytical Jacobian for your custom cost function.

Using a custom cost function, you can, for example:

Maximize profitability
Minimize energy consumption

When you specify a custom cost function for your nonlinear MPC controller, you can choose to either replace or augment the standard quadratic MPC cost function. By default, an nlmpc controller replaces the standard cost function with your custom cost function. In this case, the controller ignores the standard tuning weights in its Weights property.

To use an objective function that is the sum of the standard costs and your custom costs, set the Optimization.ReplaceStandardCost property of your nlmpc object to false. In this case, the standard tuning weights specified in the Weights property of the controller contribute to the cost function. However, you can eliminate any of the standard cost function terms by setting the corresponding penalty weight to zero. For more information on the standard MPC cost function, see Standard Cost Function.

For a multistage MPC controller, each stage has its own cost function, which is a function only of the decision variables and parameters at that stage.

Before simulating your controller, it is best practice to validate your custom functions, including the cost function and its Jacobian, using the validateFcns command.

Custom Cost Function

To configure your nonlinear MPC controller to use a custom cost function, set its Optimization.CustomCostFcn property to one of the following.

Name of a function in the current working folder or on the MATLAB^® path, specified as a string or character vector
```
Optimization.CustomCostFcn = "myCostFunction";
```
Handle to a function in the current working folder or on the MATLAB path
```
Optimization.CustomCostFcn = @myCostFunction;
```
For more information on local functions, see Local Functions.
Anonymous function
```
Optimization.CustomCostFcn = @(X,U,e,data,params) myCostFunction(X,U,e,data,params);
```
For more information on anonymous functions, see Anonymous Functions.

Note

Only functions defined in a separate file in the current folder or on the MATLAB path are supported for C/C++ code generation. Therefore, specifying state, output, cost, or constraint functions (or their Jacobians) as local or anonymous functions is not recommended.

Your custom cost function must have one of the following signatures.

If your controller does not use optional parameters:
```
function J = myCostFunction(X,U,e,data)
```
If your controller uses parameters. Here, params is a comma-separated list of parameters:
```
function J = myCostFunction(X,U,e,data,params)
```

This table describes the inputs and outputs of this function, where:

N_x is the number of states and is equal to the Dimensions.NumberOfStates property of the controller.
N_u is the number of inputs, including all manipulated variables, measured disturbances, and unmeasured disturbances, and is equal to the Dimensions.NumberOfInputs property of the controller.
p is the prediction horizon.
k is the current time.

Argument Input/Output Description

X Input State trajectory from time k to time k+p, specified as a (p+1)-by-N_x array. The first row of X contains the current state values, which means that the solver does not use the values in X(1,:) as decision variables during optimization.

U Input Input trajectory from time k to time k+p, specified as a (p+1)-by-N_u array. The final row of U is always a duplicate of the preceding row; that is, U(end,:) = U(end-1,:). Therefore, the values in the final row of U are not independent decision variables during optimization.

e

Input

Slack variable for constraint softening, specified as a nonnegative scalar. e is zero if there are no soft constraints in your controller.

If you have nonlinear soft constraints defined in your inequality constraint function (Model.CustomIneqConFcn), use a positive penalty weight on e and make them part of the cost function.

data

Input

Additional signals, specified as a structure with the following fields:

Field	Description
`Ts`	Prediction model sample time, as defined in the `Ts` property of the controller
`CurrentStates`	Current prediction model states, as specified in the `x` input argument of `nlmpcmove`
`LastMV`	MV moves used in previous control interval, as specified in the `lastmv` input argument of `nlmpcmove`
`References`	Reference values for plant outputs, as specified in the `ref` input argument of `nlmpcmove`
`MVTarget`	Manipulated variable targets, as specified in the `MVTarget` property of an `nlmpcmoveopt` object
`PredictionHorizon`	Prediction horizon, as defined in the `PredictionHorizon` property of the controller
`NumOfStates`	Number of states, as defined in the `Dimensions.NumberOfStates` property of the controller
`NumOfOutputs`	Number of outputs, as defined in the `Dimensions.NumberOfOutputs` property of the controller
`NumOfInputs`	Number of inputs, as defined in the `Dimensions.NumberOfInputs` property of the controller
`MVIndex`	Manipulated variables indices, as defined in the `Dimensions.MVIndex` property of the controller
`MDIndex`	Measured disturbance indices, as defined in the `Dimensions.MDIndex` property of the controller
`UDIndex`	Unmeasured disturbance indices, as defined in the `Dimensions.UDIndex` property of the controller

params

Input

Optional parameters, specified as a comma-separated list (for example p1,p2,p3). The same parameters are passed to the prediction model, custom cost function, and custom constraint functions of the controller. For example, if the state function uses only parameter p1, the constraint functions use only parameter p2, and the cost function uses only parameter p3, then all three parameters are passed to all of these functions.

If your model uses optional parameters, you must specify the number of parameters using the Model.NumberOfParameters property of the controller.

J Output Computed cost, returned as a scalar

Your custom cost function must:

Be a continuous, finite function of U, X, and e and have finite first derivatives
Increase as the slack variable e increases or be independent of it

To use output variable values in your cost function, you must first derive them from the state and input arguments using the prediction model output function, as specified in the Model.OutputFcn property of the controller. For example, to compute the output trajectory Y from time k to time k+p, use:

p = data.PredictionHorizon;
for i=1:p+1
    Y(i,:) = myOutputFunction(X(i,:)',U(i,:)',params)';
end

For more information on the prediction model output function, see Specify Prediction Model for Nonlinear MPC.

Typically, you optimize control actions to minimize the cost function across the prediction horizon. Since the cost function value must be a scalar, you compute the cost function at each prediction horizon step and add the results together. For example, suppose that the stage cost function is:

$J = 10 u_{1}^{2} + 5 x_{2}^{3} + x_{1}$

That is, you want to minimize the difference between the first output and its reference value, and the product of the first manipulated variable and the second state. To compute the total cost function across the prediction horizon, use:

p = data.PredictionHorizon;
U1 = U(1:p,data.MVIndex(1));
X1 = X(2:p+1,1);
X2 = X(2:p+1,2);
J = 10*sum(U1.^2) + 5*sum(X2.^3) + sum(X1);

In general, for cost functions, do not use the following values, since they are not part of the decision variables used by the solver:

U(end,:) — This row is a duplicate of the preceding row.
X(1,:) — This row contains the current state values.

Since this example cost function is relatively simple, you can specify it using an anonymous function handle.

For relatively simple costs, you can specify the cost function using an anonymous function handle. For example, to specify an anonymous function that implements just the first term of the preceding cost function, use:

Optimization.CustomCostFcn = @(X,U,data) 10*sum(U(1:end-1,data.MVIndex(1)).^2);

Cost Function Jacobian

To improve computational efficiency, it is best practice to specify an analytical Jacobian for your custom cost function. If you do not specify a Jacobian, the controller computes the Jacobian using numerical perturbation. To specify a Jacobian for your cost function, set the Jacobian.CustomCostFcn property of the controller to one of the following.

Name of a function in the current working folder or on the MATLAB path, specified as a string or character vector
```
Jacobian.CustomCostFcn = "myCostJacobian";
```
Handle to a local function, or a function defined in the current working folder or on the MATLAB path
```
Jacobian.CustomCostFcn = @myCostJacobian;
```
For more information on local functions, see Local Functions.
Anonymous function
```
Jacobian.CustomCostFcn = @(X,U,e,data,params) myCostJacobian(X,U,e,data,params)
```
For more information on anonymous functions, see Anonymous Functions.

Note

Your cost Jacobian function must have one of the following signatures.

If your controller does not use optional parameters:
```
function [G,Gmv,Ge] = myCostJacobian(X,U,e,data)
```
If your controller uses parameters. Here, params is a comma-separated list of parameters:
```
function [G,Gmv,Ge] = myCostJacobian(X,U,e,data,params)
```

The input arguments of the cost Jacobian function are the same as the inputs of the custom cost function. This table describes the outputs of the Jacobian function, where:

N_x is the number of states and is equal to the Dimensions.NumberOfStates property of the controller.
N_mv is the number of manipulated variables.
p is the prediction horizon.

Argument Description

G Jacobian of the cost function with respect to the state trajectories, returned as a p-by-N_x array, where $G (i, j) = \partial J / \partial X (i + 1, j)$ . Compute G based on X from the second row to row p+1, ignoring the first row.

Argument	Description
`G`	Jacobian of the cost function with respect to the state trajectories, returned as a p-by-N_x array, where $G (i, j) = \partial J / \partial X (i + 1, j)$ . Compute `G` based on `X` from the second row to row p+1, ignoring the first row.
`Gmv`	Jacobian of the cost function with respect to the manipulated variable trajectories, returned as a p-by-N_mv array, where $Gmv (i, j) = \partial J / \partial U (i, M V (j))$ and MV(j) is the jth MV index in `data.MVIndex`. Since the controller forces `U(p+1,:)` to equal `U(p,:)`, if your cost function uses `U(p+1,:)`, you must include the impact of both `U(p,:)` and `U(p+1,:)` in the Jacobian for `U(p,:)`.
`Ge`	Jacobian of the cost function with respect to the slack variable, `e`, returned as a scalar, where $Ge = \partial J / \partial e$ .

Gmv

Jacobian of the cost function with respect to the manipulated variable trajectories, returned as a p-by-N_mv array, where $Gmv (i, j) = \partial J / \partial U (i, M V (j))$ and MV(j) is the jth MV index in data.MVIndex.

Since the controller forces U(p+1,:) to equal U(p,:), if your cost function uses U(p+1,:), you must include the impact of both U(p,:) and U(p+1,:) in the Jacobian for U(p,:).

Ge Jacobian of the cost function with respect to the slack variable, e, returned as a scalar, where $Ge = \partial J / \partial e$ .

To use output variable values and their Jacobians in your cost Jacobian function, you must first derive them from the state and input arguments. To do so, use the Jacobian of the prediction model output function, as specified in the Jacobian.OutputFcn property of the controller. For example, to compute the output variables Y and their Jacobians Yjacob from time k to time k+p, use:

p = data.PredictionHorizon;
for i=1:p+1
    Y(i,:) = myOutputFunction(X(i,:)',U(i,:)',params)';
end
for i=1:p+1
    Yjacob(i,:) = myOutputJacobian(X(i,:)',U(i,:)',params)';
end

Since prediction model output functions do not support direct feedthrough from inputs to outputs, the output function Jacobian contains partial derivatives with respect to only the states in X. For more information on the output function Jacobian, see Specify Prediction Model for Nonlinear MPC.

To find the Jacobians, compute the partial derivatives of the cost function with respect to the state trajectories, manipulated variable trajectories, and slack variable. For example, suppose that your cost function is as follows, where u₁ is the first manipulated variable.

$J = 10 u_{1}^{2} + 5 x_{2}^{3} + x_{1}$

To compute the Jacobian with respect to the state trajectories, use the following. Recall that you compute G based on X from the second row to row p+1, ignoring the first row.

p = data.PredictionHorizon;
Nx = data.NumOfStates;
U1 = U(1:p,data.MVIndex(1));
X2 = X(2:p+1,2);

G = zeros(p,Nx);
G(1:p,1) = 1;
G(1:p,2) = 15*X2.^2;

To compute the Jacobian with respect to the manipulated variable trajectories, use:

Nmv = length(data.MVIndex);

Gmv = zeros(p,Nmv);
Gmv(1:p,1) = 20*U1;

In this case, the derivative with respect to the slack variable is Ge = 0.

Note

For multistage nonlinear MPC, you can automatically generate a stage cost Jacobian function using generateJacobianFunction.

Related Examples

Economic MPC Control of Ethylene Oxide Production