[ANN] ModelPredictiveControl.jl

[langestefan]

Thanks. I guess I missed this sentence.

all the keyword arguments above but with a first capital letter, except for the terminal constraints, e.g. Ymax or C_Δumin: for time-varying constraints (see Extended Help)

A nicer interface would be something akin to setconstraint!(variable, index, lb, ub).

Or just a mini DSL that can formulate things like x + y \geq a, or x + y = a.

For example for my problem I also need to be able to define constraints of the form:

\begin{align} x_t &\geq a_t - z_t\\ y_t &\geq z_t - a_t \\ x_t, y_t &\geq 0 \end{align}

Is that possible?

For context, I am working on the ‘real’ version of Economic MPC for residential HVAC · DyadControlSystems by Fredrik

[franckgaga]

What do you mean by “variable” ? Decision variables, outputs, states ? And “index” ? The discrete time step ?

What is x_t, y_t and z_t here ? States, manipulated inputs or outputs ?

[langestefan]

Was trying to avoid some jargon but the form is:

\begin{align} \underline{u}_t &\geq T^\mathrm{set}_t - T^i_t\\ \overline{u}_t &\geq T^i_t - T^\mathrm{set}_t\\ \underline{u}_t, \overline{u}_t &\geq 0 \end{align}

Where \underline{u}_t and \overline{u}_t are the lower and upper comfort penalties. T^i_t is a state variable, and T^\mathrm{set}_t is a setpoint. Additional slack can be added to create a ‘comfort band’ that a user can define.

I also have constraints of the form:

\begin{align} Q_t^\mathrm{HP} &\leq b_0 + b_1 T^a_t + b_2 T^s_t \end{align}

Where T^s_t is a state variable and Q_t^\mathrm{HP} an input variable. Output thermal power limits typically depend on ambient and water temperatures.

[franckgaga]

No, this structure is currently not supported by LinMPC object. It may be supported by the LinearMPC.jl package however, see Constraints and add_constraint!. The documentation is scarce but you are interested in Ax, Ar and ks arguments (x for state, r for setpoint and ks for the discrete time steps, starting at 1 for the present time step).

My main priority is to provide a user interface that is easy to use, meaning it is more restrictive in what can be accomplished in terms of objectives and constraints. Currently, if you have complex constraint or objective structures, there is always the NonLinMPC object that allows custom nonlinear constraints (gc and nc kwargs) and objectives (Ewt and JE kwargs). But it will obviously not be the most computationally efficient algorithm if these structures are truly linear and quadratic. And code generation is not supported for NonLinMPC.

That being said, my second priority is to be as flexible as possible. A feature that I could introduce is the possibility to add custom linear inequality constraints in the form of (with \mathbf{\tilde{Z}} being the decision variables augmented with the scalar slack variable \epsilon, if activated):

\mathbf{A_c} \mathbf{\tilde{Z}} \le \mathbf{b_c}

But there are two issues with this:

• It needs a lot of work from the user (as in, you) to compute the \mathbf{A_c} matrix and \mathbf{b_c} vector.
• For your specific structure, the \mathbf{b_c} vector will need to be re-computed at each control period k by the user (since the setpoints may change). This feature is not supported by code generation (i.e. there is no C code equivalent of calling setconstraint! online).

A second option would be to introduce an API similiar to LinearMPC’s add_constraint! in ModelPredictiveControl.jl, if it fits your need. And it would obviously support code generation.

[langestefan]

@franckgaga

My main priority is to provide a user interface that is easy to use, meaning it is more restrictive in what can be accomplished in terms of objectives and constraints.

The expressiveness of the linear MPC framework is already very limited and the current syntax limits it further. As long as my expression is mathematically valid within the linear MPC framework I think I should be able to define it.

@franckgaga

Currently, if you have complex constraint or objective structures, there is always the NonLinMPC object that allows custom nonlinear constraints (gc and nc kwargs) and objectives (Ewt and JE kwargs). But it will obviously not be the most computationally efficient algorithm if these structures are truly linear and quadratic. And code generation is not supported for NonLinMPC.

Unfortunately for commercial embedded systems every kilobyte matters so C codegen is a must. I think it’s already a modern miracle we can do MPC with a <200 kB solver (thanks Boyd!).

@franckgaga

This feature is not supported by code generation (i.e. there is no C code equivalent of calling setconstraint! online).

Yeah this is very tricky but I think a requirement for practical applications. Between mild and cold weather with high supply temperatures maximum output power can drop from 5 to 3 kW.

Maybe you know this already, but cvxpygen will generate a function cpg_update_<param>(idx, val) to update each parameter in your problem. For that to work you would first need to know where the parameters are and there are some rules as to how parameters can enter the problem (DPP compliance).

Now for LinearMPC.jl, you could say everything that is a constant in my problem is a parameter and thus could be updated and included in the codegen. Having both an expressive API and efficient codegen is even more tricky and you’re halveway to a DSL already.

Really thinking out loud here so feel free to shoot this down, but why not use JuMP as the main problem creation interface instead? There’s a DSL, and there is MOI.Parameter.

ps. It goes without saying but I appreciate your work on these packages.

[franckgaga]

ModelPredictiveControl.jl v1.16.0

An update to annonce the support of custom linear inequality constrains in LinMPC and exporting to C code!

Thanks for the feedback @langestefan

edit: oopsie, I forgot to merge the branch with the C codegen support . It will be available in v1.16.1 before the end of the day.

[franckgaga]

ModelPredictiveControl.jl v2.0.0

This is major release to introduce the new built-in OrthogonalCollocation transcription. There are also two breaking changes:

1. The oracle keyword argument is removed from NonLinMPC and MovingHorizonEstimator constructor. The nonlinear constraints now always rely on the VectorNonlinearOracle. The set is now supported by most solvers incl. Ipopt.jl, MadNLP.jl, UnoSolver.jl and KNITRO.jl, plus the legacy splatting syntax is slow and no longer fully tested, so deleting this codebase is worthwhile. It was also easy to forget setting oracle=true with custom optimizer.
2. The slack variable \epsilon is added in the signature of the custom economic function J_E of NonLinMPC. It can be useful for non-quadratic weighting. The new signature is JE(Ue, Ŷe, D̂e, p, ϵ) -> cost.

About the OrthogonalCollocation, the number of decision variables is even larger than MultipleShooting, but the problem is sparser and tailored for highly stiff nonlinear systems. It supports both :gaussradau (default and L-stable) and :gausslegendre (A-stable) quadratures. Multi-threading is also supported with the f_threads and h_threads keyword argument. Note that it currently supports only piecewise constant manipulated inputs \mathbf{u}.

It was slightly longer than expected to implement since I did not find a single reference that explained all the important information to code it from scratch. For example, most references skip the part about pre-computing a differentiation matrix and how to use it on the collocation points (incl. do-mpc doc and Biegler, 2010). Other references skip the part about the continuity constraint. If you are curious, I tried to summarize the procedure as concisely as possible in the internals, see:

• the related init_orthocolloc method, and
• the related con_nonlinprogeq! method.

The release will be registered soon. I will work on adding piecewise linear inputs with the h=1 argument.

The change log since my last post is:

• BREAKING: removed oracle keyword from NonLinMPC and MovingHorizonEstimator
• BREAKING: added slack ϵ in the JE function signature of NonLinMPC
• added: OrthogonalCollocation with gaussradau and gausslegendre schemes
• added: custom linear constraint vector \mathbf{W} in getinfo dictionary
• debug: avoid duplicate precompile macros
• debug: InternalModel now work with all CollocationMethods
• removed: specialized obj_nonlinprog! for LinModel in NonLinMPC (slower)
• removed: deprecated preparestate! method
• doc: concise documentation of OrthogonalCollocation internals

[franckgaga]

ModelPredictiveControl.jl v2.1.0

The last release include some bugfixes related to race condition with multi-threading on the collocation methods. It it also introduces the support of piecewise linear manipulated input \mathbf{u} for OrthogonalCollocation with h=1.

I’ve run some local benchmark to compare OrthogonalCollocation to the other transcription methods on the inverted pendulum. It is indeed slightly more expensive than everything else, but nothing excessive. But it’s interesting to point out that Ipopt.jl seems to struggle with it compared to UnoSolver.jl + filtersqp, since it’s about 3 times slower. Once more, many thanks @cvanaret !

Benchmarks

Case study	v2.1.0
Ipopt – SingleShooting – Hessian	0.194 s
Ipopt – MultipleShooting – Hessian	0.431 s
Ipopt – TrapezoidalCollocation – Hessian	0.356 s
Ipopt – OrthogonalCollocation – Hessian	2.049 s
UnoSolver – MultipleShooting – Hessian	0.152 s
UnoSolver – OrthogonalCollocation – Hessian	0.615 s

[franckgaga]

ModelPredictiveControl v2.4.0

It’s been a while since my last update! First, if you miss my post on LinkedIn, I updated the arXiv paper with the new features and up-to-date benchmark results. The NonLinMPC controller is now up to 6× faster than an equivalent MATLAB implementation. The LinMPC performances also slightly improved both on OSQP.jl and DAQP.jl solvers.

The v2.4 release adds the support of custom nonlinear inequality constraint for the MovingHorzionEstimator. It was more complex than the equivalent feature for NonLinMPC for two reasons. First, the data windows are growing at the beginning when k < H_e. Second, the windows are inherently “misaligned” in terms of time steps. For example, the current manipulated input \mathbf{u}(k) is unavailable with the current form (direct=true), but all the measurements are normally available. I chose to introduce extended vectors (denoted with the subscript \mathbf{e}) to artificially align all the time series in the argument of the custom constraint function, and substituting the missing values with NaN. The arguments are detailed in the extended help.

Also, the v2.3 release introduces new information in the getinfo dictionnary: the number of colors of the sparsity patterns for all the Jacobians and Hessians in it. It can be useful to tweak the ordering algorithm in order to improve performances at sparse differentiation.

Here’s the changelog since my last post:

• added : custom nonlinear inequality constraints for MovingHorizonEstimator
• added : reduce allocations in MovingHorizonEstimator and initstate!
• added: number of colors in Jacobians and Hessians sparsity patterns of getinfo dictionnary
• added: current estimate x̂ and last input lastu in getinfo dictionnary for MPCs
• added: stochastic state defects as linear equality constraints for non-SingleShooting
• changed: exclude Inf bounds in the inequality constraint fields of getinfo
• test: add Aqua.jl to the test suite
• doc: updating MTK example to v11 in manual
• doc: various improvements

[ufechner7]

Thank you so much for the great work!

Question: Why is explicit support for DAE’s needed?

We are mostly controlling kites, and we have numeric models and MTK based models.

Can we already use your package to control our models, or do we need the aforementioned feature implemented?

[franckgaga]

Thanks!

No, except if the DAE is of index 1, which is in fact almost an ODE. Another way is by modifying the model assumptions and equations to eliminate the fast dynamics and the resulting algebraic equations by substitutions, when possible. For example, we can frequently assume that the fast reactions in a bioreactors are at steady-state (thus instantaneous) and eliminate the resulting implicit algebraic equation by direct substitution. The solving of the optimization problem will be just easier, so it’s always a good idea to do these simplifications when possible.

Else, yes I need to implement the explicit support of DAEs. This is not easy since state estimation for DAEs is not as well developed as for ODEs. The most “straightforward” approach is adding the support of orthogonal collocation to the MovingHorizonEstimator, but it still needs a lot of work since everything is hard-coded as a single shooting transcription right now. And, once more, this is more complex than the equivalent feature for NonLinMPC, mainly because the number of decision variables grows at the beginning (because of the growing windows).

[Bart_van_de_Lint]

@ufechner7 our MTK model is an ODE

[Bart_van_de_Lint]

Thanks for your great work. Very high quality code and nice paper.

[franckgaga]

Many thanks!
Yes, I try to write code as readable, maintainable and idiomatic to Julia as possible (I may abuse multiple dispatch a little, oops ). I also like when the code reads like the math (thanks to UTF-8 and VS Code LaTeX shortcuts), and the documentation just above explains the notation and the main equations, like here, for example. Especially useful with in-place matrix operations, in which we deliberately loose on the readability to favor computational efficiency.

I aim for a code quality that I consider still not achievable by any LLMs for now, even with Claude Opus 4.8, GPT-5, etc. That’s our special power…us…mere humans.

[franckgaga]

ModelPredictiveControl v2.5.0

I will release ModelPredictiveControl.jl v2.5.0 soon!

There has been a lot of work under the hood. The underlying QP and NLP problems should be more numerically robust now. All the box constraints are no longer treated as linear inequality constraints (as it was the case before, to simplify the internal logic). They are now defined as proper box constraints, which is, for the MPC:

• the slack variable \epsilon \ge 0
• the input increment bounds \mathbf{\Delta u_{min/max}}
• the terminal constraint \mathbf{\hat{x}}_i(k + H_p), for transcription other than single shooting

and, for the MHE:

• the slack variable \varepsilon \ge 0
• the estimated process noise bounds \mathbf{\hat{w}_{min/max}}
• the arrival estimate \mathbf{\hat{x}}_k(k + N_k + p)

as they should be, to leverage the selected inequality-handling method of the NLP solver, or, to exploit QP solvers that natively support them. Also, the dimensions of the optimization problem is now pretty-printed:

LinMPC controller with a sample time Ts = 4.0 s:
├ estimator: SteadyKalmanFilter
├ model: LinModel
├ optimizer: OSQP 
├ transcription: MultipleShooting
└ dimensions:
  │ ├ 10 prediction steps Hp
  │ ├  5 control steps Hc
  │ ├  2 manipulated inputs u (0 integrating states)
  │ ├  6 estimated states x̂
  │ ├  2 measured outputs ym (2 integrating states)
  │ ├  0 unmeasured outputs yu
  │ └  1 measured disturbances d
  └ optimization: # <--- new section
    ├ 71 decision variables Z̃ (1 slack variable, 21 bounds) 
    ├ 40 linear inequality constraints A (0 custom)
    └ 60 linear equality constraints Aeq

Thanks @cvanaret for the help! Here is the change log:

• added: support box constraint in LinMPC and NonLinMPC
• added: support box constraints in MovingHorizonEstimator
• added: pretty-print optimization problem dimensions for MPC and MHE
• added: covestim keyword argument in both MovingHorizonEstimator constructors
• debug: setconstraint! for MHE and softness parameter arguments now works
• debug: honor the P̂_0 argument at MovingHorizonEstimator construction
• changed: Cwt=1e4 in custom constraints case study
• added: dispatch on covestim type when inverting \mathbf{\bar{P}} in the MovingHorizonEstimator
• doc: added C code generation example in “Manual: Linear Design” (via LinearMPC.jl)
• doc: various minor improvements
• test: improve coverage of predictive controller constraint violation
• bench: MovingHorzizonEstimator benchmark with SteadyKalmanFilter for the arrival covariance
• bench: removed all LinMPC and MHE with Ipopt
• bench: MovingHorizonEstimator with covestim=SteadyKalmanFilter(...)
• bench: add_bridges=true for OSQP
• bench: increase samples for LinMPC

[franckgaga]

ModelPredictiveControl.jl v2.6.0

A small update to announce a new feature similar to the previous one, specializing on the constraint types. The continuity constraints of the OrthogonalCollocation transcription is now treated as linear equality constraints. I’m not sure if I’m the first one to notice that they are structurally linear, presumably not, but I saw no mention of this in all the textbooks, papers and online documentation that I read.

The main gain from this is reducing the size of the Jacobian for the automatic differentiation (AD). The performances gains at NLP solving itself is presumably negligible, except if the solver treats linear constraints separately in a special path. Note there are no significant improvements in the benchmarks for the simple inverted pendulum case study.

On this subject, the defects of the stochastic states (the integrators associated to the unmeasured disturbance model) are also treated as linear equality constraints since v2.1.2, even if the plant model is nonlinear, also to reduce the size of the Jacobian.

Here’s the changelog since the last post:

• added: treat continuity constraints as linear in OrthogonnalCollocation to exclude them from AD
• changed: new API for FastGaussQuadrature and bump its compat entry
• changed: hessian=true default to dense ForwardDiff.jl for MovingHorizonEstimation
• test: improve RungeKutta and LinearMPCext coverage

[franckgaga]

ModelPredictiveControl v2.8.0

The last two updates introduce practical tools for real engineering applications: gracefully handling NaN values in the measurement feedback \mathbf{y^m}. In process control, it’s common that the measured output vector \mathbf{y^m} is a mix of composition/concentration and temperature quantities. The reliability of the former is frequently lower than the latter. For example, online spectrometer will fail if there is dirt accumulation on the sapphire window and return no measurement at all.

When it occurs, all the Kalman filters and the Luenberger observer will skip the correction step altogether and print a warning. This is the equivalent of temporarily opening the control loop. The MovingHorizonEstimator is more sophisticated since it supports partial correction out-of-the-box. The NaN values will be simply ignored when evaluating the objective function, that is, treat them as missing measurements. The InternalModel also supports partial correction by assigning a zero in the estimated stochastic output vector \mathbf{\hat{y}_s^m}.

Here’s the changelog since my last post:

• added: handle NaN values in measured output ym for all state estimators
• added: handle NaN values in measured outputs ym for MovingHorizonEstimator
• debug: truncate :V̂ and :X̂ fields in getinfo(estim::MovingHorizonEstimator)
• debug: do not fill unused decision variables in MovingHorizonEstimator after solve
• doc: update con_nonlinprogeq! method docstring for OrthogonalCollocation
• doc: details on NaN handling in MovingHorizonEstimator and InternalModel

[franckgaga]

ModelPredictiveControl.jl v2.9.0

Huge update today! I am happy to announce that the MovingHorizonEstimator now fully supports MultipleShooting transcriptions!

This was an absolute beast of a feature to implement , see for yourself (colorblind mode FTW):

It was essentially 4 massive jobs rolled into one humongous update:

1. Added support for linear equality constraints in the MHE: This was honestly the hardest part. Getting the math (i.e. the prediction matrices + the state defect matrices) and dimensions right was incredibly tricky mainly because the number of rows in \mathbf{A_{eq} \tilde{Z}} = \mathbf{b_{eq}} dynamically grows at the beginning, when N_k < H_e.

2. Added support for nonlinear equality constraints in the MHE: The implementation relies on the VectorNonlinearOracle of JuMP.jl and DifferentiationInterface.jl for the Jacobian and Hessian.

3. Multiple shooting transcription for linear plant models (QP): Should perform better with unstable plant model compared to the default single shooting, plus the state estimate constraints become box constraints instead of linear.

4. Multiple shooting transcription for nonlinear plant models (NLP): Ironically, setting this up for the nonlinear models actually ended up being easier than the linear ones, except for the multi-threading features (f_threads and h_threads keyword argments) and the potential race conditions.

This is a huge step towards supporting DAEs, since the heavy lifting for the collocation constraints is now established, and a MHE with an OrthogonalCollocation transcription is how I will support DAE systems, at least as a first step. The next step in the short term will be however implementing something similar as v2.1.2 but for the MHE: treating the states of the unmeasured disturbance model as linear equality constraints in the optimization problem.

I will register the update soon.

Happy estimation!

Changelog since the last post:

• added: MultipleShooting support for the MovingHorizonEstimator
• debug: AutoSparse hessian with MovingHorizonEstimator now works
• debug: print warning about NaN values in ym when direct=false
• doc: clearer J prediction matrices (related to measured disturbances d) for MHE
• debug: always update linear eq. constraints for OrthogonalCollocation

[franckgaga]

ModelPredictiveControl v2.11.0

Small update to announce the support of the TrapezoidalCollocation transcription in the MovingHorizonEstimator. It was wayyyy easier than expected . The PR of v2.9.0 really did lay the bricks for any transcription methods in the MHE.

Next step is adding OrthogonalCollocation (should be quite easy) and, just after, the support of DAE systems.

The change log since my last post:

• added: TrapezoidalCollocation for MovingHorizonEstimator
• added: stochastic defects treated as linear equality constraints in MovingHorizonEstimator
• debug: update state constraints in MovingHorizonEstimator based on NonLinModel and MultipleShooting
• debug: update defect matrices in setmodel! for MovingHorizonEstimator and LinModel
• debug: reset prepared flag in initstate!(::StateEstimator) method
• doc: MovingHorizonEstimator data widows in Internals → State Estimators
• doc: prediction and defect matrices of MultipleShooting for MovingHorizonEstimator
• bench: adding new MovingHorizonEstimator cases with MultipleShooting

[franckgaga]

ModelPredictiveControl.jl v2.13.0

Ok, enough about AI debates! Let’s talk about human creations, documentation and API consumed by proper human developers.

Knowing that the v2.12.0 update added the support of OrthogonalCollocation for MovingHorizonEstimator, all the main blocks are now in place for supporting differential and algebraic equations (DAEs). Following up, the v2.12.0 update introduces a new NonLinModelDAE model type

That being said, the package only supports open-loop simulations of NonLinModelDAE objects for now through evaloutput, updatestate! and sim! methods, hence:

DISCLAIMER

NonLinMPC and MovingHorizonEstimator currently do not support NonLinModelDAE objects.

I still wanted to make an announcement here since this is a huge step toward supporting acausal modeling tool like ModelingToolkit.jl more broadly, since they typically produces DAEs instead of ODEs. I would mention here that ODEs are insanely simpler in terms of algorithmic complexity, so a NonLinModel object for ODEs should be used when possible to reduce the dimensions and complexity (thanks to multiple dispatch!).

I also want some feedback on the API for this new type, see the link above. I want to know if you think it is clear and flexible enough. In short, the constructor expect this semi-explicit DAE:

\begin{aligned} \mathbf{\dot{x}}(t) &= \mathbf{f}\Big( \mathbf{x}(t), \mathbf{a}(t), \mathbf{u}(t), \mathbf{d}(t), \mathbf{p} \Big) \\ \mathbf{0} &= \mathbf{q}\Big( \mathbf{x}(t), \mathbf{a}(t), \mathbf{u}(t), \mathbf{d}(t), \mathbf{p} \Big) \\ \mathbf{y}(t) &= \mathbf{h}\Big( \mathbf{x}(t), \mathbf{a}(t), \mathbf{d}(t), \mathbf{p} \Big) \end{aligned}

in which \mathbf{x} is the differential state, \mathbf{a} is the algebraic variable (avoided “z” since already reserved for the decision variable \mathbf{Z}), \mathbf{u} is the manipulated input, \mathbf{d} is the measured disturbance, \mathbf{y} the output and \mathbf{p} the parameter. The function \mathbf{f} and \mathbf{q} are combined into one single Julia function fq!(ẋ, res, x, a, u, d, p) -> nothing to avoid redundant computations. Like all the other model types, the new type only support strictly proper systems w.r.t. \mathbf{u}, since otherwise there is a causality paradox in the flow of information for discrete control. That is why \mathbf{u} is not an argument of \mathbf{h}. I currently use JuMP.jl and Ipopt.jl as a root solver for the open-loop simulations since I did not want to add new unnecessary dependencies. That’s not the scope of the package, it’s closed-loop computations.

The changelog since my last post:

• added: NonLinModelDAE for differential and algebraic equations + sim! support
• added: OrthogonalCollocation for MovingHorizonEstimator
• added: slightly reduce precompilation time
• debug: correctly fill nonlinear equality constraints for non-SingleShooting MovingHorizonEstimator
• doc: using the new DocumenterCodeBlocks.jl plugin
• test: migration to the new test item framework CI workflow
• test: added CI testing on MacOS AArch64
• test: added CI testing on Linux x86_64 (32-bit)