New jump nlp syntax, merge into main?

[franckgaga]

JuMP.jl is introducing a new syntax for NLP, see https://jump.dev/JuMP.jl/stable/manual/nonlinear/ warning at the top.

It is cleaner and marginally faster from my quick tests (the speed difference is negligible in fact). The drawback is that the optimizer support is partial right now. Ipopt.jl works well (the default for NonLinMPC) but NLopt.jl does not work. IMO, this meta-package include many interesting solvers for NMPC e.g. a SQP method called SLSQP.

I finished the migration to the new syntax (both for NonLinMPC and MovingHorizonEstimator). Come to think on it, maybe we should wait for at least NLopt.jl support before merging into main. What's your opinion ?

[codecov-commenter]

Codecov Report

All modified and coverable lines are covered by tests ✅

Comparison is base (28ccc09) 98.90% compared to head (60db884) 98.98%.

Additional details and impacted files

@@            Coverage Diff             @@
##             main      #31      +/-   ##
==========================================
+ Coverage   98.90%   98.98%   +0.08%     
==========================================
  Files          23       23              
  Lines        2465     2466       +1     
==========================================
+ Hits         2438     2441       +3     
+ Misses         27       25       -2     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[baggepinnen]

I don't have any strong opinions on this, just a few comments

but NLopt.jl does not work. IMO, this meta-package include many interesting solvers for NMPC e.g. a SQP method called SLSQP.

I've evaluated NLopt for MPC before and found it to be a poor fit. To solve challenging nonlinear MPC problems {large, long horizon, stiff dynamics, algebraic equations}, one usually require a multiple-shooting or direct-collocation transcription of the dynamics (as opposed to the single-shooting transcription used in this package so far). Such transcriptions result in large but sparse optimization problems, where the Lagrangian Hessian and the constraint Jacobian are both very sparse. NLopt does not support sparsity, and thus scales very poorly to large problems, whereas Ipopt is a sparse solver and handles it fine. If the intention with this package is to eventually support any such transcription methods, NLopt thus becomes irrelevant for all but very small problems.

Unfortunately, the new JuMP nonlinear interface doesn't yet support vector-valued nonlinear functions, which means that it's very awkward to integrate continuous-time dynamics with anything other than very primitive integrators like forward Euler. Forward Euler is a horrendously bad integration method that should in principle never be used, even a simple RK4 is incredibly much better. However, in JuMP, one has to integrate each equation in the dynamics scalar-wise using RK4, rather than having the integrator coded up in a function since JuMP cannot represent this function from $\mathbb{R}^{n_x}, \mathbb{R}^{n_u} \to \mathbb{R}^{n_x}$. This makes JuMP scale poorly to large problems with multiple shooting. This is a large part of why I didn't use JuMP for the MPC package in JuliaSimControl, and also why I was a it disappointed this summer when the new nonlinear interface was presented, the only thing I was missing was still missing :P

Direct collocation transcription could probably be implemented in JuMP since you do not suffer the increase in symbolic complexity from recursive application of the dynamics in the integrator, I haven't tried that but it should be relatively straightforward.

[baggepinnen]

To show precisely how bad forward Euler is, here's a comparison using the model from the first example of nonlinear design in the documentation of this package

using ModelPredictiveControl, SeeToDee

function pendulum(x, u, par, t)
    g, L, K, m = par        # [m/s²], [m], [kg/s], [kg]
    θ, ω = x[1], x[2]       # [rad], [rad/s]
    τ  = u[1]               # [N m]
    dθ = ω
    dω = -g/L*sin(θ) - K/m*ω + τ/m/L^2
    return [dθ, dω]
end
# declared constants, to avoid type-instability in the f function, for speed:
const par, Ts = (9.8, 0.4, 1.2, 0.3), 0.1
f_euler(x, u, _ ) = x + Ts*pendulum(x, u, par, 0) # Euler method
rk4 = SeeToDee.Rk4(pendulum, Ts)
f_rk4(x, u, _ ) = rk4(x, u, par, 0) # RK4 method



h(x, _ )    = [180/π*x[1]]  # [°]
nu, nx, ny = 1, 2, 1
model_euler = NonLinModel(f_euler, h, Ts, nu, nx, ny)
model_rk4 = NonLinModel(f_rk4, h, Ts, nu, nx, ny)

using Plots
u = [0.5]
N = 35
res_euler = sim!(model_euler, N, u)
res_rk4 = sim!(model_rk4, N, u)
plot(res_euler, plotu=false, legend=true, label = "Euler")
plot!(res_rk4, plotu=false, legend=true, label = "RK4", c=2)

the example simulates the pendulum using forward Euler and compares that to a simulation using RK4, the difference is rather big

If the sample time Ts is increased to Ts = 0.2, forward Euler becomes unstable, whereas RK4 remains accurate.

[franckgaga]

the example simulates the pendulum using forward Euler and compares that to a simulation using RK4, the difference is rather big

If the sample time Ts is increased to Ts = 0.2, forward Euler becomes unstable, whereas RK4 remains accurate.

I know, forward Euler is terrible 😢 . I used this method in the manual to keep the explanation as simple as possible. I wanted to put the focus on estimator and controller design, not model solving. I typically also used Runge-Kutta in the NMPC that I implemented on real systems. I could maybe add super-samples in the example (i.e. multi-step Euler), since we are near instability.

Direct collocation transcription could probably be implemented in JuMP since you do not suffer the increase in symbolic complexity from recursive application of the dynamics in the integrator, I haven't tried that but it should be relatively straightforward.

I know about direct collocation but I never dug deep in it. I will read on this tomorrow. I would be interested to add it in the package, in the short term if it's not too hard. I don't plan to add multiple shooting methods in the short term. I like the simplicity and elegance of the single shooting method. Easy to understand also means easy to troubleshoot. I understand that it scales poorly, but not everyone want to solve large problems. And its like linear algebra packages: sparse formulations are faster for large sparse problems, but they are also generally slower for small/medium problems, compared to dense formulation. The direct collocation method could maybe fill this gap in the package.

However, in JuMP, one has to integrate each equation in the dynamics scalar-wise using RK4, rather than having the integrator coded up in a function since JuMP cannot represent this function from $\mathbb{R}^{n_x}, \mathbb{R}^{n_u} \rightarrow \mathbb{R}^{n_x}$.

I'm not sure I fully understand this part. JuMP does not directly support functions with vector outputs, that's true (it is planned though), but there is a hacky workaroud. I do not work with scalars in the package, otherwise the user would not be able to change the prediction and control horizons.

[baggepinnen]

I generally prefer direct collocation over multiple shooting due to the simplicity with which one can

1. Handle stiff systems
2. evaluate any function of the state (such as constraints) at a greater temporal resolution than the time step of the MPC controller

Multiple shooting and direct collocation are otherwise very similar, direct collocation includes all the intermediate stages in an integration step as variables, whereas multiple shooting only includes the initial condition for the shooting segment as variables. The complexity in their implementation is almost identical, with the exception that you need to handle the indexing of the intermediate stages for collocation. These intermediate states is also what makes it very cheap to evaluate functions of the state, since the state along the trajectory is treated as explicit variables rather than as computed from the input.

I do not work with scalars in the package

from the perspective of the user, calling f(x::Array) where x is a JuMP variable is an array operation, but JuMP will trace through f using the scalar entries of the symbolic array x, instead of treating f(x) as a function that takes a vector. This is an important limitation since f may be very complex (e.g., containing multi-step integration or linear systems solves A\b etc.), which may cause the complexity of the symbolic representation JuMP uses to blow up, become very slow and ultimately stack overflow.

[franckgaga]

New Idea: I could just support nonlinear continuous state-space models with a new syntax, and add SeeToDee.jl as a new dependency.

[baggepinnen]

That would also work. SeeToDee only contains 2 integrators at the moment, RK4 works with JuMP wile SimpleColloc does not. The implementation of RK4 is really simple, so it would almost be better to just include this in the package rather than taking on a dependency on SeeToDee, which has a non-trivial dependency tree.

It's not quite a replacement for direct collocation though, since RK4 (and other explicit integrators) does not handle stiff equations

[baggepinnen]

By the way, have you ever found SL_SQP from NLopt outperform Ipopt? Even for very small problems, every time I've tried I've found SL_SQP to be suffering from convergence problems and taking 100x longer than Ipopt to find a good solution. I think that the problem might be that SLSQP does not make use second-order information, instead the keep a BFGS-like approximation to the Hessian. Ipopt through JuMP will use AD Hessians, which probably explains the difference.

https://nlopt.readthedocs.io/en/latest/NLopt_Algorithms/#slsqp

[franckgaga]

I only tested it on the pendulum example and, yes, it was about two times slower than Ipopt. It is safe to assume that the "exact" Hessian probably does help here. Btw, Optim.jl is now available in JuMP. There are some interesting local solvers here that use AD but I did not tested them.

[franckgaga]

@baggepinnen so I made some quick test on some NLP solvers:

• NLopt.jl: most solvers are just not able to solve the pendulum problem incl. SLSQP (seems to stuck easily on local minima, etc.). The package still does not support the new NLP syntax, but, yes, it seems that they are of minor interest for NMPC, let's just ignore it.
• Ipopt.jl: works perfectly with the new syntax
• MadNLP.jl: works perfectly with the new syntax
• KNITRO.jl: works perfectly with the new syntax

Seems like a good time to migrate to the new NLP syntax and merge into main.

[baggepinnen]

Your findings seem similar to mine 👍

Seems like a good time to migrate to the new NLP syntax and merge into main.

🎉 😃