[ANN] ModelPredictiveControl.jl

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 :slight_smile:

1 Like

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 ?

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.

1 Like

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.

1 Like

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.

2 Likes

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 :slight_smile:

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

7 Likes

:bottle_with_popping_cork::tada: ModelPredictiveControl.jl v2.0.0 :bottle_with_popping_cork::tada:

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

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
3 Likes