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
. 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