Configure Optimization Solver for Nonlinear MPC

By default, nonlinear MPC controllers solve a nonlinear programming problem using the `fmincon` function with the SQP algorithm, which requires Optimization Toolbox™ software. Alternatively, you can specify a continuation/generalized minimum residual (C/GMRES) solver or your own custom nonlinear solver, neither of which require Optimization Toolbox software.

Solver Decision Variables

For nonlinear MPC controllers at time tk, the nonlinear optimization problem uses the following decision variables:

• Predicted state values from time tk+1 to tk+p. These values correspond to rows 2 through p+1 of the `X` input argument of your cost and constraint functions, where p is the prediction horizon.

• Predicted manipulated variables from time tk to tk+p-1. These values correspond to the manipulated variable columns in rows 1 through p of the `U` input argument of your cost and constraint functions.

Therefore, the number of decision variables NZ is equal to p(Nx+Nmv) + 1, Nx is the number of states, Nmv is the number of manipulated variables, and the +1 accounts for the global slack variable.

The problem formulation used here, in which all the intermediate state variables are part of the vector of decision variables, is called a sparse formulation. It is also known as a simultaneous or multiple-shooting formulation. In this formulation, equality constraints of the type x(k+1) = f(x(k),u(k)) must be enforced. This formulation results in larger problems with many equality constraints, but it leads to a sparse, regular, structure. Such structure yields better numerical stability and conditioning for large horizons and for unstable or nonsmooth plants. It can also result in a much more efficient optimization when you use a (possibly parallel) solver that can exploit sparsity. The sparse formulation is the default used by Model Predictive Control Toolbox™ for nonlinear MPC problems (when `fmincon` is used as a solver).

The sparse formulation differs from the dense formulation, in which equality constraints such as x(k+1) = f(x(k),u(k)) are solved using numerical simulation, and the intermediate state variables are eliminated by substitution into the cost and constraint functions, thereby leaving just the manipulated variables (plus a slack variable) as decision variables of the optimization problem.

The dense formulation is also referred to as a sequential or single-shooting formulation, and is used by Model Predictive Control Toolbox for the C/GMRES solver. A dense formulation is also the default for linear MPC problems, as described in Optimization Problem.

Specify Initial Guesses

A properly configured standard linear MPC optimization problem has a unique solution. However, nonlinear MPC optimization problems often allow multiple solutions (local minima), and finding a solution can be difficult for the solver. In such cases, it is important to provide a good starting point near the global optimum.

During closed-loop simulations, it is best practice to warm start your nonlinear solver. To do so, use the predicted state and manipulated variable trajectories from the previous control interval as the initial guesses for the current control interval. In Simulink®, the Nonlinear MPC Controller block is configured to use these trajectories as initial guesses by default. To use these trajectories as initial guesses at the command line:

1. Return the `opt` output argument when calling `nlmpcmove`. This `nlmpcmoveopt` object contains any run-time options you specified in the previous call to `nlmpcmove`. It also includes the initial guesses for the state (`opt.X0`) and manipulated variable (`opt.MV0`) trajectories, and the global slack variable (`opt.Slack0`).

2. Pass this object in as the `options` input argument to `nlmpcmove` for the next control interval.

These command-line simulation steps are best practices, even if you do not specify any other run-time options.

Configure `fmincon` Options

By default, nonlinear MPC controllers optimize their control move using the `fmincon` function from the Optimization Toolbox. When you first create your controller, the `Optimization.SolverOptions` property of the `nlmpc` object contains the standard `fmincon` options with the following nondefault settings:

• Use the `SQP` algorithm (```SolverOptions.Algorithm = 'sqp'```)

• Use objective function gradients (```SolverOptions.SpecifyObjectiveGradient = 'true'```)

• Use constraint gradients (```SolverOptions.SpecifyConstraintGradient = 'true'```)

• Do not display optimization messages to the command window (`SolverOptions.Display = 'none'`)

These nondefault options typically improve the performance of the nonlinear MPC controller.

You can modify the solver options for your application. For example, to specify the maximum number of solver iterations for your application, set `SolverOptions.MaxIter`. For more information on the available solver options, see `fmincon` (Optimization Toolbox).

In general, you should not modify the `SpecifyObjectiveGradient` and `SpecifyConstraintGradient` solver options, since doing so can significantly affect controller performance. For example, the constraint gradient matrices are sparse, and setting `SpecifyConstraintGradient` to false would cause the solver to calculate gradients that are known to be zero.

Configure C/GMRES Solver Options

The C/GMRES method is another built-in option for optimization. The computational efficiency and convergence properties of the C/GMRES method are particularly useful when applied to large scale systems and systems that exhibit nonlinear and non-convex constraints.

To use the C/GMRES method to solve multistage nonlinear MPC problems, specify the `Optimization.Solver` property of your multistage nonlinear MPC object as `"cgmres"`. When you do so, the software stores a default C/GMRES options object in the `Optimization.SolverOptions` property of your multistage MPC object.

The C/GMRES options object has the following properties, which you can modify using dot notation.

• `BarrierParameter` — Parameter barrier value. This is a positive scalar representing the strength of the penalty term added to the constraints in the barrier function method for nonlinear model predictive control. The penalty term is used to penalize lower and upped bounds constraint violations. The default value is `0.1`.

• `Display` — Level of display. This is a string or character vector specified as either `"iter"` or `"none"`. Use `"iter"` to display solver information at every iteration, and use `"none"` to suppress any solver display. The default value is `"none"`.

• `FiniteDifferenceStepSize` — Step size used for finite differences. This is a positive scalar representing the step size that the forward finite difference scheme uses to approximate the Jacobian in the optimization problem. A larger step size can result in faster convergence but may also lead to numerical instability. The default value is `1e-8`.

• `MaxIterations` — Maximum number of iterations allowed. This is a positive integer representing the maximum number of inner-loop iterations allowed for the solver. The default value is `100`.

• `Restart` — Number of outer-loop iterations. This is a nonnegative integer representing the maximum number of outer-loop iterations for the solver. The default value is `10`.

• `StabilizationParameter` — Stabilization parameter. This is a positive scalar representing a forcing term that you can use to control the step size and improve the numerical stability of the solver.

In inexact Newton methods, the Newton direction is typically approximated using an iterative linear solver, such as the conjugate gradient method or the GMRES method. However, due to computational limitations or other factors, the iterative solver may not find an exact solution to the linear system within a specified tolerance. The forcing term is introduced to compensate for the inexactness of the Newton direction. It adjusts the step size taken in each iteration to ensure convergence, even when the Newton direction is not accurate enough. The forcing term can be viewed as a trust region radius that limits the step size based on the accuracy of the Newton direction. Therefore a decrease of the stabilization parameter normally leads to a decrease of the step size.

If the `StabilizationParameter` is too large, it may cause overcorrection and oscillations, leading to slow convergence or even divergence. On the other hand, if the `StabilizationParameter` is too small, it may result in slow convergence as the step size is overly conservative. To enhance algorithm convergence, it is advisable to select a `StabilizationParameter` value of 1/Ts, where Ts represents the model sampling time. This parameter should not exceed 2/Ts.

The default value is `1e3`.

• `TerminationTolerance` — Termination tolerance. This is a positive scalar representing the termination tolerance for the solver. The solver stops when the relative error between two consecutive iterations is less than the termination tolerance. The default value is `1e-6`.

For an example on how to use the C/GMRES solver, see Control Robot Manipulator Using C/GMRES Solver and Solve Fuel-Optimal Navigation Problem Using C/GMRES.

Specify Custom Solver

As an alternative to the `fmincon` function and C/GMRES method, you can specify your own custom nonlinear solver. To do so, create a custom wrapper function that converts the interface of your solver function to match the interface expected by the nonlinear MPC controller. Your custom function must be a MATLAB® script or MAT-file on the MATLAB path. For an example that shows a template custom solver wrapper function, see Optimizing Tuberculosis Treatment Using Nonlinear MPC with Custom Solver.

You can use the Nonlinear Programming solver developed by Embotech AG to simulate and generate code for nonlinear MPC controllers. For more information, see Implement MPC Controllers Using Embotech FORCESPRO Solvers.

To configure your `nlmpc` object to use your custom solver wrapper function, set its `Optimization.CustomSolverFcn` property in one of the following ways:

• Name of a function in the current working folder or on the MATLAB path, specified as a string or character vector

`Optimization.CustomSolverFcn = "myNLPSolver";`
• Handle to a function in the current working folder or on the MATLAB path

`Optimization.CustomSolverFcn = @myNLPSolver;`

Your custom solver wrapper function must have the signature:

`function [zopt,cost,flag] = myNLPSolver(FUN,z0,A,B,Aeq,Beq,LB,UB,NLCON)`

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

• NZ is the number of decision variables.

• Mcineq is the number of linear inequality constraints.

• Mceq is the number of linear equality constraints.

• Ncineq is the number of nonlinear inequality constraints.

• Nceq is the number of nonlinear equality constraints.

ArgumentInput/OutputDescription
`FUN`Input

Nonlinear cost function to minimize, specified as a handle to a function with the signature:

`[F,G] = FUN(z)`

and arguments:

• `z` — Decision variables, specified as a vector of length NZ.

• `F` — Cost, returned as a scalar.

• `G` — Cost function gradients with respect to the decision variables, returned as a column vector of length NZ, where $\text{G}\left(i\right)=\partial F/\partial z\left(i\right)$.

`z0`InputInitial guesses for decision variable values, specified as a vector of length NZ
`A`InputLinear inequality constraint array, specified as an Mcineq-by-NZ array. Together, `A` and `B` define constraints of the form $\text{A}\cdot z\le \text{B}$.
`B`InputLinear inequality constraint vector, specified as a column vector of length Mcineq. Together, `A` and `B` define constraints of the form $\text{A}\cdot z\le \text{B}$.
`Aeq`InputLinear equality constraint array, specified as an Mceq-by-NZ array. Together, `Aeq` and `Beq` define constraints of the form $\text{Aeq}\cdot z=\text{Beq}$.
`Beq`InputLinear equality constraint vector, specified as a column vector of length Mceq. Together, `Aeq` and `Beq` define constraints of the form $\text{Aeq}\cdot z=\text{Beq}$.
`LB`InputLower bounds for decision variables, specified as a column vector of length NZ, where $\text{LB}\left(i\right)\le z\left(i\right)$.
`UB`InputUpper bounds for decision variables, specified as a column vector of length NZ, where $z\left(i\right)\le \text{UB}\left(i\right)$.
`NLCON`Input

Nonlinear constraint function, specified as a handle to a function with the signature:

`[cineq,c,Gineq,Geq] = NLCON(z)`

and arguments:

• `z` — Decision variables, specified as a vector of length NZ.

• `cineq` — Nonlinear inequality constraint values, returned as a column vector of length Ncineq.

• `ceq` — Nonlinear equality constraint values, returned as a column vector of length Nceq.

• `Gineq` — Nonlinear inequality constraint gradients with respect to the decision variables, returned as an NZ-by-Ncineq array, where $\text{Gineq}\left(i,j\right)=\partial \text{cineq}\left(j\right)/\partial z\left(i\right)$.

• `Geq` — Nonlinear equality constraint gradients with respect to the decision variables, returned as an NZ-by-Nceq array, where $\text{Geq}\left(i,j\right)=\partial \text{ceq}\left(j\right)/\partial z\left(i\right)$.

`zopt`OutputOptimal decision variable values, returned as a vector of length NZ.
`cost`OutputOptimal cost, returned as a scalar.
`flag`Output

Exit flag, returned as one of the following:

• `1` — Optimization successful (first order optimization conditions satisfied)

• `0` — Suboptimal solution (maximum number of iterations reached)

• Negative value — Optimization failed

When you implement your custom solver function, it is best practice to have your solver use the cost and constraint gradient information provided by the nonlinear MPC controller.

If you are unable to obtain a solution using your custom solver, try to identify a special condition for which you know the solution, and start the solver at this condition. If the solver diverges from this initial guess:

• Check the validity of the state and output functions in your prediction model.

• If you are using a custom cost function, make sure it is correct.

• If you are using the standard MPC cost function, verify the controller tuning weights.

• Make sure that all constraints are feasible at the initial guess.

• If you are providing custom Jacobian functions, validate your Jacobians using `validateFcns`.