6. MathWorks Nonlinear MPC Plugin

6.1. Introduction

As a result of a long-term collaboration, MathWorks Inc. and Embotech AG have extended the Model Predictive Control Toolbox™ with a plugin for the FORCESPRO nonlinear solvers. Users are now able to use the FORCESPRO nonlinear interior-point (IP) and sequential quadratic programming (SQP) solvers in MATLAB® and Simulink® from within the MATLAB® Model Predictive Control Toolbox within the nonlinear MPC API. This plugin leverages the powerful design capabilities of the Model Predictive Control Toolbox™ and the computational performance of FORCESPRO. FORCESPRO extends the Model Predictive Control Toolbox with code-generated IP and SQP solvers that are not based on finite-difference derivatives computation, resulting in faster convergence. Thanks to FORCESPRO, the nonlinear API now comes with two classes of nonlinear solvers compatible with code generation that can be deployed to various real-time targets.

The FORCESPRO nonlinear MPC plugin consists of the following two API methods, which are covered in details later:

  • nlmpcToForces generates a FORCESPRO nonlinear solver from a nonlinear MPC (NMPC) object designed by the Model Predictive Control Toolbox

  • nlmpcmoveForces calls the generated solver on a specific NMPC problem instance

The nonlinear plugin also comes with a Simulink® library that enables users to run the FORCESPRO solvers from within their Simulink® models. It is compatible with MATLAB R2020a and R2020b.

This interface is provided with Variant L and partially with Variant M of FORCESPRO.

6.2. Defining a nonlinear model

In order to call the FORCESPRO code generation, a nonlinear MPC object needs to be built from a Model object. The process is essentially the same as the one described here. The model dynamics and output functions need to be provided as MATLAB functions via the properties Model.StateFcn and Model.OutputFcn of a nonlinear MPC object. However one should note that the FORCESPRO code generation ignores the jacobian functions that may be provided in Jacobian.StateFcn and Jacobian.OutputFcn, since these will be automatically generated by the automatic differentiation tool CasADi. Moreover, the following requirements on the fields Model.StateFcn and Model.OutputFcn need to be fulfilled for the plugin to work seamlessly:

  • they must be the name of a function file, not an anonymous functions

  • they must be compatible with MATLAB code generation

  • they must follow CasADi coding conventions. Most importantly, the state derivative dxdt has to be built explicitly, as shown below.

dxdt = [expression; expression;]

As a word of caution, the following code snippet will result in an undesired behaviour from CasADI.

dxdt = x;            % Do not write this, CasADI takes it as reference !
dxdt(1,1) = a1*x(1) + a2*x(2) + b1*u(2);
dxdt(2,1) = a3*x(1) + a4*x(2) + b2*u(2);
dxdt(3,1) = x(2)*x(1) + x(4);
dxdt(4,1) = (1/tau)*(-x(4) + u(1));
dxdt(5,1) = x(1) + x(3)*x(6);
dxdt(6,1) = x(2) - 0*x(3);

FORCESPRO calls the model functions from its own objects, which follow an assigment by reference convention, hence the assignement dxdt = x is made by reference. This implies that updating dxdt also changes x, which builds the wrong symbolic dynamics.

If the model contains a parameter, it must be a single vector parameter. In other words, users need to set nlobj.Model.NumberOfParameters = 1 and at run-time write onlinedata.Parameter = value where value is a column vector.

6.3. Generating an NLP solver from an “nlmpc” object

The MATLAB nonlinear MPC API can now be set to use the FORCESPRO code generation. The main difference compared to the existing nonlinear MPC from The MathWorks based on the fmincon solver from the Optimization Toolbox is a code generation step that takes the nonlinear MPC object as argument. This is needed in order to build a mex interface for a FORCESPRO nonlinear solver that is customized to the model provided by the user.

Given an NLMPC object created by the nlmpc command, users can generate an IP or SQP nonlinear solver tailored to their specific problem via the following command:

% nlobj is the output of nlmpc(...)
% options is the output of nlmpcToForcesOptions(...)

[coredata, onlinedata] = nlmpcToForces(nlobj, options);

Two types of nonlinear solvers can be generated via nlmpcToForces: a nonlinear interior-point solver and a sequential quadratic programming solver whose features are covered in details in Sequential quadratic programming algorithm.

The nlmpcToForces API is described in more details in the tables below. The nlmpcToForces command expects an NLMPC object nlobj and a structure options as arguments. Is also has a few limitations as it currently does not support custom cost and constraints. It also requires double precision.

Table 6.1 nlmpcToForces arguments

Input

Description

nlobj

NMPC object constructed by Model Predictive Control Toolbox (see here)

options

Object that provides solver generation options.

The outputs of nlmpcToForces consist of two structures coredata, a structure containing the constant NLMPC information used by nlmpcmoveForces and onlinedata, a structure that allows you to specify online signals such as x, lastMV, ref, MVTarget, md as well as weights or bounds used by nlmpcmoveForces.

In order to provide the solver options to nlmpcToForces, the user needs to run the command nlmpcToForcesOptions. The options structure contains the following fields:

  • SolverName. This is the solver name used by MEX and C files. Its default value is myForcesNLPSolver.

  • SolverType. This option specifies which FORCES nonlinear programming solver to use. Its default value is InteriorPoint. To use the FORCESPRO SQP algorithm set the value to SQP.

  • SkipSolverGeneration. This option indicates whether nlmpcToForces should generate the custom NLP solver. When true, nlmpcToForces will return structures without regenerating the MEX and C files. Its default value is false.

  • Server. This option specifies the FORCES server address for remote solver generation. Its default value is https://forces.embotech.com.

  • PrintLevel. This option specifies the amount of information displayed in the solver log.

    • \(0\): no output will be written

    • \(1\): summary line of each solve

    • \(2\): summary line of each iteration

    Its default value is \(0\).

  • IntegrationNodes. This option specifies the the number of intermediate points between \(t\) and \(t+Ts\) during numerical integration of a continuous time model. Use larger values when the plant is stiff at the price of computational efficiency. Its default value is \(1\). The approach used here is refered to as direct multiple shooting.

  • x0. This option is used to create initial guess of optimal state trajectory at cold start. It must be a column vector of nx-by-1. The typical value should be the initial state of the prediction model. If it is left empty, zeros will be used for cold start. Its default value is [].

  • mv0. This option is used to create initial guess of optimal manipulated variable trajectory at cold start. It must be a column vector of nmv-by-1. The typical value should be the last known control action. If it is left empty, zeros will be used for cold start.

  • Parameter. This option should be specified if the prediction model has a parameter. It must be a column-vector and it can be updated at run-time. Its default value is false.

  • UseMVTarget. This option enables/disables MV reference signal. When equal to true, the MV reference signal is provided via the onlinedata structure. In this case, MV weights should be positive for proper tracking. When equal to false, the MV reference is 0 by default. In this case, MV weights should be zero to avoid unexpected behavior. Default value is false.

  • UseOnlineWeightOV. This option enables/disables online OV weight change. When equal to true, OV weight needs to be provided via onlinedata structure. Its default value is false.

  • UseOnlineWeightMV. This option enables/disables online MV weight change. When equal to true, MV weight needs to be provided via onlinedata structure. Its default value is false.

  • UseOnlineWeightMVRate. This option enables/disables online MVRate weight change. When equal to true, MVRate weight needs to be provided via onlinedata structure. Its default value is false.

  • UseOnlineWeightECR. This field enables/disables online ECR weight change. When equal to true, ECR weight needs to be provided via onlinedata structure. Its default value is false.

  • UseOnlineConstraintStateMax. This option enables/disables online state upper bound change. When equal to true, state upper bound needs to be provided via onlinedata structure. Its default value is false.

  • UseOnlineConstraintStateMin. This field enables/disables online state lower bound change. When equal to true, state lower bound needs to be provided via onlinedata structure. Its default value is false.

  • UseOnlineConstraintOVMax. This field enables/disables online OV upper bound change. When equal to true, OV upper bound needs to be provided via the onlinedata structure. Its default value is false.

  • UseOnlineConstraintOVMin. This option enables/disables online OV lower bound change. When equal to true, OV lower bound needs to be provided via the onlinedata structure. Its default value is false.

  • UseOnlineConstraintMVMax. This field enables/disables online MV upper bound change. When equal to true, MV upper bound needs to be provided via the onlinedata structure. Its default value is false.

  • UseOnlineConstraintMVMin. This field enables/disables online MV lower bound change. When equal to true, MV lower bound needs to be provided via the onlinedata structure. Its default value is false.

  • UseOnlineConstraintMVRateMax. This option enables/disables online MVRate upper bound change. When equal to true, MVRate upper bound needs to be provided via the onlinedata structure. Its default value is false.

  • UseOnlineConstraintMVRateMin. This option enables/disables online MVRate lower bound change. When equal to true, MVRate lower bound needs to be provided via the onlinedata structure. Its default value is false.

The following set of options are specific to the nonlinear interior point solver:

  • IP_MaxIteration. This field specifies the maximum number of iterations in the interior point solver. When the maximum number of iterations is reached (i.e. ExitFlag is \(0\)), the NLP solver aborts calculations and the result should be discarded. Default value is \(200\).

  • IP_Mu0. This field specifies initial barrier parameter. It must be positive and its default value is \(0.1\).

  • IP_BarrierStrategy. This option specifies the strategy used to update the barrier parameter at every iteration of the nonlinear interior point solver. It needs to be either monotone or loqo. logo often leads to faster convergence, while monotone may help convergence for difficult problems. Default value is loqo.

  • IP_LinearSolver. This option sets the linear solver. It must be either normal_eqs, symm_indefinite, or symm_indefinite_fast. With normal_eqs, the KKT system is solved in normal equations form. With symm_indefinite, the KKT system is solved using block-indefinite factorizations. With symm_indefinite_fast, the KKT system is solved in symmetric indefinite form, using regularization and positive definite Cholesky factorizations only. Default value is normal_eqs.

  • IP_EqualityTolerance. This option specifies the tolerance on the nonlinear equality constraints used by the nonlinear interior point solver. It must be positive. Default value is \(10^{-6}\).

  • IP_InequalityTolerance. This field specifies the tolerance on the nonlinear inequality constraints used by the interior-point solver. It needs to be positive and its default value is \(10^{-6}\).

  • IP_StationarityTolerance. This option specifies the tolerance on the stationarity measure used in the nonlinear interior point solver. It needs to be positive and its default value is \(10^{-5}\).

The following set of options are specific to the sequential quadratic programming solver:

  • SQP_MaxIteration. This field specifies the maximum number of iterations used by the inner QP solver. Its default value is \(50\).

  • SQP_MaxQPS. This enables the SQP solver to solve a fixed amount of quadratic approximations at every call to the solver. In general, the more quadratic approximations are solved, the more accurate control performance is achieved. The tradeoff is that the solvetime also increases. The default value is \(1\).

  • SQP_RegHessian. This field stands for the level of regularization of the hessian approximation. Increasing this parameter may help if the SQP solver returns exitflag \(-8\) on your problem. The default value is \(5\cdot{}10^{-9}\).

  • SQP_EqualityTolerance. This option specifies the tolerance on the nonlinear equality constraints. It must be positive and its default value is \(10^{-6}\).

  • SQP_InequalityTolerance. This option specifies the tolerance on the linear inequality constraints. It must be positive and its default value is \(10^{-6}\).

  • SQP_StationarityTolerance. This field specifies the tolerance on stationarity. It must be positive and its default value is \(10^{-5}\).

6.6. Examples

Two examples illustrating the FORCESPRO nlmpc plugin are desribed below. The first example is entirely run by a MATLAB script, whereas the second one is based on the FORCESPRO nlmpc Simulink block.

6.6.1. Controlling a CSTR reactor

In this example we create a nonlinear MPC controller for a CSTR reactor using the MathWorks Nonlinear MPC Plugin. The objective is to control the concentration \(CA\) of reagent \(A\).

You can download the MATLAB code used in this example to try it out for yourself by clicking here along with the dynamics here and the output function here.

Click here for a detailed description of the model. The state of our plant will be denoted by \(x\), while our control input will be denoted by \(u\).

\begin{align*} x_1 &: \text{Reactor temparature}\ (K) \\ x_2 &: \text{Concentration of $A$ in reactor tank}\ \left( \frac{kgmol}{m^3} \right) \\ u_1 &: \text{Jacket coolant temperature} \ (K) \\ u_2 &: \text{Concentration of A in inlet feed stream} \ \left( \frac{kgmol}{m^3} \right) \\ u_3 &: \text{Inlet feed stream temperature} \ (K) \\ \end{align*}

The system dynamics are given by the following first order differential equation

\begin{align*} \dot{x_1} &= (u_3 - x_1) + 0.3 \cdot (u_1 - x_1) + 11.92 \cdot 27944640 \cdot \exp( \tfrac{-5894.14}{x_1} ) \cdot u_2\\ \dot{x_2} &= (u_2 - x_2) - 27944640 \cdot \exp( \tfrac{-5894.14}{x_1}) \cdot u_2\\ \end{align*}

For the purpose of this demonstration the MATLAB function describing the state dynamics will be denoted by exocstrStateFcnCT. Our output \(y\) is simply given by the concentration of \(A\):

\begin{align*} y = x_2 \end{align*}

6.6.1.1. Creating an NLMPC object

The MATLAB function implementing this output will be denoted by exocstrOutputFcn. With the implemented exocstrStateFcnCT and exocstrOutputFcn MATLAB functions at hand we can go ahead create our NLMPC object. The following code-snippet constructs the NLMPC object and specifies our model.

nx = 2;
ny = 1;
nu = 3;
nlobj = nlmpc(nx,ny,'MV',1,'MD',[2 3]);
Ts = 0.5;
nlobj.Ts = Ts;
nlobj.PredictionHorizon = 6;
nlobj.ControlHorizon = [2 2 2];
nlobj.MV.RateMin = -5;
nlobj.MV.RateMax = 5;
nlobj.Model.StateFcn = 'exocstrStateFcnCT';
nlobj.Model.OutputFcn = 'exocstrOutputFcn';

6.6.1.2. Specifying solver options

The followowing specifies the code options specific to FORCESPRO’s MathWorks Nonlinear MPC Plugin:

options = nlmpcToForcesOptions();
options.SolverName = 'CstrSolver';
options.SolverType = 'SQP';
options.IntegrationNodes = 5;
options.SQP_MaxQPS = 5;
options.SQP_MaxIteration = 500;
options.x0 = [311.2639; 8.5698];
options.mv0 = 298.15;

6.6.1.3. Generating the NLP solver

Once we have our NLMPC object and our options we can generate an NLP solver through the nlmpcToForces function:

[coredata, onlinedata] = nlmpcToForces(nlobj,options);

6.6.1.4. Calling the solver

This will generate our NLP solver named CstrSolver. We can call this solver in two different ways:

  • Through the generic nlmpcmoveForces function which comes with the FORCESPRO MathWorks Nonlinear MPC Plugin

  • Or through the generated MEX function nlmpcmove_CstrSolver (the name of the MEX is always “nlmpc_<solver name>”). In general it is advantagous from a performance perspective to use the MEX over the generic nlmpcmoveForces function.

Calling the NLP solver through the generic nlmpcmoveToForces can be done as in the following code-snippet:

onlinedata.md = [10 298.15];
[mv, onlinedata, info] = nlmpcmoveForces(coredata,x,mv,onlinedata);

And the MEX can be called as follows:

[mv, onlinedata, info] = nlmpcmove_CstrSolver(x,mv,onlinedata);

6.6.1.5. Results

The NLP solver generated through the above code-snippets were applied in a simulation for 200 seconds. As can be seen in the plots Figure 6.1, Figure 6.2 and Figure 6.3 the generated solver succeeds in controlling the CSTR reactor with a very fast solvetime while the output stays close to the reference.

../_images/cstr_cost.png

Figure 6.1 Cost as a function of time.

../_images/cstr_solvetime.png

Figure 6.2 Solve time as a function of simulation time.

../_images/cstr_concentration.png

Figure 6.3 Concentration of \(A\) as a function of simulation time.

6.6.3. Deploying the Lane Following Model on Speedgoat

The lane following model in Figure Figure 6.8 can be easily deployed on Speedgoat platforms by means of the code below.

% Choose Speedgoat x86 platform to run FORCESPRO solver
options.ForcesTargetPlatform = 'Speedgoat-x86';
% x0 and u0 are used to create a primal initial guess
options.x0 = x0;
options.mv0 = u0;
% Generate FORCESPRO solver
tm = tic;
[coredata, onlinedata] = nlmpcToForces(nlobj,options);
tBuild = toc(tm);

%%
% Start code generation for Speedgoat x86
mdl = 'LaneFollowingNMPC_Speedgoat_x86';
open_system(mdl);                               % Open Simulink(R) Model
load_system(mdl);                               % Load Simulink(R) Model
rtwbuild(mdl);                                  % Start Code Generation

% Deploy application from the start
tg = slrt;
if(~strcmpi(tg.Application, 'loader'))
    tg.unload();
end
tg.load(mdl);

% Execute application
tg.start();
while(strcmpi(tg.Status, 'running'))
    pause(Ts);
end
scope1 = tg.getscope(1);
scope2 = tg.getscope(2);
scope3 = tg.getscope(3);
../_images/laneFollowSpeedgoatModel.PNG

Figure 6.8 Simulink Real-Time Lane Following model for Speedgoat deployment.

All the files necessary to run this example can be downloaded here.