[ANN] ModelPredictiveControl.jl

Alright, I see how it could be useful.

As I said, since \mathbf{\hat{y}}(k) = \mathbf{h} \big( \mathbf{\hat{x}}_k(k) \big), it is equivalent as bounding the estimated sensor noise \mathbf{\hat{v}}(k) = \mathbf{y}(k) - \mathbf{h} \big( \mathbf{\hat{x}}_k(k) \big) in the MHE. I will move the definition of \mathbf{\hat{V}} and \mathbf{\hat{W}} vectors in the extended help of MovingHorizonEstimator constructor to clarify that point.

I found this great lecture that mentions bounds on \mathbf{\hat{x}}, \mathbf{\hat{v}} and \mathbf{\hat{w}}. I will implement the three types of bounds. It will probably leads to infeasibilities in some case so I will also add soft constraints later.

1 Like

ModelPredictiveControl v0.15.0

The MovingHorizonEstimator types now support hard constraints on the estimated state \mathbf{\hat{x}}, process noise \mathbf{\hat{w}} and sensor noise \mathbf{\hat{v}}. Additionally, MHEs based on LinModels are now treated as quadratic programs (QP) solved with OSQP.jl by default instead of nonlinear problems. From my quick tests, it’s about 500 times faster than v0.14! I never did the exercise for the MHE, and it was surprisingly a long and hard job to construct all the prediction matrices that appear in the quadratic objective and the constraints :laughing:. I was also surprised when I understood that, contrarily to LinMPC, the Hessian matrix of the QP is not constant, because of the time-varying value of the arrival covariance \mathbf{\bar{P}}. This means that solvers need to re-factorize the QP matrices at each discrete time k. This is still drastically cheaper than treating the optimization problem as nonlinear.

The changelog is:

  • added: MHE hard constraints on estimated process and sensor noises
  • added: MHEs based on LinModel are now treated as a quadratic optimization instead of nonlinear
  • added: Plots recipe avoid repeat xlabels for MIMOs
  • added: getinfo method for MovingHorzionEstimator for troubleshooting
  • debug: update MHE arrival covariance with I/Os at arrival
  • improved coverage

The next step is supporting soft constraints for the MHE.


ModelPredicitveControl v0.16.0

I added the support for soft-constraints in the MovingHorizonEstimator. I chose to disable this feature by default for two reasons:

  • most people will only constraints the state estimate. It should not lead to major infeasibility if the bounds are correctly chosen.
  • from my tests, soft-constraining NonLinModel estimations leads to problems that Ipopt has difficulties to solve, because of the weird scaling of the decision variables (no issue with LinModel and OSQP tho). Indeed, with the C weight around 10^5, the solution for \epsilon is typically around 10^{-5}, while the other decision variables follow the scale of the states, generally larger in magnitude. Reducing C and nlp_scaling_max_gradient option of Ipopt can mitigate this.

The changelog is:

  • added: MHE soft constraints with a single slack variable ϵ (a.k.a. equal concern for relaxation, disabled by default thus only hard constraints).
  • changed: more concise MPC logging in REPL.
  • removed: deprecated arguments for setconstraint! and NonLinMPC.
  • doc: hide extended help with details admonition.
  • doc: more concise setconstraint! and ExplicitMPC documentation.

I will try to add an example with the MHE in the manual.


ModelPredictiveControl v.0.19.0

It’s been a while since my last post! Many new features were added since then, with a breaking change: NonLinModel constructor now supposes continuous dynamics by default. Use solver=nothing for discrete-time models.

Here is the changelog since my last update:

  • added: support for NonLinModel with continuous dynamics
  • added: 4th order Runge-Kutta solver with 0 allocation
  • added: support for mutating NonLinModel (in-place) to reduce the allocations
  • added: support for non-diagonal M_Hp, N_Hc and L_Hp weight, allowing terminal costs and LQR tuning
  • doc: simple MHE example on the CSTR (manual)
  • doc: info on terminal weights, continuous NonLinModel and mutating syntax
  • doc: the pendulum example now use the built-in Runge-Kutta solver
  • added some tests with continuous NonLinModel, both mutating and non-mutating

Also, the package is now part of JuliaControl organization! :partying_face: Thanks @baggepinnen for the help!

BTW @odow is there any timeline for when NLopt.jl will support the new NLP syntax ? I finished the migration but I would like at least two FOSS NLP solvers that works (i.e. other than Ipopt.jl) before merging into main.


Let me take a look

1 Like

I will briefly present the content of ModelPredictiveControl.jl in the beautiful Montréal City this summer. See you there!

JuMP-dev 2024 | JuMP for details.


ModelPredictiveControl v.0.21.0

I’m super excited to announce two new major features:

  1. linear model adaptation of controller/estimator at runtime
  2. model linearization at non-equilibrium points

These features allow adaptive MPC based on online parameter estimation (e.g. recursive system identification). Furthermore, successive linearization MPC can now be implemented with minimal effort. It’s a great feature since e.g. MATLAB also supports it, but the analytical expressions of the Jacobians are required (hidden in the “Successive Linearizer” block). This is not the case here with the AD-based linearize function. Also, similar functionalities are available for the MovingHorizonEstimator.

A massive and long refactoring of the operating points and the mathematical notation was required to support linearization at non-equilibrium points. It was totally worth it since the successive linearization MPC on the pendulum is about 125 times faster than the NMPC counterpart, with similar closed-loop performances! It’s a great application of Julia’s AD tools and its first class linear algebra library!

All the keyword arguments related to initial values e.g. σP0, x0 and x̂0 now require an underscore e.g. σP_0, x_0, x̂_0 (to clearly differentiate from deviation vectors around operating points)

  • Added: setmodel! for runtime model adaptation of controller/estimator based on LinModel
  • Added: linearize and setop! now support non-equilibrium points
  • Added: successive linearization MPC with the new setmodel! and linearize functions
  • Added: successive linearization MHE with the new setmodel! and linearize functions
  • Added: linearize! method for in-place model linearization (to reduce allocations)
  • Added: 6 args. LinModel constructor now support scalars (similarly to ss function)
  • Added: ExtendedKalmanFilter now compute the Jacobians in-place (to reduce allocations)
  • Changed: struct state data x and state estimate renamed to x0 and x̂0
  • Debug: ExplicitMPC with non-Float64 now works
  • Debug: accept integers in linearize arguments
  • Debug: call empty! on JuMP.Model to support re-construction of MPC instances
  • Doc: new setmodel!, setop! and linearize function documentation
  • Doc: example of model adaptation with successive linearization on the pendulum (very efficient!)

other changes, since my last post on discourse:

  • Added: print info on controller and estimator constraint softening (slack var. ϵ)
  • Added: custom estimator for the approximation of the arrival covariance in the MHE
  • Added: improve performance of LinMPC and MHE with new JuMP batch update methods
  • Various allocation reductions
  • Changed: MHE now default to an UnscentedKalmanFilter for the arrival covariance estimation (instead of an EKF)
  • Debug: MHE with NonLinModel update covariance with the correct state estimate

One MPC problem that I’ve never seen anyone really address…

  • MPC algorithms tend to assume that each optimization takes zero time, in other words: that the computed optimal solution u*[t_k] for u[t_k] can be injected into the system at time t = t_k.
  • In practice, computation takes some time, so the computed value u*[t_k] is actually injected at a delayed time t = t_k+tau, i.e., u[t_k+tau] = u*[t_k], while for t_k-1 + tau < t < t_k+tau, u[t] = u*[t_k-1]

In other words: in real life, MPC adds a time delay to the system – an this time delay may actually vary – depending on constraints in the system, current state, etc. This may or may not be a problem.

One way around this is to assume a maximum time delay and then include this time delay in the model so that the computation of u*[t_k] already is aware of this time delay.

Anyways, I don’t know whether your MPC package can address such problems?

1 Like

The approach of modeling the computational delay as an input delay in the plant is of course available, even if the delay is a fraction of a sample (using Thiran approximation). When using linear MPC, neglecting the computational delay is not always such a bad idea, it takes a couple of microseconds to solve even reasonably large problems (in state dimension and prediction horizon). For computationally more demanding problems, real-time iteration (RTI) is often worthwhile, this is something that isn’t yet supported, but would be good to get in there eventually. With RTI, you don’t bother iterating each problem until convergence, and instead favor to get new measurement data at a higher rate such that later iterations are made using the latest available data. The reasoning being that performing optimizer iterations with “old data” is wasteful if new data is available.


I agree on linear MPC, although it may depend on whether you allow for hard constraints on inputs, and more importantly on outputs.

Already 20+ years ago, Biegler and co-workers solved linear MPC with 200 inputs and 800 outputs (?) in … 10 or 100 ms or so; they considered only input constraints.

So if one has n states and a horizon of N, this gives in the order of n\cdot N unknowns (+ inputs etc., etc.). It then probably is important to take advantage of the structure of the resulting problem. Without constraint, this can be solved by handling the n\cdot N \times n\cdot N matrix involved (ok: really larger, by including inputs, outputs, etc.). Alternatively, one can instead solve this using the Riccati-equation formulation, which reduces the complexity to handling N problems of size n\times n. At the outset, this only works for unconstrained problems. But I seem to recall that in a PhD thesis (or licenciate thesis?) from Uppsala (?) some 15-20 years ago, someone (don’t recall the name) reformulated linear quadratic problems with input constraints into a Riccati-like form. [This could also be interesting for receding horizon state estimation.]

With output constraints, QP type algorithms are needed, which may slow the solution down considerably if the problem in certain operating points is infeasible or near infeasible. Output constraints are important in some applications, e.g., electricity generation where there are (semi-) hard constraints on grid frequency, etc.

But if the system is nonlinear, or worse: distributed (spatial distribution, population balance, etc.), the problem I raise is more interesting. Of course, for systems that are slow to solve, in practical implementations, one may resort to surrogate models, etc.

1 Like

ModelPredictiveControl v0.22.0

A quick update on the last improvements. I first updated the internal code to the new NLP syntax of JuMP.jl. It is a slightly breaking change for some custom optimizers. There is only NLopt.jl that does not support it yet. I tested Ipopt.jl, MadNLP.jl and KNITRO.jl and they work well.

Also, following the model adaptation feature, it is now possible to modify the weight matrices in the objection function at runtime using setmodel!. Same thing is possible with the covariance matrices of the Kalman Filter and moving horizon estimator.

Release notes:


Migration to the new nonlinear programming syntax of JuMP. Some optimizers may not support it yet.

  • added: online modification of objective function weights and covariance matrices with setmodel!
  • added: accept indices and ranges in plot keyword arguments
  • changed: migrate the NLP code to the new syntax
  • doc: added plot result previews with static images
  • new tests for new plot recipes and weights/covariances modification

and, the changelog since my last post:

  • added: non-Unicode alternative keyword arguments in public functions
  • added: modify plots Y-axis labels with setname!
  • added: show sim! % progress and ETA in VS Code status bar
  • added: plot estimated state constraints x̂min/x̂max for moving horizon estimator
  • added: plot recipe documentation
  • debug: debug: plots Ru instead of Ry for u setpoint (recipe)
  • various doc improvements

The package will be registered soon.