305 solve documentation

docs/Project.toml

@@ -24,10 +24,10 @@ Suppressor = "fd094767-a336-5f1f-9728-57cf17d0bbfb"
 [compat]
 ADNLPModels = "0.8"
 CTBase = "0.16"
-CTDirect = "0.15"
 CTFlows = "0.8"
 CTModels = "0.6"
-CTParser = "0.5"
+CTParser = "0.6"
+CTDirect = "0.16"
 CommonSolve = "0.2"
 DataFrames = "1"
 Documenter = "1.8"

docs/make.jl

@@ -111,7 +111,7 @@ cp("./docs/Project.toml", "./docs/src/assets/Project.toml"; force=true)
 repo_url = "github.com/control-toolbox/OptimalControl.jl"
 
 makedocs(;
-    draft=true, # if draft is true, then the julia code from .md is not executed # debug
+    draft=false, # if draft is true, then the julia code from .md is not executed # debug
     # to disable the draft mode in a specific markdown file, use the following:
     # ```@meta
     # Draft = false

docs/src/manual-abstract.md

@@ -22,7 +22,7 @@ A variable (only one is allowed) is a finite dimensional vector or reals that wi
 end
 ```
 
-!!! caveat
+!!! warning
     Note that the full code of the definition above is not provided (hence the `...`) The same is true for most examples below (only those without `...` are indeed complete). Also note that problem definitions must at least include definitions for time, state, control, and dynamics.
 
 
@@ -192,7 +192,7 @@ F₁(x) = [0, 1]
 !!! note
     The vector fields `F₀` and `F₁` can be defined afterwards, as they only need to be available when the dynamics will be evaluated.
 
-Currently, it is not possible to declare the dynamics component after component, but a simple workaround is to use *aliases* (check the relevant [aliases](@ref manual-abstract-aliases) section below):
+While it is also possible to declare the dynamics component after component (see below), one may equivalently use *aliases* (check the relevant [aliases](@ref manual-abstract-aliases) section below):
 
 ```julia
 @def damped_integrator begin
@@ -207,7 +207,30 @@ Currently, it is not possible to declare the dynamics component after component,
 end
 ```
 
-## Constraints
+## [Dynamics (coordinatewise)](@id manual-abstract-dynamics-coord)
+
+```julia
+:( ∂($x[$i])($t) == $e1 ) 
+```
+
+The dynamics can also be declared coordinate by coordinate. The previous example can be written as
+
+```julia
+@def damped_integrator begin
+    tf ∈ R, variable
+    t ∈ [0, tf], time
+    x = (q, v) ∈ R², state
+    u ∈ R, control
+    ∂(q)(t) == v(t)
+    ∂(v)(t) == u(t) - c(t)
+    ...
+end
+```
+
+!!! warning
+    Declaring the dynamics coordinate by coordinate is **compulsory** when solving with the option `:exa` to rely on the ExaModels modeller (check the [solve section](@ref manual-solve)), for instance to [solve on GPU](@ref manual-solve-gpu).
+
+## [Constraints](@id manual-abstract-constraints)
 
 ```julia
 :( $e1 == $e2        ) 
@@ -218,7 +241,7 @@ end
 ```
 
 Admissible constraints can be
-- five types: boundary, control, state, mixed, variable,
+- of five types: boundary, variable, control, state, mixed (the last three ones are *path* constraints, that is constraints evaluated all times)
 - linear (ranges) or nonlinear (not ranges),
 - equalities or (one or two-sided) inequalities.
 
@@ -265,7 +288,7 @@ end
 end
 ```
 
-!!! caveat
+!!! warning
     Write either `u(t)^2` or `(u^2)(t)`, not `u^2(t)` since in Julia the latter means `u^(2t)`. Moreover,
     in the case of equalities or of one-sided inequalities, the control and / or the state must belong to the *left-hand side*. The following will error:
 
@@ -286,7 +309,7 @@ using OptimalControl
 end
 ```
 
-!!! caveat
+!!! warning
      Constraint bounds must be *effective*, that is must not depend on a variable. For instance, instead of
 ```julia
 o = @def begin
@@ -318,6 +341,9 @@ o = @def begin
 end
 ```
 
+!!! warning
+    When solving with the option `:exa` to rely on the ExaModels modeller (check the [solve section](@ref manual-solve)), for instance to [solve on GPU](@ref manual-solve-gpu), it is **compulsory** that *nonlinear* constraints (not ranges) are *scalar*, whatever the type (boundary, variable, controle, state, mixed).
+
 ## [Mayer cost](@id manual-abstract-mayer)
 
 ```julia                                      
@@ -423,7 +449,7 @@ Quite readily, Mayer and Lagrange costs can be combined into general Bolza costs
 end
 ```
 
-!!! caveat
+!!! warning
     The expression must be the sum of two terms (plus, possibly, a scalar factor before the integral), not *more*, so mind the parentheses. For instance, the following errors:
 
 ```julia
@@ -469,7 +495,7 @@ The single `=` symbol is used to define not a constraint but an alias, that is a
 end
 ```
 
-!!! caveat
+!!! warning
     Such aliases do *not* define any additional function and are just replaced textually by the parser. In particular, they cannot be used outside the `@def` `begin ... end` block.
 
 !!! hint
@@ -487,7 +513,7 @@ end
 end true;
 ```
 
-!!! caveat
+!!! warning
     The dynamics of an OCP is indeed a particular constraint, be careful to use `==` and not a single `=` that would try to define an alias:
 
 ```@repl main-repl
@@ -528,4 +554,4 @@ end
 
 ## Known issues
 
-- [Constants and (reverse over forward) AD](https://github.com/control-toolbox/OptimalControl.jl/issues/481)
+- [Reverse over forward AD issues with ADNLP](https://github.com/control-toolbox/OptimalControl.jl/issues/481)

docs/src/manual-solve-gpu.md

@@ -4,10 +4,78 @@
 CollapsedDocStrings = false
 ```
 
-In this manual, we explain how to use the [`solve`](@ref) function from [OptimalControl.jl](https://control-toolbox.org/OptimalControl.jl) on GPU. We rely on [MadNLP](https://github.com/MadNLP/MadNLP.jl) and currently only provide support for NVIDIA thanks to [CUDA.jl](https://github.com/JuliaGPU/CUDA.jl).
+In this manual, we explain how to use the [`solve`](@ref) function from [OptimalControl.jl](https://control-toolbox.org/OptimalControl.jl) on GPU. We rely on [ExaModels.jl](https://exanauts.github.io/ExaModels.jl/stable) and [MadNLPGPU.jl](https://github.com/MadNLP/MadNLP.jl) and currently only provide support for NVIDIA thanks to [CUDA.jl](https://github.com/JuliaGPU/CUDA.jl). Consider the following simple Lagrange optimal control problem:
 
-```@docs; canonical=false
-solve(::CTModels.Model, ::Symbol...)
+ ```julia
+using OptimalControl
+using MadNLPGPU
+using CUDA
+
+ocp = @def begin
+    t ∈ [0, 1], time
+    x ∈ R², state
+    u ∈ R, control
+    v ∈ R, variable
+    x(0) == [0, 1]
+    x(1) == [0, -1]
+    ∂(x₁)(t) == x₂(t)
+    ∂(x₂)(t) == u(t)
+    0 ≤ x₁(t) + v^2 ≤ 1.1
+    -10 ≤ u(t) ≤ 10
+    1 ≤ v ≤ 2
+    ∫(u(t)^2 + v) → min
+end
 ```
 
-TBC
+!!! note
+    We have used MadNLPGPU instead of MadNLP, that is able to solve on GPU (leveraging [CUDSS.jl](https://github.com/exanauts/CUDSS.jl)) optimisation problems modelled with ExaModels.jl. As a direct transcription towards an `ExaModels.ExaModel` is performed, there are limitations on the syntax:  
+    - dynamics must be declared coordinate by coordinate (not globally as a vector valued expression)
+    - nonlinear constraints (boundary, variable, control, state, mixed ones, see [Constraints](@ref manual-abstract-constraints) must also be scalar expressions (linear constraints *aka.* ranges, on the other hand, can be vectors)
+    - all expressions must only involve algebraic operations that are known to ExaModels (check the [documentation](https://exanauts.github.io/ExaModels.jl/stable)), although one can provide additional user defined functions through *registration* (check [ExaModels API](https://exanauts.github.io/ExaModels.jl/stable/core/#ExaModels.@register_univariate-Tuple%7BAny,%2520Any,%2520Any%7D)) 
+
+Computation on GPU is currently only tested with CUDA, and the associated backend must be passed to ExaModels as is done below (also note the `:exa` keyword to indicate the modeller, and `:madnlp` for the solver):
+
+```julia
+julia> sol = solve(ocp, :exa, :madnlp; exa_backend=CUDABackend())
+▫ This is OptimalControl version v1.1.0 running with: direct, exa, madnlp.
+
+▫ The optimal control problem is solved with CTDirect version v0.16.0.
+
+   ┌─ The NLP is modelled with ExaModels and solved with MadNLP.
+   │
+   ├─ Number of time steps⋅: 250
+   └─ Discretisation scheme: trapeze
+
+▫ This is MadNLP version v0.8.7, running with cuDSS v0.4.0
+
+Number of nonzeros in constraint Jacobian............:     2506
+Number of nonzeros in Lagrangian Hessian.............:     2006
+
+Total number of variables............................:      754
+                     variables with only lower bounds:        0
+                variables with lower and upper bounds:      252
+                     variables with only upper bounds:        0
+Total number of equality constraints.................:      504
+Total number of inequality constraints...............:      251
+        inequality constraints with only lower bounds:        0
+   inequality constraints with lower and upper bounds:      251
+        inequality constraints with only upper bounds:        0
+
+iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
+   0  1.0200000e+00 1.10e+00 1.00e+00  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
+  ...
+  26  9.8902986e+00 2.22e-16 7.11e-15  -9.0 1.32e-04    -  1.00e+00 1.00e+00h  1
+
+Number of Iterations....: 26
+
+                                   (scaled)                 (unscaled)
+Objective...............:   9.8902986337530514e+00    9.8902986337530514e+00
+Dual infeasibility......:   7.1054273576010019e-15    7.1054273576010019e-15
+Constraint violation....:   2.2204460492503131e-16    2.2204460492503131e-16
+Complementarity.........:   4.8363494304578671e-09    4.8363494304578671e-09
+Overall NLP error.......:   4.8363494304578671e-09    4.8363494304578671e-09
+
+...
+
+EXIT: Optimal Solution Found (tol = 1.0e-08).
+```

docs/src/manual-solve.md

@@ -48,7 +48,7 @@ This is because the default method uses a direct approach, which transforms the
 \text{minimize}\quad F(y), \quad\text{subject to the constraints}\quad g(y) \le 0, \quad h(y) = 0. 
 ```
 
-!!! caveat
+!!! warning
 
     Calling `solve` without loading a NLP solver package first will notify the user:
 
@@ -89,7 +89,7 @@ solve(ocp, :direct, :adnlp, :ipopt)
 
 !!! warning
 
-    The dynamics must be defined coordinatewise to use ExaModels.jl (`:exa`).
+    Dynamics and nonlinear constraints must be defined coordinatewise to use ExaModels.jl (`:exa`). Check [the problem definition section](@ref manual-abstract-syntax) for more information.
 
 For instance, let us try MadNLP solver with ExaModel modeller.