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:
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.
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.
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.
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!).
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.
This is major release to introduce the new built-in OrthogonalCollocation transcription. There are also two breaking changes:
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.
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 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 !
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
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).
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.
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
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
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
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:
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.
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.
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.
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
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
. It will be available in v1.16.1 before the end of the day.
ModelPredictiveControl.jl v2.0.0 
). 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 
, see for yourself (colorblind mode FTW):
. The PR of v2.9.0 really did lay the bricks for any transcription methods in the MHE.
