Sparsity pattern

[0Yassine0]

@amontoison @frapac @ocots @jbcaillau @PierreMartinon
There might be an issue due to the order in which we set up variables with OC.
The sparsity pattern and x0 that we obtain with OC differ from those we get with JuMP.

In this file, we can see that for the rocket problem, we have:
With JuMP : x0 = [x1[1],x2[1],x3[1],x1[2],x2[2],x3[2],u,v]
With OC : x0 = [x1[1],x1[2],x2[1],x2[2],x3[1],x3[2],u,v].

The problem is that once the problem is defined, I don't think that we can change the order of the variables.
We can't inject JuMP's x0 into the nlp of OC because that would provide the wrong initial guess for the variables.

[frapac]

I agree that the ordering of the variables is of primal importance.
Can you point to me how OC is ordering its variables internally?

[0Yassine0]

I'm not sure about it in the new organization of control-toolbox but probably somewhere around this function:

https://github.com/control-toolbox/CTBase.jl/blob/b3f2c3756c75c38087eb7f20891e6a9e85b2cba2/src/onepass.jl#L227

[frapac]

The variables' layout is defined here for the direct model: https://github.com/control-toolbox/CTDirect.jl/blob/main/src/problem.jl#L2

See also the following code, where the initial point is being set: https://github.com/control-toolbox/CTDirect.jl/blob/main/src/problem.jl#L205-L222

[ocots]

I think it is here: https://github.com/control-toolbox/CTDirect.jl/blob/main/src/problem.jl

@PierreMartinon Right?

[ocots]

Ici on peut voir comment récupérer le NLP : https://control-toolbox.org/OptimalControl.jl/stable/tutorial-nlp.html

Je sais pas si on peut ensuite manipuler les variables et autres pour rérdonner.

[PierreMartinon]

Hi everyone,

Up to the current release (0.9.6), the layout of variables is as follows:
[X0,X1..,XN, U0,U1,..,UN, V]
On main branch, and next release, the layout is:
[X0,U0, X1,U1, .. XN,UN, V]
This gives a more diagonal jacobian, with a small improvement in performance (less than 10% but it does seem a bit faster). The 0.10.0 release with the new layout should be out this week hopefully.

Although if I understand correctly, what you have is more an issue of columns / rows for the state variables. I don't think you can alter the ordering once the NLP is built, but you can certainly reorder the initial guess.

[frapac]

I am not surprised that the new ordering helps.
I think the ordering is of primal importance, and should be better documented in OC.

We should also have a closer look at what's going on inside the linear solver, and look at the KKT system once permuted by the linear solver.

It is well known that for optimal control problem, the KKT linear system is banded. However, this structure is not exploited when we use a generic linear solver. Long term, we should investigate if using a tailored solution (like Ricatti recursion) help in the solution.

For the reference, HPIPM and Fatrop are using Ricatti recursions internally.

[jbcaillau]

Thanks @frapac for the feedback. A few things:

• not clear to me which ordering would be better (if not the best); if sth is in order (conversely, to avoid), we can (and must) take care of it in CTDirect
• never heard of HPIPM (though we obviously know some of the people involved), Fatrop; just a bit of TrajectoryOptimisation / Altro
• these packages heavily rely on fast solvers for LQR subproblems, with MPC as a target (but not only, as this can be a general approach to solve a nonlinear problem iteratively): must have a look
• never heard either of Riccati recursion (see indeed Section 3.3 in this paper), but this further from my expertise; relevant, anyway
• yes, a tailored approach (including AD + sparsity treatment) is in order to exploit the structure and the appropriate linear algebra
• in addition, there is also the previously mentioned fact that all evaluations rely by design on a very limited set of functions (evaluated on a large number of grid points)

[PierreMartinon]

Note: we could have several variables layouts in CTDirect for testing purposes. The solution we return is the functional one anyway, so it does not depend on the NLP layout or even transcription formula. The NLP variables are only manipulated by getters/setters, so different layouts can be implemented in these functions without changing the rest of the code.

Things may be more complicated if we provide sparsity information at some point.

[jbcaillau]

@PierreMartinon what do you mean by layout? ordering of the NLP variables?

[jbcaillau]

@0Yassine0 can you please re-run the benchmark and let us know about speed-up regarding ADNLModels last update?
control-toolbox/CTDirect.jl#111 (comment)

[0Yassine0]

@jbcaillau
I pushed the re-run using the latest ADNLModels and OptimalControl updates ( see outputs folder).
I haven't noticed a huge improvement speed wise compared to the other version.

[jbcaillau]

thanks @0Yassine0

• just had a look at https://github.com/0Yassine0/COTS.jl/blob/main/Benchmark/outputs/model_Benchmark_file.pdf
• would be nice to build a summary of JuMP / OptimalControl.jl on more problems (than just 2)
• the array format is not so readable: better to gather performances on the same problem for both solves (on the same line... one line per problem and size, one column by solver...)
• I see a 3x factor in favour of JuMP on rocket, acceptable ; strictly less than 2x would be nice (and < 1x perfect, our target with Preparing for inplace NLP CTBase.jl#232 🙂)

[0Yassine0]

@jbcaillau
I added another version of results display with only the total time spend per solver applied on 4 models: https://github.com/0Yassine0/COTS.jl/blob/main/Benchmark/outputs/Model_Benchmark_TTonly_file.pdf
I hope it's better to evaluate performances.

[jbcaillau]

@jbcaillau I added another version of results display with only the total time spend per solver applied on 4 models: https://github.com/0Yassine0/COTS.jl/blob/main/Benchmark/outputs/Model_Benchmark_TTonly_file.pdf I hope it's better to evaluate performances.

very nice (and easier to read) 👍🏽 even better if you add tests for

• N = 500, 1000, 2000, 5000
• other models

[0Yassine0]

very nice (and easier to read) 👍🏽 even better if you add tests for

• N = 500, 1000, 2000, 5000

• other models

I added some problems and N values (500, 1000).
I couldn't add more because the execution took 30 min already to give results.
link

[jbcaillau]

very nice (and easier to read) 👍🏽 even better if you add tests for

• N = 500, 1000, 2000, 5000
• other models

I added some problems and N values (500, 1000). I couldn't add more because the execution took 30 min already to give results. link

👍🏽 @0Yassine0 could be nice to add the number of iterations (same solver = Ipopt for both solves). might explain some differences.

[0Yassine0]

👍🏽 @0Yassine0 could be nice to add the number of iterations (same solver = Ipopt for both solves). might explain some differences.

@jbcaillau
I added the number of iterations and changed the used models to the ones that gives the same results with both solvers ; here

[jbcaillau]

👍🏽 @0Yassine0 could be nice to add the number of iterations (same solver = Ipopt for both solves). might explain some differences.

@jbcaillau I added the number of iterations and changed the used models to the ones that gives the same results with both solvers ; here

thanks @0Yassine0 # of iterations confirm what I thought; does not explain why there are such differences, though. to be continued.

[jbcaillau]

hi @0Yassine0! assuming you're still around 🤞🏽, one last thing (apple mode off 🙂): could you please add the number of allocations / memory usage in the benchmark stats? will help to see what follows this upcoming update: control-toolbox/CTDirect.jl#188

[0Yassine0]

hi @0Yassine0! assuming you're still around 🤞🏽, one last thing (apple mode off 🙂): could you please add the number of allocations / memory usage in the benchmark stats? will help to see what follows this upcoming update: control-toolbox/CTDirect.jl#188

Yess @jbcaillau , I'm always available! The results are in the same file
I used @allocated for memory usage instead of running @benchmark from BenchmarkTools as usual.
If it's an issue I can change it.

[jbcaillau]

@0Yassine0 thanks! yes, @allocated that returns the number of bytes allocated (not the number of allocations) is fine.