Docs questions

[baggepinnen]

The following is a short list of questions I have tried to answer by reading the docs, but haven't found a clear answer to

• Can I set a unique terminal cost that differs from the running cost? If I don't choose a terminal cost, what cost will be used by default? In particular, I'd like o be able to set the solution to the Ricatti equation of the corresponding LQR problem as the terminal cost matrix.
• Can I set a terminal state constraint that differs from the state constraints during the rest of the trajectory?
• Can I use non-diagonal cost matrices? In particular, my weight matrices may come from an initial LQR design, which may or may not be optimization based or use loop-transfer recovery, resulting in non-diagonal matrices.

[franckgaga]

That's a great summary of the current blind spot in the package. Thanks for that careful revision :)

• No, it's not supported right now (i.e. the terminal cost are always 0). I never had to use terminal cost in my experience. When I had problem with stability, I just increased Hp and it was sufficient. In some rare problematic occasion I used terminal constraints. That is why I skip this part. That's only my own opinion and, for sure, some people would like to use them. I agree that it should be there. I will add this feature as planned and work on the implementation ASAP. Out of curiosity, did you personally used terminal cost in the MPCs that you designed ?
• Yes it's supported, see the last inequality in setconstraint!(::PredictiveController) method (x̂min et x̂max arguments). The trajectory state constraints are in fact the output constraints ymin and ymax.
• No, it's not supported right now. The weight matrices are currently stored as M_Hp::Diagonal{T, Vector{T}}. I could change their type. Should I restrict them to hermitians ? Or just Matrix{T} ? If I allow any kind of matrices, is there a risk that the final QP problem is non-convex (i.e. Hessian with negative eigenvalues) ?

[baggepinnen]

Terminal cost is very useful since it allows the MPC controller to exactly match an infinite-horizon LQR controller as long as no constraints are active. This means that you can use the MPC as a drop-in replacement for LQR, but a replacement that also handles constraints, should they become active. It also means that you can use existing tuning methodologies for LQR controllers and apply them to the MPC controller.

Terminal constraints are also interesting in several formulations of robust MPC under uncertainty, where you want to make sure the system lands in some form of safe set etc.

Should I restrict them to hermitians ? Or just Matrix{T} ? If I allow any kind of matrices, is there a risk that the final QP problem is non-convex (i.e. Hessian with negative eigenvalues) ?

Hermitian is fine. The QP can be nonconvex already when you store the diagonal only, if the diagonal entries are negative. What's required is that the matrices are positive semi definite for the state penalties, and positive definite for the input penalties.

[franckgaga]

Damn that seem great at least on paper. But yes combining LQR and MPC constraints seems quite powerful Edit: alright I read a little bit on MPCs that act as LQR but with constraints. Seems really powerful, I'll clearly add these features in the next release.

[franckgaga]

I was thinking about how to allow terminal cost for "LQR tunings". I'm hesitating between two options:

1. Only modify the code to allow completely custom M_Hp and L_Hp weights (restricted to hermitians). The last block in these two block matrices correspond to the terminal step. This way, the terminal state x would not explicitly appears in the cost function, only implicitly through the output y and the manipulated inputs u. From my understanding, this is how Matlab does it, and it needs some clever but not obvious tricks to replicate a LQR, see https://www.mathworks.com/help/mpc/ug/using-terminal-penalty-to-provide-lqr-performance.html
2. Also modify the code to allow completely custom M_Hp and L_Hp weights, but also add a new x'*T*x term in the cost function, where x is the terminal state. This would simplify LQR tuning, but the user would need to remember to put the last blocks in M_Hp and L_Hp matrices to zero. I could possibly issue a warning if not.

Do you have a preference, or other ideas ?

[baggepinnen]

Typically, the terminal cost and state constraints are handled slightly differently to the rest of the trajectory, and we do not consider a terminal penalty on the control input. That is, the penalized state trajectory is of length N+1 while the penalized input trajectory is of length N. This is done simply because at the terminal time, we are "done" and do not perform any action. This is also why the Bellman cost-to-go depends on the state only, i.e. $x^T P x$, and it is this cost-to-go that we use as terminal cost in order to have the linear quadratic MPC exactly match LQR.

In JuliaSimControl, we take this approach, which most closley resembles your second option above, with the difference that there is no need to set any blocks to zero, since the terminal cost relates to the state one time step after the trajectory that relates to M_Hp, L_Hp. This is likely how matlab also implements it, even if it's not obvious from the docs.

You'll the same approach taken in the docs example for OSQP, where the implementation is simple enough to inspect https://osqp.org/docs/release-0.6.3/examples/mpc.html notice how the quadratic part in the QP cost function is

P = blockdiag(kron(speye(N), Q), QN, kron(speye(N), R))

i.e., there are N terms for both state and control, plus the addition of the QN term for the terminal state in the middle

[franckgaga]

Alright thanks for the details! :)

I inspected in details the solution of JuliaSimControl and OSQP. It is in fact the first approach that I mentioned, but my explanation was a bit shaky XD. If you look carefully in the doc (at the beginning) the output tracking term (with M_Hp weight) goes one time step further than the input (with L_Hp weight) and input increment (with N_Hc) terms (strickly speaking, the input increment term "finishes" even sooner because of the control horizon Hc). If we allow a completely custom M_Hp, the last block can be modified to include the terminal weight, leading to the objective function of JuliaSimControl. I'll implement that.

[franckgaga]

I just registered a new version that allows non-diagonal weights, terminal cost and LQR tuning. I added a comment on terminal costs in LinMPC doc (there is an issue with doc deployement, see #26): I also added a unit test on a second order system that compares the closed-loop trajectories of a LQR controller and a LinMPC controller with terminal weights computed by the Ricatti equation. This update should cover all your questions raised here.

[franckgaga]

New version should solve the issues.