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 🤔

[baggepinnen]

Yeah, the state estimator would need a root solver. Open loop simulation could use the same transcription as when doing optimization, but with x as the only desicion variables and u fixed. Ipopt could then solve this as a pure constraint satisfaction problem. You could in principle do this in the state estimator as well, but it would perhaps feel a bit heavyweight for a single step integration

[franckgaga]

@1-Bart-1 you can try to pull the dev branch. I did not had the time to implement and implicit trapezoidal method, but I improved the RungeKutta(4) method a little bit, and I added a ForwardEuler solver. It's not a stiff solver at all (lol), but in my experience it is sometimes more efficient than RK4 for approximate simulations, because of its uber-simplicty (please use the supersample argument haha 👍 )

[franckgaga]

So adding the trapezoidal solver is much more involving than I first thought. From my understanding, solving the ODE in ipopt equality constraint means that it needs to be in the multiple-shooting form, since the decision variable are now both the manipulated inputs and the states over the prediction horizon. Am I right @baggepinnen ? If so, I would add the support of multiple shooting method first.

[baggepinnen]

since the decision variable are now both the manipulated inputs and the states over the prediction horizon

yes

multiple-shooting form

I would call this the "direct collocation form", since this is how you would implement a collocation method (of which trapezoidal integration might be the simplest). Multiple shooting is typically associated with a stepper like RK4, but breaks the trajectory up in multipel segments. They are thus similar, but not identical.

[franckgaga]

Ok but independantly from the ODE solver, the MPC optimization problem can be expressed with or without the state in the decision variables. How would you call these two forms ? The are called sparse and dense formulation in MATLAB doc:

I need to give it a name for the argument in NonLinMPC constructor, a name independent of the specified solver in the plant model. How would you call these two forms ?

[baggepinnen]

sparse/dense does not distinguish between multiple shooting and direct collocation, which both lead to sparse problems

"single shooting", "multiple shooting" and "direct collocation" is usually use to distinguish the three most common transcription methods.

Collocation methods differ from MS in general since they include also the intermediate stage variables as optimization variables, i.e., an RK4 method implemented as direct collocation would include 4*nx optimization variables per time step, while MS would include nx variables and single shooting would include no state variables at all. Trapezoidal integration being the simplest one happen to need no more variables than MS, but that's a special case

[franckgaga]

Another question, I want to fully understand the difference between the 3 common transcription method before starting the implementation, and direct collocation methods are still not clear to me:

I would call this the "direct collocation form", since this is how you would implement a collocation method (of which trapezoidal integration might be the simplest). Multiple shooting is typically associated with a stepper like RK4, but breaks the trajectory up in multipel segments. They are thus similar, but not identical.

By stepper you meant an explicit ODE solver ? So, from my understanding:

• Explicit ODE solver + state in the decision variables : multiple shoothing transcription
• Implicit ODE solver + state in the decision variables : direct collocation transcription
  Do you agree ?

[baggepinnen]

Multiple shooting is essentially calling solve(ode) on multiple segments of the trajectory, with constraints making sure each segment ends where the next starts, and each segment start becomes an optimization variable. Direct collocation adds the internal stages of the ODE solver as optimization variables and there is no call to any solver left anymore, e.g., the following rk4 implementation translated to a collocation method

function rk4(f::F, Ts) where {F}
    # Runge-Kutta 4 method
    function (x, p, t)
        f1 = f(x, p, t)
        f2 = f(x + Ts / 2 * f1, p, t + Ts / 2)
        f3 = f(x + Ts / 2 * f2, p, t + Ts / 2)
        f4 = f(x + Ts * f3, p, t + Ts)
        x += Ts / 6 * (f1 + 2 * f2 + 2 * f3 + f4)
        return x
    end
end

would include not only x as an optimization variable, but also f1,f2,f3,f4 (alternative but equivalent formulations are a bit more common). I'd recommend section 8.5 and 8.5.3 in particular in
"Model Predictive Control: Theory, Computation, and Design, 2nd Edition", Rawlings, Mayne, Diehl
for all the juicy details :)

Chapter 10 in "Nonlinear Programming, Concepts, Algorithms, and Applications to Chemical Processes", Lorenz T. Biegler
is also helpful

[franckgaga]

Ooooh yes I see, no solver is called anymore in direct collocation. Thanks for the ref, I'll take a look at them!

[franckgaga]

Just a quick follow-up : since everything is "hard-coded" as a single shooting transcription in the package right now, the implementation of other transcription method is a big task. I'm working on implementing the multiple shooting method first, since it's a bit easier than direct collocation.

Edit: both LinMPC and NonLinMPC will support the multiple shooting transcription, since the OSQP solver can also leverage its sparse structure. That's also why the task is big, it's almost two distinct implementations.