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.
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;
Anonymous function
Optimization.CustomCostFcn = @(X,U,e,data,params) myCostFunction(X,U,e,data,params);
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 commaseparated 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)byN_{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)byN_{u} array.
The final row of U is always a duplicate of the preceding row;
that is, U(end,:) = U(end1,:) . 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. If you have nonlinear soft constraints defined in your
inequality constraint function (  
data  Input  Additional signals, specified as a structure with the following fields:
 
params  Input  Optional parameters, specified as a commaseparated list (for example
If your model uses optional
parameters, you must specify the number of parameters using the
 
J  Output  Computed cost, returned as a scalar 
Your custom cost function must:
Be a continuous, finite function of
U
,X
, ande
and have finite first derivativesIncrease 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(sum(U1.^2)) + 5*sum(sum(X2.^3) + sum(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(sum((U(1:end1,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 function in the current working folder or on the MATLAB path
Jacobian.CustomCostFcn = @myCostJacobian;
Anonymous function
Jacobian.CustomCostFcn = @(X,U,e,data,params) myCostJacobian(X,U,e,data,params)
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 commaseparated 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 pbyN_{x} array,
where $$\text{G}\left(i,j\right)=\partial \text{J}/\partial \text{X}\left(i+1,j\right)$$. 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
pbyN_{mv} array,
where $$\text{Gmv}\left(i,j\right)=\partial \text{J}/\partial \text{U}\left(i,MV\left(j\right)\right)$$ and MV(j) is the
jth MV index in
Since the controller forces

Ge  Jacobian of the cost function with respect to the slack variable,
e , returned as a scalar, where $$\text{Ge}=\partial \text{J}/\partial \text{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_{1} 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
.