Added: adaptation of controller/estimator based on `LinModel`

[franckgaga]

@baggepinnen New model adaptation feature, related to #48 request.

I don't ask to review everything, there is way too much modifications! Here's a quick summary:

For the support of adaptive linear MPC and state estimators, I completely refactored the handling of the operating points (and also my notation in the doc/struct). The old implementation with uop, yop and dop only was too simplistic. It works well when the model comes from system identification at an equilibrium point, but it is not sufficient for linearizing a model at any points, including non-equilibrium. The xop and fop constants are now needed. I also tried to combine xop and fop into a single constant for simplification, but it's impossible: changing the model at runtime implies a rescaling of the current state estimate with the old and new xop value.

See https://courses.engr.illinois.edu/ece486/fa2021/documentation/lectures/slides/Lecture25B_Linearization.pdf (slide 10) for the procedure on continuous dynamics at equilibrium (discrete case is similar).

I added some details in setop! docstrings:

linearize function:

and setmodel! function:

I added an example of model adaptation based on successive linearization on the pendulum in the manual. The closed-loop response is very similar to the nonlinear MPC, and the computations are about 125 times faster, impressive! Note that MATLAB also support it (https://www.mathworks.com/help/mpc/ug/adaptive-mpc-control-of-nonlinear-chemical-reactor-using-successive-linearization.html) but inside the "Successive Linearizer" block is in fact the analytical derivative of the CSTR model. It's way more simple with ForwardDiff.jacobian!

I also added simple tests for setmodel! on the estimator/controller that support it.

[franckgaga]

P.S. A new constant vector $\mathbf{B}$ / $\mathbf{b_\hat{x}}$ was also needed in the predictions :

similar vectors is needed for the MHE. I did not have the time to implement it. For now setmodel!(::MovingHorizonEstimator, ::LinModel) return an error message that it's not implemented yet.

[franckgaga]

I enabled push_preview in deloy_doc. Here's the peview for the SLMPC on the pendulum : https://juliacontrol.github.io/ModelPredictiveControl.jl/previews/PR54/manual/nonlinmpc/#Adapting-the-Model-via-Successive-Linearization

[franckgaga]

I'll just merge into main and you can comment here for advices/corrections. They will be included in the next version.

[baggepinnen]

Hey! Sorry it takes me a while to get to things on github at the moment, I'm on parental leave so I have a new full-time job. :)

I was initially a bit confused about what exactly "adaptation" referred to in this context, this is not adaptation as in adaptive control (example), but a heuristic way to achieve nonlinear MPC by solving linear MPC problems and using a new linearization at each time step? In any case, it sounds like a feature worth having!

Having figured that out, it looks like the new feature is reached by

linmodel = linearize(nonlinmodel; u, x=x̂[1:2])
setmodel!(mpc, linmodel)

i.e., the user may call linearize, or figure out a new linear model in any other way they see fit, and then pass this to mpc to be used in the next solve. That seems like a nice interface, it allows "custom" adaptation that may be different from calling linearize on a nonlinear model :+1:

Questions:

• What does the $p$ here do?

[franckgaga]

No worry, in the end open source development is voluntary work :). My working week is 4 days and I don't have kids, so more time available.

The "p" should be in subscript. The correct additive constant to update the deviation vector of the state $\mathbf{x0}$ should be $\mathbf{f{op} - x_{op}}$ (see eq. (11) here). It's corrected in the latest documentation:

A mistake in $\mathbf{B}$ was also corrected in the prediction matrices:

FYI, MATLAB call the two operating points $\bar{\Delta x} = \mathbf{f{op} - x{op}}$ and $\bar{x} = \mathbf{x{op}}$. I think it's clearer with my notation, since we have directly $\mathbf{f{op} = f ( x{op}, u{op}, d_{op} )}$

There is two approaches for adaptation of a linear plant model. To quote MATLAB doc (the link with the CSTR example above):

If a linear plant model can be obtained at run time, you should use Adaptive MPC Controller block to achieve nonlinear control. There are two typical ways to obtain a linear plant model online:

(1) Use successive linearization as shown in this example. Use this approach when a nonlinear plant model is available and can be linearized at run time.

(2) Use online estimation to identify a linear model when loop is closed. See Adaptive MPC Control of Nonlinear Chemical Reactor Using Online Model Estimation for more details. Use this approach when linear plant model cannot be obtained from either an LPV system or successive linearization.

I applied the first approach in the manual since it was straightforward, and it shows how to correctly call setmodel! in closed-loop.

The setmodel! can be used for any adaptation scheme of LinModel, either a linearization or an online parameter estimation algorithm like a Kalman Filter augmented with the time-varying parameters, or e.g. recursive ARX.

P.S. In your adaptive MPC example, I don't understand why and how the OSQP solver (quadratic programming) is used to solve the optimization problem if the plant model is nonlinear (it should be a nonlinear program). You also perform succesive linearization ?

[baggepinnen]

In your adaptive MPC example, I don't understand why and how the OSQP solver (quadratic programming) is used to solve the optimization problem if the plant model is nonlinear (it should be a nonlinear program). You also perform succesive linearization ?

I use SQP (sequential quadratic programming), so I linearize around the trajectory, solve QP, linearize again and iterate

[franckgaga]

Oh I see, you implemented it yourself ? Impressive! A FOSS SQP method is lacking in the current JuMP ecosystem.

If I understand correctly, my sequential linearization MPC is equivalent to your solution with one iteration.

[baggepinnen]

Oh I see, you implemented it yourself ? Impressive! A FOSS SQP method is lacking in the current JuMP ecosystem.

Yeah, but I also wish there was a robust solution available in the ecosystem, my solver is not battle tested and likely not as robust as I'd like.

If I understand correctly, my sequential linearization MPC is equivalent to your solution with one iteration.

Not quite, you linearize in a single point before each solve, and before the next solve linearize in a new point. SQP instead linearizes around the trajectory of the previous solution, i.e., before each solve, the model is linearized in N different points for a prediction horizon N. This scheme can be iterated until convergence, since even if you solve the QP to convergence, you haven't solved the original problem to convergence due to the linear time-varying approximation made in the trajectory linearizaiton. In practice, one may choose to not perform the outer iterations to convergence, and instead opt for a faster sample rate of the controller. The reasoning here being that all but the first iterations are solved using "old data", and it would be better to take a new measurement and solve the next iteration with an updated initial state.