Optimisers

[jbcaillau]

@ocots @frapac @joseph-gergaud @PierreMartinon

Bonjour Olivier, Jean-Baptiste,

et merci pour ce retour ! Je pense qu'on se complémente bien en effet : nous sommes de notre coté justement intéressé par avoir une base de benchmark.

Parfait

Pour la suite, je pense que le plus simple est que je prenne en main OptimalControl.jl et essaie d'utiliser MadNLP pour > résoudre les problèmes, pour voir. À moyen terme, je me demande si nous pourrions essayer de monter un stage en > commun pour essayer d'implémenter un solveur KKT dédié aux problèmes de controle optimal, ala Stephen Wright.

Oui : OptimalControl.jl peut être simplement utilisé comme une interface de modélisation ; que ce soit avec la syntaxe "abstraite" ¹ ou "fonctionnelle" ², le pb de contrôle optimal est traduit

• tout d'abord en une liste de contraintes NLP sous la forme de contraintes linéaires élémentaires (boites sur un range d'indices) ou non-linéaires (boites sur des fonctions), voir nlp_constraints
• ensuite sous la forme d'un ADNLPModel ; cette traduction est pour l'instant interne à la fonction CTDirect.solve mais a vocation à être extraite et généralisée à d'autres modeleurs / solvers. Voir cette issue.

Volontiers pour dév. commun, par exemple en s'appuyant sur un stage. Calons un point ensemble pour préciser.

Footnotes

1. voir cet exemple ↩

2. voir cet autre exemple ↩

[ocots]

Dans un premier temps, on aimerait avoir le problème de Goddard et d'autres cas d'usage sous les différentes formes¹ données par le diagramme ci-dessous.

Légende. Quand je mets .jl dans le diagramme c'est que l'on utilise le package mais je ne détaille pas comment. Sinon, je parle directement du type. OCP c'est le problème de contrôle optimal donné par sa formulation mathématique.

stateDiagram-v2
	OCP --> JuMP.jl
	OCP --> ExaModels.jl
	OCP --> OptimalControl.jl
	OptimalControl.jl --> ADNLPModel
	JuMP.jl --> Ipopt.jl
	ExaModels.jl --> ExaModel
	ADNLPModel --> AbstractNLPModel 
        ExaModel --> AbstractNLPModel 
        AbstractNLPModel --> NLPModelsIpopt.jl
	AbstractNLPModel --> MadNLP.jl
        NLPModelsIpopt.jl --> Ipopt.jl

Loading

• ADNLPModels.jl : depuis OptimalControl.jl, on définit un ADNLPModel qui est un sous-type de AbstractNLPModel.
• ExaModels.jl : le type ExaModel est un sous-type de AbstractNLPModel. Pour qu'un solveur puisse gérer un nouveau type, sous-type de AbstractNLPModel, il suffit de surcharger quelques fonctions de NLPModels.jl qui retourne la fonction, son gradient... C'est le multiple dispatch qui gère le truc. En effet, puisqu'il n'y a pas d'héritage, on passe par le multiple dispatch.
• Ipopt.jl
• JuMP.jl
• MadNLP.jl : il peut résoudre un problème défini par un AbstractNLPModel donc un ExaModel ou un ADNLPModel par exemple. Remarque. Il peut aussi se brancher à JuMP et donc résoudre un problème défini par JuMP.
• NLPModels.jl : défini AbstractNLPModel.
• NLPModelsIpopt.jl : ce n'est pas pour modéliser mais c'est un wrapper d'ipopt capable de résoudre un problème défini par un NLPModel.
• OptimalControl.jl

Footnotes

1. vous pouvez faire des modifications sur le diagramme. ↩

[ocots]

Voulons-nous :

flowchart TB
B(OCP) --> OptimalControl.jl --> Packages
Types --> AbstractNLPModel
AbstractNLPModel --> Wrappers
subgraph Modelisation
direction TB
subgraph Packages
ADNLPModels.jl
ExaModels.jl
JuMP.jl
end
subgraph Types
ADNLPModel
ExaModel
NLPModelsJuMP
end
ADNLPModels.jl --> ADNLPModel
ExaModels.jl --> ExaModel
JuMP.jl --> NLPModelsJuMP
AbstractNLPModel
end
subgraph Resolution
subgraph Wrappers
NLPModelsIpopt.jl
NLPModelsKnitro.jl
end
subgraph Solvers
Ipopt.jl
KNITRO.jl
MadNLP.jl
end
NLPModelsIpopt.jl --> Ipopt.jl
NLPModelsKnitro.jl --> KNITRO.jl
end
AbstractNLPModel --> MadNLP.jl
JuMP.jl --> Solvers

Loading

• KNITRO.jl : Remarque. KNITRO peut aussi se brancher sur JuMP bien entendu.
• NLPModelsKnitro.jl : c'est un wrapper à KNITRO pour résoudre un problème défini par un AbstractNLPModel.
• NLPModelsJuMP

[ocots]

Ce serait bien d'avoir un benchmark de ces méthodes sur une liste de problèmes de contrôle optimal.

[jbcaillau]

@ocots Yes! Going towards a list of control problems (not their discretisation as math. programs or so) using a high level description (e.g. OptimalControl abstract syntax 🙂) would be a great outcome of CTProblems

[ocots]

We need to define an AbstractOCPModel or AbstractOptimalControlModel and methods on it, used by our solvers and functors such as CTDirect. If everything is written taking as argument an AbstractOCPModel, then it becomes easier to implement different OCP modelisations. Then, from our AbstractOCPModel we get an NLP thanks for instance to CTDirect:

using OptimalControl
ocp = Model()      # ocp is a OCPModel <: AbstractOCPModel
nlp = Direct(ocp, :adnlp) # nlp is a AbstractNLPModel

or

ocp = Model()
ocp = Model(:jump) # we define the model à la JuMP?
ocp = Model(:exa) # ...

and

nlp = Direct(ocp, :adnlp)
nlp = Direct(ocp, :jump)
nlp = Direct(ocp, :exa)

[ocots]

Maybe we need also a type: AbstractDOCPModel for abstract discretized optimal control problem.

using OptimalControl
ocp = Model()
docp = Direct(ocp, steps=100, scheme=:trapeze)
nlp = NLP(docp, :adnlp) # or nlp = NLP(docp, :exa)
typeof(nlp) <: AbstractNLPModel # return true

and we can solve nlp by any solver we want. Other possibilities:

using OptimalControl
ocp = Model()
docp = Direct(ocp, steps=100, scheme=:trapeze)
solve(docp, :adnlp, :ipopt)

or

using OptimalControl
ocp = Model()
solve(ocp, :direct, :adnlp, :ipopt, steps=100, scheme=:trapeze)

[PierreMartinon]

This could be a good idea. As a side note, when writing the transcription function (currently takes an ocp ans returns the corresponding nlp) I encountered the small problem of the initial guess: ADnlp typically expects the starting point to be given, while for me this is not strictly part of the discretized problem but more a piece of information passed to the solve() function.

Back to your first example block, we could for instance decide that the initial guess would be normally passed to the NLP() function, not the Direct() function. More generally, I feel it should be possible to call solve() with a specific initial point, so we would need a function that would set/update the initial point in the NLP problem anyway. Worst case scenario, we'll need to create a new nlp for each initial guess.

[ocots]

On remarque deux choses :

• Pour NLPModels ils font eux-mêmes les wrappers aux solveurs là où pour JuMP, ce sont les solveurs qui font en sorte d'être compatibles en passant par MathOptInterface.
• En ce qui concerne MadNLP, ils se sont rendus compatibles à la fois à JuMP et à NLPModels.

[jbcaillau]

@frapac @ocots @PierreMartinon @joseph-gergaud

• MadNLP.jl: perf. vs. Ipopt : approx. x1.2, interface with NLPModels.jl
• ExaModels.jl: no user defined functions (as JuMP), but allows vector constraints (JuMP doesn't)
• for the record: ADNLPModels.jl can use several backends (ForwardDiff.jl, in practice) for AD of full Julia code
• motivation to allow full Julia expressions: avoid code duplication between several DSLs (the common base should be Julia not yet another DSL)
• done in OptimalControl.jl, including no "world age pb" (check this test)
• idea: allow ExaModels expressions (it is Julia code, indeed) inside OptimalControl @def? (Implies a specific post-treatment)
• @frapac @ocots internship (M1/M2 level): use case (Goddard and more) to compare JuMP and MadNLP
• related topic: test with Knitro.jl (cf. trust-region vs. interior point)
• conferences mentioned by @gergaud: SMAI-MODE (deadline 18/12/2023, no MS), FGS'2024 (MS deadline 14/02/2024)
• @frapac please send a pointer towards the power grid optim. benchmark (Los Alamos) you indicated

[frapac]

the benchmark of the different AD backends on the ACOPF problems can be found here: https://discourse.julialang.org/t/ac-optimal-power-flow-in-various-nonlinear-optimization-frameworks/78486/59

We should keep in mind that JuMP remains slower than AMPL by a factor of two. In that regards, ADNLPModels.jl is not tractable for us (yet).

[frapac]

has someone tried FastDifferentiation.jl? The original package is around since the 2000s, and has been ported in Julia recently. I would be curious to know the performance of that other AD backend.

[jbcaillau]

@frapac not heard of it, having a look (-> added in this issue)

has someone tried FastDifferentiation.jl? The original package is around since the 2000s, and has been ported in Julia recently. I would be curious to know the performance of that other AD backend.

[jbcaillau]

PS. Congrats @frapac and co-authors !!!

https://www.linkedin.com/posts/fran%C3%A7ois-pacaud-99b01487_github-madnlpmadnlpjl-a-solver-for-nonlinear-activity-7120440464666701825-_evI?utm_source=share&utm_medium=member_ios

[frapac]

thank you :-)

I am glad MadNLP and ExaModels got some recognition from COIN-OR. I hope this will give some momentum to the MadNLP ecosystem!

[ocots]

Ideas for an internship:

In a first place it could be great to compare some pairs of (model, solve) for some case studies: Goddard... See the list of problems here.

The list of pairs of models and solvers (to be completed, updated) :

Model	Solver
JumP	Ipopt
JumP	KNITRO
JumP	MadNLP
ADNLPModel	Ipopt
ADNLPModel	KNITRO
ADNLPModel	MadNLP
ExaModel	Ipopt
ExaModel	KNITRO
ExaModel	MadNLP
Casadi	Ipopt
Casadi	KNITRO
Bocop	Ipopt

Then, it could be great to define a process for benchmarking with online results and with the possibility for users to make their own benchmarks.

[ocots]

We need to define in which terms we make the comparisons. We should also compare the cost of a call to the derivatives.

[ocots]

WORHP ?

[ocots]

Pour le choix des problèmes, voir benchmark COPS et son implémentation JumP, cf. package COPSBenchmark.jl.

[ocots]

Les différents problèmes intéressants : Goddard, un problème lisse, du chattering, de nombreuses commutations (transfert orbital), contraintes d'état. Contraintes MPEC, problèmes hybrides (au delà du stage).

[jbcaillau]

@ocots @PierreMartinon @joseph-gergaud check this post on JuMP and related things (thanks @tmigot for pointing to ct). Most notably, function tracing now available in JuMP: @tmigot @amontoison does this entail full AD support (though not sparse, I guess) for the expressions generated by overloading on JuMP abstract expression type?

[amontoison]

I didn't find time yet @jbcaillau to check the new NLP API of JuMP.
Maybe this summer for the JuMP-dev at Montréal.

[ocots]

Stage M1 2024 :

Livrables :

• liste des problèmes tests dans un dépôt Github :
  • Priorité 1 : définis sous JumP et OptimalControl.jl
  • Priorité 2 : Casadi
  • Priorité 3 : Bocop (voir avec Pierre)
• performance profile et tableau instance vs comparaison des outils (temps dans la diff auto, temps total, nombre itérations, sparse pattern) sur les exemples tests pour les paires (Model, Solver) suivantes :
  • Task 1 : (JuMP, Ipopt, JuMPDefault backend, _linear_solver), _linear_solver in {MUMPS, HSLMA57}
  • Task 2 : (JuMP, Ipopt, _backend, best_linear_solver), _backend in {JuMPDefault, SymbolicAD, ExaModel}
  • Task 3 : (JuMP, _solver, best_backend, best_linear_solver), _solver in {Ipopt, Knitro_IPM, Knitro_SQP, MadNLP}
  • Task 4 : (ADNLPModel, Ipopt, best_linear_solver) : via OptimalControl.jl cf. remarque
  • Task 5 : OptimalControl.jl -> JuMP (cf. point plus haut)

Remarque : ADNLPModel est utilisé via OptimalControl.jl. Pour le moment OptimalControl.jl ne fournit que ADNLPModel + Ipopt.

Perspectives / Bonus :

• ADNLPModel avec les solveurs Knitro et MadNLP si intéressant
• Voir Casadi.jl
• Attention aux schémas de discrétisation, pour le moment seul les trapèzes sont disponibles dans la control-toolbox.

[jbcaillau]

@ocots @frapac trouvé un stagiaire M1 / M2 sur cette action opti ?

[frapac]

Nous en avions un mais il nous a fait faux bond. Aucun autre candidat pour le moment.

[ocots]

On devrait avoir quelqu'un ! Je vous tiens au courant.

[jbcaillau]

@amontoison @tmigot rappel des infos sur JuMP discutées dans ce post plus haut avec

• function tracing (essentiellement surcharge des opérateurs classiques, puis ReverseAD sur les expressions générées dans le DSL JuMP)
• user defined operators (and falling back to ForwardDiff.jl)

[amontoison]

Merci @jbcaillau!
Je comprends comment ça fonctionne maintenant.
Ça ne doit pas être bien compliqué à interfacer dans NLPModelsJuMP.jl.

D'un point de vue utilisateur, il faut néanmoins définir des fonctions qui peuvent prendre un Vector{VariableRef} comme argument.
C'est ce qu'on a fait avec Tangi dans ADNLPModels.jl.
La différence c'est qu'on utilise Symbolics.jl et qu'on peut définir f(x::Vector{Float64}).
Symbolics.jl arrive à nous donner le "sparsity pattern" du Jacobien / Hessien.

Je pense que JuMP arrive à faire une analyse plus en détails et isoler l'arbre d'expression de f.
Je vais en parler à Oscar à l'occasion.

[jbcaillau]

Le truc c'est de pouvoir utiliser des fonctions définies en julia standard, et non pas nécessairement des expressions bâties sur les briques de base de JuMP. Si leur surcharge (aka "function tracing") fait le job de façon raisonnablement large c'est OK. À voir.

[ocots]

@frapac For the benchmark we should have a look to:

• https://jso.dev/OptimizationProblems.jl/stable/
• https://github.com/JuliaSmoothOptimizers/SolverBenchmark.jl

[ocots]

AD benchmark : JuliaSmoothOptimizers/ADNLPModels.jl#241

[amontoison]

L'équipe de contrôle optimal, vous utilisez quoi comme backend pour le Hessian creux dans ADNLPModels.jl?
Je viens d'y penser ce soir en répondant à un autre message, vous utilisez SparseADHessian ou SparseReverseADHessian?
Le second fait du "forward over reverse" et devrait être optimal si ReverseDiff arrive à manger votre fonction.
Par défaut c'est l'autre backend SparseADHessian qui fait du "forward over forward". Il est plus robuste parce que ForwardDiff.jl mange tout.

[ocots]

@PierreMartinon Une idée ?

[PierreMartinon]

Pour le moment on utilise l'option par defaut, donc SparseADHessian a priori, mais on pourrait tester le reverse oui !
On pourrait passer le choix de backend a solve/directTranscription (qui genere le ADNLPModel) et aussi avoir une fonction qui modifie un probleme existant du genre setADbackend(docp, ...), comme pour l'initialization.

Edit: ah en fait on utilise backend=:optimized

[jbcaillau]

merci @amontoison sur ce sujet, voir ici :

• il y a bien un set_AD_backend, voir ici
• l'idée étant par ex. de pouvoir faire du reverse over forward 🤞🏽

@PierreMartinon just added this as an issue

[ocots]

@frapac @0Yassine0 : Hessienne.