Flexible linear constraints for MPC and codegen

[franckgaga]

Following the subject [ANN] ModelPredictiveControl.jl with @langestefan and @darnstrom:

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.

Agree. I would like to support more flexible linear constraints in ModelPredictiveControl.jl.

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

And that’s also why I think running MPC on traditional embedded systems (e.g. not Raspberry Pi or Nvidia Jetson) is not ideal. Even by working very hard to make C codegen fully feature-complete, there is still one unsolvable issue: it’s another form of the two language problem. A major problem that comes with it is reproducibility (the numerical results will be presumably different, and possibly drastically different for ill-posed corner cases). But that’s another subject.

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.

Since right now I rely on LinearMPC.jl for codegen, I need to include @darnstrom in the loop. My understanding is online updates of e.g. the lb and ub arguments of add_constraint! within the C code is not possible right now, am I right? Would it be something that is possible to add in LinearMPC.jl ?

About this @langestefan, do you think that the add_constraint! API provided by LinearMPC.jl is flexible enough for your specific constraint structure? Let’s assume that you are also able to update the value of lb and ub online in the generated C code.

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.

It would not help here since I’m using LinearMPC.jl for C codegen, and JuMP is not used at all in this package.

[darnstrom]

Actually, I created an interface to JuMP for ParametricDAQP.jl a month ago (see example in README). What it does is that it creates a multi-parametric quadratic program based on a JuMP model, which it then used to create an explicit solution.

My intention was also to add support for codegen also (i.e., given an mpQP → generate code for DAQP, which is basically what goes on underneath the hood of LinearMPC.jl.) I have a local branch with some progress here, but it was down-prioritized due to other projects. It would, however, not be too much work to finalize it. If there is interest from others (e.g., @franckgaga and @langestefan) I could bump the priority of that action point.

[langestefan]

@franckgaga

And that’s also why I think running MPC on traditional embedded systems (e.g. not Raspberry Pi or Nvidia Jetson) is not ideal.

Not ideal, but definitely feasible! In my case I am trying to add MPC to an existing product that is currently using basic PI control. Which does the job, but we can do better. I think there are hundreds of companies out there in similar situations. And personally I think the whole field of doing complicated stuff on limited hardware is very exciting and rewarding.

@franckgaga

Even by working very hard to make C codegen fully feature-complete,

So that already exists: cvxpygen does an amazing job. If you look at the source code it’s really not that complicated, it’s essentially just a wrapper around cvxpy that packages your solver of choice with a nice interface. It is ofcourse a generic tool and not purpose-made for MPC applications, which is where we can do a better job

@franckgaga

there is still one unsolvable issue: it’s another form of the two language problem.

Yes it is. I think in most situations you will have to resort to C codegen because the target MCU is also doing 10 other things.

@franckgaga

A major problem that comes with it is reproducibility (the numerical results will be presumably different, and possibly drastically different for ill-posed corner cases). But that’s another subject.

Not sure I am following you here, reproducibility of what?

@franckgaga

About this @langestefan, do you think that the add_constraint! API provided by LinearMPC.jl is flexible enough for your specific constraint structure? Let’s assume that you are also able to update the value of lb and ub online in the generated C code.

Yes, that looks quite flexible. Not entirely sure what this part means:

(additional terms Ar rₖ, Aw wₖ, Ad dₖ, Aup u⁻ₖ are possible)

But I will have a look at the source code.

[franckgaga]

Interesting, but would it solve the main issue discussed here? That is, would it be possible to change the upper bound and lower bound of the linear inequality constraints online in the generated C code (based on the DAQP solver)?

[langestefan]

@darnstrom

Actually, I created an interface to JuMP for ParametricDAQP.jl a month ago (see example in README). What it does is that it creates a multi-parametric quadratic program based on a JuMP model, which it then used to create an explicit solution.

If I understand correctly we would be doing explicit MPC? My experience with explicit methods is that they only work for (very) small problems, the codegen quickly explodes in size. Happy to be proven wrong.

[franckgaga]

@langestefan

Yes, that looks quite flexible. Not entirely sure what this part means:

(additional terms Ar rₖ, Aw wₖ, Ad dₖ, Aup u⁻ₖ are possible)

The supported structure of add_constraint! is (using non-boldface font for scalars and boldface font for vectors and matrices):

\underline{b_{k}} - \epsilon_k \le \mathbf{A_x} \mathbf{x}_{k} + \mathbf{A_u} \mathbf{u}_{k} + \mathbf{A_{up}} \mathbf{u}_{k-1} + \mathbf{A_r} \mathbf{r}_{k} + \mathbf{A_d} \mathbf{d}_{k} + \mathbf{A_w} \mathbf{w}_{k} \le \overline{b_k} + \epsilon_k

in which :

• \underline{b_{k}} is the lower bound (lb argument)
• \overline{b_{k}} is the upper bound (ub argument)
• \epsilon is the slack for softening (can be disabled)
• \mathbf{x}_{k} is the state at time step k
• \mathbf{u}_{k} is the input at time step k
• \mathbf{u}_{k-1} is the input at time step k-1
• \mathbf{r}_{k} is the setpoint at time step k
• \mathbf{d}_{k} is the known disturbance at time step k (if applicable, for feedforward comprensation)
• \mathbf{w}_{k} is the unknown disturbances at time step k (if applicable, for robust MPC)

You can add an “infinite amount” of these constraints (with distinct lb, ub, Ax, Au, … values), and for each of them you select the k value that applies (the ks argument. The present time step is 1 and the terminal time step is Np+1). Is it clearer ?

[darnstrom]

Any parameter that enters linearly in the objective / constraint will be handled, which is the case for upper/lower bounds.

[darnstrom]

The current code for ParametricDAQP.jl is for explicit solutions to mpQPs, but what I was trying to say is that I also aim to add support for an implicit solver (DAQP). With that support, you will be able to generate solver code (implicit or explicit) for any JuMP model that is of the form of a multi-parameteric QP.

[langestefan]

Wow! I think that is super exciting. Would that potentially support the B&B method as well?

[darnstrom]

Yes!

[franckgaga]

So FYI I’m working on adding the support of custom linear inequality constraint in ModelPredictiveControl.jl. It will support C codegen, except for online modification of \mathbf{G_{min}} and \mathbf{G_{max}}, at least for now. The API will looks like this (see the last inequality):