Benchmarking

[ocots]

@jbcaillau @joseph-gergaud @PierreMartinon @frapac @tmigot @gdalle @amontoison @0Yassine0

I just want to tell to everyone that we are doing some benchmarking with @frapac and @0Yassine0. The starting point is this discussion: https://github.com/orgs/control-toolbox/discussions/47#discussioncomment-8554416

The idea is to compare some solvers first on the Goddard problem. Then, on the optimal control problems of COPSBenchmark.jl.

We want to play on the modeler (JuMP, ADNLPModel), the NLP solver (ipot, knitro, madlnp), the AD backend, the linear solver... We want to compare them in terms of sparsity, times...

Feel free to have a look at the repo COTS that will be forked into our control-toolbox organisation at the end and will be of course renamed :-)

[gdalle]

You're gonna have a much easier time playing with the AD backend once ADNLPModels switches to DifferentiationInterface. The work on that hasn't started, but DI is nearly ready for it.

[ocots]

Thanks for the info.

[ocots]

@PierreMartinon @jbcaillau @frapac @joseph-gergaud @0Yassine0

We should make a point about benchmarking.

Currently, we have some problems defined here:

• COTS.jl: we can find the model in JuMP and OptimalControl. Each problem comes with a notebook to compare both modelisation. Besides, we can find some benchmarking to compare linear solvers, nlp solvers, AD backends... The problems are rather complex and non-academic.
• CTProblems.jl: here is provided an OptimalControlModel with an OptimalControlSolution either analytical or numerical for each problem. The problems are small ones but cover the different possibilities: regular, bang, singular, boundary arcs. A list of the problems is given here with a way to extract some problems to make some studies.
• CTDirect.jl: here is provided an OptimalControlModel with the optimal objective value.

We should take some time to decide what to do about a repo of problems for benchmarking.

Remark. The package CTProblems.jl was made to have a list of problems to do some testing but only OptimalControl.jl uses it. There are only two tests on the double integrator and on goddard, both only on the objective value.

[PierreMartinon]

Hi @ocots
Indeed, some housekeeping may be in order :p

As for me, I stopped using CTProblems because it pulled the whole indirect part for the analytical solutions, which slowed down the tests considerably. This is why the problems in CTDirect do not include a complete solution but only the objective.

Maybe we will have to keep different sets of problems ?

• more complete and documented examples in CTproblems, to be used for general benchmarking
• internal problems for specific package testing, with only the relevant information

Also, there is the related question of applications. If I remember correctly, we were going for individual repos for applications, with a summary/link somewhere in OptimalControl. These applications would be typically linked with papers or conferences, and could provide material for both CTproblems and the internal testing problems.

Add to the next CT meeting ?

[ocots]

Do we have a meeting this week? Because i have sent an email to organise one.

[ocots]

I am not sure we need two sets for benchmarking and testing. For testing, it may be a small subset of the one for benchmarking. And maybe even for benchmarking, we should not do something to complex. Maybe in the doc, just the list of problems, no more.

For the applications, we have now individual repos. Each repo can have its own structure. We can mix math and simulations. As you said, I see these repos as extensions to articles or an article itself or a beginning. They could provide material for the set of problems for benchmarking.

[ocots]

Réunion mercredi 17/07 à 11:00. Ordre du jour.

J'aimerais que l'on discute rapidement des points suivants :

Point 1

• rappel des différents repo qui contiennent des problèmes
• choix de la structure du futur repo de benchmarking qui remplacera CTProblems.jl

Point 2

• retours sur le benchmarking fait par @0Yassine0 et @PierreMartinon : quelles sont les grosses différences entre JuMP et OptimalControl ?
• effort à faire pour la suite : prise en compte de la structure et calcul à la main des dérivées, cibler des amélioration des modeleurs, gpu...

[ocots]

Sur le point 1

On pourrait simplement faire un repo (nom à déterminer) :

• OptimalControlBenchmark
• OCBenchmark
• OptimalControlProblems
• OCProblems

qui contiendrait une liste de problèmes modélisés :

• JuMP avec règle des trapèzes
• OptimalControl
• ...

Puis pour récupérer un problème, par exemple celui de Goddard, on ferait simplement :

using OptimalControlBenchmark
ocp = Goddard() # ou Goddard(model=:OptimalControl), Goddard(model=:JuMP)

• Fait-on une documentation ?
• Ajoute-t-on des métadonnées à chaque problème ? Comme la valeur de l'objectif, la satisfaction des contraintes pour obtenir cet objectif, le nombre de variables d'état, de contrôle...

Remarque. Intégrons-nous la modélisation JuMP dans OptimalControl ? ADNLPModel $\Rightarrow$ JuMP ?

[amontoison]

Pourquoi ne pas les ajouter à OptimizationProblems.jl directement?
Ça donnera plus de visibilité à vos problèmes (à mon sens).

[ocots]

A discuter car vous avez bien des problèmes modélisés en JuMP et ADNLPModels ce qui est en soi comme ce que nous prévoyons. Mais nos problèmes de contrôle optimal seraient noyés dans la masse des problèmes d'optimisation. Je n'ai pas d'avis tranché.

[ocots]

Sur le point 2

Un effort est encore à faire sur le calcul des dérivées. Jacobienne et hessienne du lagrangien. Des efforts ont été fait du côté de ADNLPModels mais d'autres restent à faire notamment pour le calcul de la hessienne. Etant donnée la structure très particulière de notre problème, nous pourrions améliorer nos performances en utilisant cette information. En effet, nous rappelons sur toute la grille de discrétisation, les mêmes fonctions de dynamique et contraintes simplement en des valeurs différentes de leurs arguments.

• Est-ce que l'on essaie d'écrire nous-même les dérivées et d'optimiser ces calculs pour voir les effets sur les résultats ? Il y a la possibilité d'utiliser DifferentationInterface qui permet de "préparer" la dérivée puis de l'appliquer ce qui est plus efficace si on rappelle de nombreuses fois la dérivée de la même fonctionnelle sur des valeurs différentes.
• Pourquoi JuMP est-il plus performant alors qu'a priori on ne lui donne pas plus d'information ?
• Est-il envisageable, si ceci est crucial, que le modeleur lui-même gère cela ?
• Est-ce que l'on doit jouer plutôt sur la structure même de la hessienne en choisissant d'ordonner différemment les inconnues du vecteurs à optimiser pour le problème NLP ? Il me semble que non, que c'est déjà le bon choix que l'on fait. Dans ce cas, pourquoi JuMP est plus efficace, alors qu'a priori les variables sont ordonnées différemment.

Remarque. Il serait intéressant suite à ces études :

• de rédiger un article pour faire la synthèse de ce travail
• d'écrire un tutoriel du style "How to solve faster" qui explique comment on peut jouer sur les solveurs, les backends, etc, pour être plus efficace. Sans être catégorique, ni proposer quelque chose de systématique, mais plutôt des conseils et présenter les possibilités. Ce tutoriel pourrait ne pas faire tourner de code s'il faut des licences spéciales, etc.

[gdalle]

votre réunion demain c'est sur quel fuseau horaire ? je peux vous faire le topo DifferentiationInterface si vous voulez

[ocots]

C'est demain 11:00 sur le fuseau horaire de la France :-)

[gdalle]

si tu m'envoies un lien par email je peux faire une apparition sur Zoom

[ocots]

J'ai envoyé le lien Zoom par mail.

[jbcaillau]

https://inp-toulouse-fr.zoom.us/j/8963153648

[ocots]

Et pour info, OptimalControl.jl est sur le registre général depuis 1h environ !

[jbcaillau]

Etant donnée la structure très particulière de notre problème, nous pourrions améliorer nos performances en utilisant cette information. En effet, nous rappelons sur toute la grille de discrétisation, les mêmes fonctions de dynamique et contraintes simplement en des valeurs différentes de leurs arguments [...]

Pour compléter le message d' @ocots : nous devons absolument améliorer drastiquement les performances en utilisant deux avantages inhérents à la discrétisation directe d'un problème de contrôle :

• le caractère creux
• le potentiel de parallélisme SIMD

Un problème d'opti typique issu de ce type de discrétisation a la forme suivante :

$$ \begin{align} & g(x_0, x_N) \to \min,\\ & f(x_i, x_{i+1}, u_i, u_{i+1}) = 0,\quad i = 0,N-1,\\ & h_l \leq h(x_0, x_N) \leq h_u,\\ & c_l \leq c(x_i, u_i) \leq c_u,\quad i = 0,N,\\ & x_l \leq x_i \leq x_u,\quad i=0,N,\\ & u_l \leq u_i \leq u_u,\quad i=0,N \end{align} $$

où les inconnues sont $X = (x_0,\dots,x_N) \in \mathbf{R}^{(N+1)n}$ et $U = (u_0,\dots,u_N) \in \mathbf{R}^{(N+1)m}$ (les entiers $n$ et $m$ sont respectivement les dimensions de l'état et du contrôle, $N$ correspondant au nombre de points de discrétisation). Il s'agit donc d'évaluer seulement quatre fonctions, $f$, $c$, $g$ et $h$, sur un grand nombre de points, avec essentiellement des liaisons $i$, $i+1$.

On souhaite que ces fonctions puissent être du code Julia arbitraire, mais dans la plupart des cas il s'agit d'expressions algébriques (ou une classe à peine plus grande : $\log$, $\exp$, fonctions trigos...)

[PierreMartinon]

Hello a tous,

Je pense que nous pouvons en effet ameliorer le calcul des contraintes (et par suite des derivees).
Par exemple sur les allocations, j'ai essaye de reduire mais sans grand succes. Il faut dire que je ne suis pas un specialiste de Julia...

Pour le parallelisme, nous l'avions effectivement fait dans bocop, avec le calcul de toutes les dynamiques en une seule fois par exemple, parallelise de facon elementaire en openMP. Ceci pourrait sans doute etre etendu aux contraintes. A tester sur trapezes notamment qui s'y prete bien. J'avais commence, mais coince sur un probleme de diff auto en essayant de remplir le vecteur de dynamiques, on pourra regler ca ;-)

En passant, il devrait etre possible de faire des benchmarks sur l'ordonnancement des variables / contraintes, le code est relativement modulaire de ce point de vue. Le script de benchmark comporte une dizaine de problemes et devrait commencer a etre exploitable.

A tout a l'heure !

[ocots]

Mais moi la question que je me pose, c'est ok nous pouvons faire pleins d'amélioration en amont en exploitant la particularité d'un problème de contrôle optimal discrétisé, mais quand on modélise avec JuMP on ne fait pas ça et pourtant il semblerait que la différence est grande. Pourquoi ?

[gdalle]

Probablement parce que le solveur appelé par JuMP fait ça lui-même ? Quand on modélise avec un langage algébrique comme JuMP, les solveurs ont beaucoup plus d'informations sur la structure du problème, et sont capables de l'exploiter plus facilement qu'avec une simple fonction Julia.

[gdalle]

https://modernjuliaworkflows.github.io/optimizing/

[gdalle]

https://github.com/JuliaLinearAlgebra/BlockBandedMatrices.jl

[amontoison]

https://github.com/JuliaSmoothOptimizers/ManualNLPModels.jl

[gdalle]

Exemples de préparation pour l'autodiff :

• https://juliadiff.org/ForwardDiff.jl/stable/user/api/#Preallocating/Configuring-Work-Buffers
• https://docs.sciml.ai/FiniteDiff/dev/api/#FiniteDiff.GradientCache
• https://gdalle.github.io/DifferentiationInterface.jl/DifferentiationInterface/stable/operators/#Preparation

[amontoison]

Tutorial for ManualNLPModel:
https://jso.dev/tutorials/create-a-manual-model/

[amontoison]

@gdalle
JuliaSmoothOptimizers/ADNLPModels.jl#271

[amontoison]

@PierreMartinon https://github.com/amontoison/Workshop-GERAD

[amontoison]

@PierreMartinon @ocots @jbcaillau @0Yassine0
The AD backend of JuMP that can take into account redundant patterns in the Jacobian:
https://github.com/lanl-ansi/MathOptSymbolicAD.jl

[jbcaillau]

nice, thanks @amontoison (had missed this comment). this is a typical and crucial feature of discretised OCPs.

[ocots]

@0Yassine0 Have you tried it?

[frapac]

Pour reboucler sur notre discussion, Casadi s'est posé les mêmes questions sur comment gérer les Jacobiens par blocks dans la différentiation automatique. Voir notamment :
https://link.springer.com/article/10.1007/s12532-018-0139-4

[gdalle]

Apparemment la coloration symmétrique des matrices Hessiennes est un game-changer, on a un speedup x100 sur certaines instances du benchmark ici : JuliaSmoothOptimizers/ADNLPModels.jl#271 (comment)

[ocots]

Very impressive: JuliaSmoothOptimizers/ADNLPModels.jl#271 (comment)

[amontoison]

I released a new version of Ipopt_jll.jl tonight.
It contains the latest release of Ipopt and the dependencies (MUMPS and SPRAL).

• Ipopt 3.14.14 -> 3.14.16
• MUMPS 5.6.2 -> 5.7.3
• SPRAL 2024.1.18 -> 2024.5.8
• Add a compat entry for the upcoming release of HSL_jll.jl

The release v1.6.4 of Ipopt.jl will automatically install this new Ipopt_jll.jl.

[jbcaillau]

@amontoison any reason to expect improved performance?

[amontoison]

@jbcaillau I discussed with the MUMPS developers during SIAM LA24 and they added many features ans improved performances in the version 5.7.0.
The issue is that we use the sequential version in Ipopt so I don't know how it impacts us directly.

[jbcaillau]

@0Yassine0 if sth is available, can you please send a pointer towards the current synthesis of the comparisons between JuMP and OptimalControl?

[0Yassine0]

I have some results in the PDF files in this folder.

As discussed with @amontoison, there might be a problem either due to the order in which we set our variables or the initial guess, which I am currently investigating.

[jbcaillau]

All very nice, thanks @0Yassine0 👍🏽 Questions:

• given the comparison on nnz's, I guess you are using the last release of NLPModels (0.8.3)?
  
• is there a file merging and comparing JuMP and OptimalControl (nnz, total time, at least) on more use cases than the one below?
  

[0Yassine0]

• Yes, I'm using the latest release of every package.
• I just modified my function and now we can give a list of models, I ran it for Hang glider problem and Goddard rocket problem (see model_Benchmark_file).
  

[jbcaillau]

very nice. good to have as many comparisons as possible 🤞🏾

[jbcaillau]

@frapac am I missing sth obvious here to reproduce the examples in the doc? exanauts/ExaModels.jl#104

@0Yassine0 have you had a look at these ExaModels examples?

[jbcaillau]

ongoing: exanauts/ExaModels.jl#105 and exanauts/ExaModels.jl#106

[amontoison]

We fixed the issue with the new release 0.21.3 of NLPModels.jl.

[jbcaillau]

@amontoison oh, thanks. wasn’t aware JuliaSmoothOptimizers/NLPModels.jl#475

[jbcaillau]

@rveltz AFAIR you've made some tests on GPUs on Inria cluster (nef?) Can you please point me towards a worked out example / framework that I could start from to do some benchmarking?

[rveltz]

HiI am away till next week, i ll send you a link then. If it is more pressing contact Aurora, I got the trick from her.Romain -------- Message d'origine --------De : Jean-Baptiste Caillau ***@***.***> Date : 06/08/2024 09:49 (GMT-07:00) À : "control-toolbox/control-toolbox.github.io" ***@***.***> Cc : Romain Veltz ***@***.***>, Mention ***@***.***> Objet : Re: [control-toolbox/control-toolbox.github.io] Benchmarking (Discussion #56) @rveltz AFAIR you've made some tests on GPUs on Inria cluster (nef?) Can you please point me towards a worked out example / framework that I could start from to do some benchmarking? —Reply to this email directly, view it on GitHub, or unsubscribe.You are receiving this because you were mentioned.Message ID: ***@***.***>

[jbcaillau]

thanks @rveltz currently looking at nef FAQ Will email aurora (she's actually named after a supercomputer - should do it 🤞🏾)