Multiple shooting transcription

[1-Bart-1]

It would be amazing if we could implement the algorithms for Real-Time Iterations for Nonlinear MPC and MHE from this paper: https://scholar.google.no/citations?view_op=view_citation&hl=en&user=38fYqeYAAAAJ&citation_for_view=38fYqeYAAAAJ:RYcK_YlVTxYC
From the paper:
"Real-time methods for MPC and MHE such as the RTI exploit
the similarity of the NLPs underlying the MPC and MHE from one
sampling time to the next. Indeed, for a reasonably high sampling
frequency, the parameters (estimated states and parameters) en-
tering the NLPs do not change significantly from one time sample
to the next, and the resulting solutions to the NLPs are very similar.
The solution of the NLP at a sampling time T_i is therefore used as
an initial guess for the solution of the NLP at the next time instant
T_(i+1) with the aim to maintain a fast rate of convergence at all time
instants. In that context, the RTI scheme relies on taking a single
Newton-type iteration at every time instant, always based on the
latest system information. Consequently, the method produces
locally sub-optimal solutions. Careful initialization strategies with
shifting and initial value embedding ensure that the good con-
traction properties of the Newton-type iterations are preserved"

This solution solves the problem where Nonlinear MPC is slower than the system's sample time, which is a common challenge in many control applications (for example: https://discourse.julialang.org/t/how-fast-does-a-model-has-to-be-for-nmpc/120694).

[baggepinnen]

Implementing RTI with a generic choice of solver is very difficult. In particular, it is not feasible when using Ipopt. If you read the papers you're citing closely, you see that they use a very particular and purpose-built solver.

[franckgaga]

RTI would be in long term goals. And yeah it's not clear how to implement this feature using the generic solver-independent interface of JuMP.

In the midterm and shorterm, there is other enhancement that should improve the NMPC speed:

1. When JuMP will support vectors in user-defined operators. It is planned but it's a humongous job for them. It should improve the performance a little bit since right now there is conversion between tuples and vectors for the decision variable.
2. Multiple shooting methods. I only support single shooting right now.
3. Better ODE solvers (will work on that in the short-term)

[1-Bart-1]

I can try to work on this myself, but it would be much appreciated if I could get some help pointing me in the right direction.

I can implement the RTI specific solver. But then it maybe doesn't make sense to register this solver to JuMP, because the RTI isn't solver-independent anyways?

[baggepinnen]

I would not attempt to implement such a solver, it is an absolutely enormous task if you'd like it to be anything but a toy. If you need this, I'd use acados instead. There's a proof of concept interface here
https://github.com/baggepinnen/AcadosInterface.jl
but you could also just implement the model in Python-casadi directly.

[1-Bart-1]

Thanks, I will try using acados!

[baggepinnen]

Do you know what aspect of your problem makes it too slow? Is it very large state dimension, is it stiff, does it have complicated constraints etc.?

[1-Bart-1]

It is a stiff ODE (stiff tether segments simulation).

[baggepinnen]

In that case it might be sufficient to use an integration method for stiff solvers, which would avoid having to take super small steps. Oftentimes, trapezoidal integration is sufficient, the simplest method that's able to handle stiff systems.

[franckgaga]

Yeah that's what Matlab uses by default for nlmpc, an implicit trapezoidal solver. @baggepinnen do you have a suggestion for a good lightweight nonlinear root solving package? I used NLsolve.jl in the past and it was working well and also lightweight.

[baggepinnen]

You need a way to differentiate through the root find for it to be efficient. SimpleNonlinearSolve.jl defines the required chain rules for this to work using the implicit function theorem

[franckgaga]

Alright thanks I'll try to implement a new solver during the Christmas holiday!

[baggepinnen]

The most efficient way to implement trapezoidal discretization tends to be to implement it as constraints to IPOPT. This way, you don't need any external root finder, Ipopt will handle everything through the constraints
$$F(x, u) = 0$$
This can also cut down the number of function evaluations by a factor of 2 (provided that you accept FOH instead of ZOH interpolations of the inputs), since you only need to evaluate the dynamics once per time step, whereas if you implement it as a time-stepping integrator with ZOH interpolation, you cannot exploit the fact that you have already evaluated $f(x_k, u_k)$ when you compute the step between $f(x_k, u_k)$ and $f(x_{k+1}, u_{k+1})$ when you computed the step $f(x_{k-1}, u_{k-1})$ to $f(x_k, u_k)$

[franckgaga]

Oh I did not though of that! That's also great since no other additional dependency. But that would also means that open-loop simulations the a NonLinModel instantiated with this ODE solver would not be possible. Only closed-loop in the MPC.

While I'm at it, I would need to use another ODE solver for the StateEstimator...

edit: maybe adding a lightweight nonlinear root solver package as a new dependency, but using it only for open loop simulations and for the state estimtator 🤔