🚀 added: C code generation via LinearMPC.jl

[franckgaga]

Closes #124, for LinMPC controllers.

I implemented the conversion code as a package extension.

It comes with some limitations, as described in the docstring:

• the solver is limited to DAQP.
• the transcription method must be SingleShooting.
• the state estimator must be a SteadyKalmanFilter with direct=true.
• only block-diagonal weights are allowed.
• the constraint relaxation mechanism is different, so a 1-on-1 conversion of the soft constraints is impossible (use Cwt=Inf to disable relaxation).

I also mentioned in the docstring that the reverse is also true: LinearMPC.jl offers some exclusive features that ModelPredictiveControl.jl do not support e.g. binary manipulated inputs and constrained explicit MPC.

[darnstrom]

@darnstrom Are you planning to register a new release of LinearMPC.jl in the short term ? I cannot set a proper [compat] bound on dev version, it needs a proper SemVer number. Thanks!

I'm writing some tests today and I will be ready to release the interface afterward.

I aim to release v0.8.0 of LinearMPC.jl during this week!

[darnstrom]

Btw, you should be able to set rate constraints quite easily with:

add_constraint!(mpc; Au=I(nu), Aup=-I(nu), lb=dumin, ub=dumax, ks=1:mpc.Nc)

where dumin and dmax are your upper and lower bounds, respectively, on $\Delta u$

[franckgaga]

Oh thanks, I did not noticed the Aup argument. The "p" stand for previous I assume ? So the equation is:

$$\Delta u_{min} \le A_u u_k + A_{up} u_{k-1} \le \Delta u_{max}$$

is that it ?

[darnstrom]

$$\Delta u_{min} \le A_u u_k + A_{up} u_{k-1} \le \Delta u_{max}$$

is that it ?

Correct! It is not that well documented, but you can also include the reference and disturbance in the constraints with Ar and Ad, respectively.

[franckgaga]

Btw, you should be able to set rate constraints quite easily with:

add_constraint!(mpc; Au=I(nu), Aup=-I(nu), lb=dumin, ub=dumax, ks=1:mpc.Nc)

where dumin and dmax are your upper and lower bounds, respectively, on Δ u

Ok I added to conversion of input increment constraint. Thanks for the tips.

But I noticed something else. Quoting the new text above, in the first comment:

• the constraint relaxation mechanism is different, so a 1-on-1 conversion of the soft constraints is impossible (use Cwt=Inf to disable relaxation).

I understood the you rely on the internal constraint softening mechanism of DAQP, and the soft_weight is simply used to configure the rho_soft option of DAQP.

In MPC.jl I explicitly introduce a new slack variable ϵ in the optimization problem. In other words, if I use the DAQP optimizer through JuMP.jl, there will be two slack variables: one introduced by MPC.jl, the other introduced internally by DAQP.jl. I did this since I wanted to support constraint softening even if the QP solver does not support it. I also needed an approach that is agnostic of the solver and its settings.

Knowing this, I deduced that it will be simply impossible to perfectly convert LinMPC controller to LinearMPC.MPC controller if softening is activated somewhere. So I simply added this @warn in such cases:

@warn "The LinearMPC conversion applies an approximate conversion " *
      "of the soft constraints.\n You may need to adjust the soft_weight "*
      "field of the LinearMPC.MPC object to replicate behaviors."

Do you feel like it's a good solution ?

[franckgaga]

BTW this is from my experience: I'm able to reproduce the results of LinearMPC.jl up to the 10th digits if softening is disable. Otherwise, the results can be vastly different if a constraint is activated, and manually playing with soft_weight after conversion can mitigate this problem.

[franckgaga]

About #307 (comment)
Do you have any better idea @baggepinnen ? I'm not entirely happy with this solution, but I don't see any better workaround...

[baggepinnen]

How about turning off DAQP softness and only use the softness implemented here?

[franckgaga]

How about turning off DAQP softness and only use the softness implemented here?

Thanks for the idea. AFAIK, this is impossible for 2 reasons:

1. LinearMPC.jl would need to support introducing a new variable in the QP problem (the slack). This is not supported right now. The decision variables are the manipulated inputs over the horizon, nothing else. Maybe it's possible but I think it would require a lot of work from @darnstrom, since he would need to implement this feature both in Julia and the generated C code.
2. I don't think we can completely disable the internal relaxation mechanism in DAQP solver.

[darnstrom]

Are you using a single slack variable for all soft constraints @franckgaga?
DAQP's handling of soft constraints is equivalent to a slack variable for each soft constraints. Although, such variables are never introduced explicitly but handled implicitly (reducing computations.) For some details see Section 3.1 in this preprint: https://arxiv.org/pdf/2512.18458

[franckgaga]

Yes, and that's the same approach than MATLAB (equal concern for relaxation). There is a single slack $\epsilon$ for all the constraints. And then, the softness parameters $\mathbf{c_{(\bullet)}}$ can be set to 0 to specify a hard bound, 1 to specify a soft bound, 0.5 to specify a "somewhat soft" bound, etc.

It's entirely possible that the way you handle soft bounds is more efficient, thus more adapted to low-power embedded platform. It's another sign of the tradeoff mentionned in #124. Full reproducibility is hard with code generation, especially if the targeted platform is vastly different in term of computational power.

So I think that the best solution is me trying to find a good heuristic for a conversion factor between Cwt in MPC.jl and soft_weight in LinearMPC.jl (by trial-and-error), and showing the warning to notice the user that it may still be not a good heuristic. Thanks for the explanations on DAQP internals. It will help in finding this approximate conversion factor.

[darnstrom]

Ok! Then the difference make sense to me.

A problem with one slack variable is that some constraints can be violated even though it is not necessary: if one constraint needs to be violated, resulting in epsilon > 0, this will allow all other soft constraint to be violated without extra cost. This is corrected by having a slack variable for each soft constraint. The obvious drawback with one slack variable for each soft constraint is that a lot more extra decision variables are required (although, as mentioned above, this is not the case for DAQP due to some tricks internally.)

If one is forced to add slack variables explicitly (which is the case for MPC.jl that wants to support any generic QP solver) it make sense to just use one slack variable for computational reasons!

[franckgaga]

A problem with one slack variable is that some constraints can be violated even though it is not necessary: if one constraint needs to be violated, resulting in epsilon > 0, this will allow all other soft constraint to be violated without extra cost. This is corrected by having a slack variable for each soft constraint. The obvious drawback with one slack variable for each soft constraint is that a lot more extra decision variables are required (although, as mentioned above, this is not the case for DAQP due to some tricks internally.)

Super interesting and well-explained 👏 !

[franckgaga]

Okay some final comments:

1. You are entirely right @darnstrom, my solution tends to violate more the soft constraints, even if it would not be necessary for feasibility.

2. By trial-and-error on a SISO system, I found:

   conversion_factor = 1/2/n_soft
   soft_weight = conversion_factor*Cwt

   It approximately replicates the solutions in various scenario of constraint softening. Since DAQP implicitly introduces a new variable for each softened bounds, it make sense that the conversion factor is inversely proportional to the number of softened bounds. This is the best I could find. For the rest there is the warning to notice the user that it may not be a good choice.

3. The issue that we are discussing is less apparent when the softened constraints are the output variable $\mathbf{y}$ only, which is the most common case in practice, IMO.

[franckgaga]

Quick update, by playing a little bit more with my SISO case, I found that:

soft_weight = 10*sqrt(nsoft*Cwt)

is a better conversion, that is, the solutions are more similar in closed-loop (up to the 3rd digits, and it was up the the 1st digit before this change). And it does make sense to include a sqrt here since the $\rho$ weight is squared, according to eq. (10) in "Prioritized Constraints in Optimization-Based Control", which is not the case for Cwt weight in MPC.jl.