From d33a0a480c529a5bedb856622efbecad168274ca Mon Sep 17 00:00:00 2001
From: Jean-Baptiste Caillau <jean-baptiste.caillau@univ-cotedazur.fr>
Date: Fri, 2 Aug 2024 00:44:48 +0200
Subject: [PATCH 01/17] first commit

---
 docs/src/juliacon2024.md | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/docs/src/juliacon2024.md b/docs/src/juliacon2024.md
index f8804275..7907aa53 100644
--- a/docs/src/juliacon2024.md
+++ b/docs/src/juliacon2024.md
@@ -32,8 +32,7 @@ plus boundary, control and state constraints
 
 - [Basic example](tutorial-basic-example.html)
 - [Goddard problem](tutorial-goddard.html)
-- [Orbit transfer](application-orbit.html)
-- [Solar sailing](application-sail.html)
+- [Orbit transfer](http://control-toolbox.org/kepler/stable)
 
 ## Wrap up
 

From 0e7210e8d94337d661f78da23874a08cb4dff3b8 Mon Sep 17 00:00:00 2001
From: Jean-Baptiste Caillau <jean-baptiste.caillau@univ-cotedazur.fr>
Date: Fri, 2 Aug 2024 00:45:35 +0200
Subject: [PATCH 02/17] first commit

---
 docs/make.jl | 1 +
 1 file changed, 1 insertion(+)

diff --git a/docs/make.jl b/docs/make.jl
index 74850b8f..d0188e24 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -40,6 +40,7 @@ makedocs(
             "tutorial-lqr-basic.md",
             "tutorial-iss.md",
             "tutorial-goddard.md",
+            "tutorial-abstract.md",
         ],
         "API" => [
             "api-optimalcontrol.md",

From 4eab766d378de8249c086f76b09ab071a2a4b2b2 Mon Sep 17 00:00:00 2001
From: Jean-Baptiste Caillau <jean-baptiste.caillau@univ-cotedazur.fr>
Date: Fri, 2 Aug 2024 00:58:20 +0200
Subject: [PATCH 03/17] doc update

---
 docs/src/tutorial-abstract.md | 38 +++++++++++++++++++++++++++++++++++
 1 file changed, 38 insertions(+)
 create mode 100644 docs/src/tutorial-abstract.md

diff --git a/docs/src/tutorial-abstract.md b/docs/src/tutorial-abstract.md
new file mode 100644
index 00000000..3b840075
--- /dev/null
+++ b/docs/src/tutorial-abstract.md
@@ -0,0 +1,38 @@
+# [Abstract syntax](@id abstract)
+
+Todo:
+- give the full grammar
+- thanks MLStyle
+- examples for most features
+- link to example + API for functional syntax
+
+Let us consider a wagon moving along a rail, whom acceleration can be controlled by a force $u$.
+We denote by $x = (x_1, x_2)$ the state of the wagon, that is its position $x_1$ and its velocity $x_2$.
+
+```@raw html
+<img src="./assets/chariot.png" style="display: block; margin: 0 auto 20px auto;" width="300px">
+```
+
+We assume that the mass is constant and unitary and that there is no friction. The dynamics we consider is given by
+
+```math
+    \dot x_1(t) = x_2(t), \quad \dot x_2(t) = u(t),\quad u(t) \in \R,
+```
+
+which is simply the [double integrator](https://en.wikipedia.org/w/index.php?title=Double_integrator&oldid=1071399674) system.
+Les us consider a transfer starting at time $t_0 = 0$ and ending at time $t_f = 1$, for which we want to minimise the transfer energy
+
+```math
+    \frac{1}{2}\int_{0}^{1} u^2(t) \, \mathrm{d}t
+```
+
+starting from the condition $x(0) = (-1, 0)$ and with the goal to reach the target $x(1) = (0, 0)$.
+
+First, we need to import the `OptimalControl.jl` package to define the optimal control problem and `NLPModelsIpopt.jl` to solve it. 
+We also need to import the `Plots.jl` package to plot the solution.
+
+```@example main
+using OptimalControl
+using NLPModelsIpopt
+using Plots
+```

From 78115e92032d32d1420a409400015bba3c98fde1 Mon Sep 17 00:00:00 2001
From: Jean-Baptiste Caillau <jean-baptiste.caillau@univ-cotedazur.fr>
Date: Fri, 2 Aug 2024 19:08:53 +0200
Subject: [PATCH 04/17] docs

---
 docs/make.jl                                  |   2 +-
 docs/src/tutorial-abstract.md                 | 166 +++++++++++++++---
 docs/src/tutorial-goddard.md                  |   2 +-
 .../{tutorial-nlp.md => tutorial-nlp.md.foo}  |   0
 4 files changed, 146 insertions(+), 24 deletions(-)
 rename docs/src/{tutorial-nlp.md => tutorial-nlp.md.foo} (100%)

diff --git a/docs/make.jl b/docs/make.jl
index d0188e24..1a4bbb15 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -35,7 +35,7 @@ makedocs(
             "tutorial-double-integrator.md",
             "tutorial-initial-guess.md",
             "tutorial-continuation.md",
-            "tutorial-nlp.md",
+            #"tutorial-nlp.md",
             "tutorial-plot.md",
             "tutorial-lqr-basic.md",
             "tutorial-iss.md",
diff --git a/docs/src/tutorial-abstract.md b/docs/src/tutorial-abstract.md
index 3b840075..9a1bfbe8 100644
--- a/docs/src/tutorial-abstract.md
+++ b/docs/src/tutorial-abstract.md
@@ -1,38 +1,160 @@
 # [Abstract syntax](@id abstract)
 
-Todo:
-- give the full grammar
-- thanks MLStyle
-- examples for most features
-- link to example + API for functional syntax
+The full grammar of OptimalControl.jl DSL[^1] is given below. The idea is to use a syntax that is
+- pure Julia (and, as such, effortlessly analysed by the standard Julia parser),
+- as close as possible to the mathematical description of an optimal control problem. 
+
+While the syntax will be transparent to those users familiar with Julia expressions, we provide examples for every case that should be widely understandable. We rely heavily on [MLStyle.jl](https://github.com/thautwarm/MLStyle.jl) and its pattern matching abilities for the semantic pass. Abstract definitions use the macro `@def`.
 
-Let us consider a wagon moving along a rail, whom acceleration can be controlled by a force $u$.
-We denote by $x = (x_1, x_2)$ the state of the wagon, that is its position $x_1$ and its velocity $x_2$.
+## Variable declaration
 
-```@raw html
-<img src="./assets/chariot.png" style="display: block; margin: 0 auto 20px auto;" width="300px">
+```julia
+   # variable                    
+   :( $v ∈ R^$q, variable            ) 
+   :( $v ∈ R   , variable            ) 
 ```
 
-We assume that the mass is constant and unitary and that there is no friction. The dynamics we consider is given by
+A variable (only one is allowed) is a finite dimensional vector or reals that will be optimised along with state and control values. To define an (almost empty!) optimal control problem, named `ocp`, having a dimension two variable named `v`, do the following:
 
-```math
-    \dot x_1(t) = x_2(t), \quad \dot x_2(t) = u(t),\quad u(t) \in \R,
+```@example main
+using OptimalControl #hide
+@def ocp begin
+    v ∈ R^2, variable
+end
 ```
 
-which is simply the [double integrator](https://en.wikipedia.org/w/index.php?title=Double_integrator&oldid=1071399674) system.
-Les us consider a transfer starting at time $t_0 = 0$ and ending at time $t_f = 1$, for which we want to minimise the transfer energy
+Aliases `v₁` and `v₂` are automatically defined and can be used in subsequent expressions instead of `v[1]` and `v[2]`. The user can also define her own aliases for the components (one per dimension):
 
-```math
-    \frac{1}{2}\int_{0}^{1} u^2(t) \, \mathrm{d}t
+```@example main
+@def ocp begin
+    v = (a, b) ∈ R^2, variable
+end
 ```
 
-starting from the condition $x(0) = (-1, 0)$ and with the goal to reach the target $x(1) = (0, 0)$.
+## Time declaration
 
-First, we need to import the `OptimalControl.jl` package to define the optimal control problem and `NLPModelsIpopt.jl` to solve it. 
-We also need to import the `Plots.jl` package to plot the solution.
+```julia
+   :( $t ∈ [ $t0, $tf ], time        ) 
+```
+
+The independent variable or *time* is a scalar bound to a given interval. Its name is arbitrary.
 
 ```@example main
-using OptimalControl
-using NLPModelsIpopt
-using Plots
+tf = 5
+@def ocp begin
+    t ∈ [0, tf], time
+end
+```
+
+One (or even the two bounds) can be variable, typically for minimum time problems:
+
+```@example main
+@def ocp begin
+    v = (T, λ) ∈ R^2, variable
+    t ∈ [0, T], time
+end
+```
+
+## State declaration
+```julia
+   :( $x ∈ R^$n, state               ) 
+   :( $x ∈ R   , state               ) 
+```
+
+As for the variable, there are automatic aliases (`x₁` for `x[1]`, *etc.*) and the used can define her own aliases (one per scalar component of the state):
+
+XXXXX
+
+## Control declaration
+```julia
+   :( $u ∈ R^$m, control             ) 
+   :( $u ∈ R   , control             ) 
 ```
+
+As for the variable and the state, automatic or user-defined aliases are available for control components.
+
+## Dynamics
+```julia
+   :( ∂($x)($t) == $e1               ) 
+   :( ∂($x)($t) == $e1, $label       ) 
+```
+- add x\dot (automatic alias)
+- no line to line / component declaration
+- labeled or not
+
+## Constraints
+- boundary, control, state, mixed
+- ranges (linear) of general (*a priori* nonlinear)
+- equalities, one or two-sided inequalities
+- labeled or not (future use: retrieve associated multipliers on the discretised problem)
+
+```julia
+   :( $e1 == $e2                     ) 
+   :( $e1 == $e2, $label             ) 
+   :( $e1 ≤  $e2 ≤  $e3              ) 
+   :( $e1 ≤  $e2 ≤  $e3, $label      ) 
+   :(        $e2 ≤  $e3              ) 
+   :(        $e2 ≤  $e3, $label      ) 
+   :( $e3 ≥  $e2 ≥  $e1              ) 
+   :( $e3 ≥  $e2 ≥  $e1, $label      ) 
+   :( $e2 ≥  $e1                     ) 
+   :( $e2 ≥  $e1,        $label      ) 
+```
+
+# Mayer cost
+ ```julia                                      
+   :( $e1                      → min ) 
+   :( $e1                      → max ) 
+```
+
+# Lagrange cost
+```julia
+   :(             ∫($e1)       → min ) 
+   :(           - ∫($e1)       → min ) 
+   :(       $e1 * ∫($e2)       → min ) 
+                                       
+   :(             ∫($e1)       → max ) 
+   :(           - ∫($e1)       → max ) 
+   :(       $e1 * ∫($e2)       → max ) 
+```
+- caveat:  ... + ... + ...
+
+# Bolza cost
+```julia
+   :( $e1 +       ∫($e2)       → min ) 
+   :( $e1 + $e2 * ∫($e3)       → min ) 
+                                       
+   :( $e1 -       ∫($e2)       → min ) 
+   :( $e1 - $e2 * ∫($e3)       → min ) 
+                                       
+   :( $e1 +       ∫($e2)       → max ) 
+   :( $e1 + $e2 * ∫($e3)       → max ) 
+                                       
+   :( $e1 -       ∫($e2)       → max ) 
+   :( $e1 - $e2 * ∫($e3)       → max ) 
+                                       
+   :(             ∫($e2) + $e1 → min ) 
+   :(       $e2 * ∫($e3) + $e1 → min ) 
+                                       
+   :(             ∫($e2) - $e1 → min ) 
+   :(       $e2 * ∫($e3) - $e1 → min ) 
+                                       
+   :(             ∫($e2) + $e1 → max ) 
+   :(       $e2 * ∫($e3) + $e1 → max ) 
+                                       
+   :(             ∫($e2) - $e1 → max ) 
+   :(       $e2 * ∫($e3) - $e1 → max ) 
+ ```
+
+## Aliases
+
+- order: declaration first, then constraint and cost (no ordering for these two)
+- examples for most features
+- caveat's (check isssue) (case by base)
+- error should be OK (give an example)
+- expressions should evaluate at run
+- aliases (vars, and in general)
+- link to example + API for functional syntax
+- point towards examples for further use
+
+[^1]: Domain Specific Language
\ No newline at end of file
diff --git a/docs/src/tutorial-goddard.md b/docs/src/tutorial-goddard.md
index 9892af35..61b05d19 100644
--- a/docs/src/tutorial-goddard.md
+++ b/docs/src/tutorial-goddard.md
@@ -103,7 +103,7 @@ nothing # hide
 We then solve it
 
 ```@example main
-direct_sol = solve(ocp; grid_size=100)
+direct_sol = solve(ocp; grid_size=100, linear_solver="mumps")
 nothing # hide
 ```
 
diff --git a/docs/src/tutorial-nlp.md b/docs/src/tutorial-nlp.md.foo
similarity index 100%
rename from docs/src/tutorial-nlp.md
rename to docs/src/tutorial-nlp.md.foo

From 0842f8eee75c72a252b41325f16e3fb422cd6d6b Mon Sep 17 00:00:00 2001
From: Jean-Baptiste Caillau <jean-baptiste.caillau@univ-cotedazur.fr>
Date: Sat, 3 Aug 2024 11:54:10 +0200
Subject: [PATCH 05/17] docs

---
 docs/src/{tutorial-nlp.md.foo => tutorial-nlp.md} | 0
 1 file changed, 0 insertions(+), 0 deletions(-)
 rename docs/src/{tutorial-nlp.md.foo => tutorial-nlp.md} (100%)

diff --git a/docs/src/tutorial-nlp.md.foo b/docs/src/tutorial-nlp.md
similarity index 100%
rename from docs/src/tutorial-nlp.md.foo
rename to docs/src/tutorial-nlp.md

From 02439332a2a37f96e50ab53b727309e2e17b8603 Mon Sep 17 00:00:00 2001
From: Jean-Baptiste Caillau <jean-baptiste.caillau@univ-cotedazur.fr>
Date: Sat, 3 Aug 2024 16:24:28 +0200
Subject: [PATCH 06/17] docs

---
 docs/Project.toml             |  1 +
 docs/make.jl                  |  2 +-
 docs/src/tutorial-abstract.md | 40 +++++++++++++++++++++++++++++------
 3 files changed, 35 insertions(+), 8 deletions(-)

diff --git a/docs/Project.toml b/docs/Project.toml
index 94aa6a73..f4abba17 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -13,6 +13,7 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MINPACK = "4854310b-de5a-5eb6-a2a5-c1dee2bd17f9"
 NLPModelsIpopt = "f4238b75-b362-5c4c-b852-0801c9a21d71"
 NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
+OptimalControl = "5f98b655-cc9a-415a-b60e-744165666948"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 Percival = "01435c0c-c90d-11e9-3788-63660f8fbccc"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
diff --git a/docs/make.jl b/docs/make.jl
index 1a4bbb15..d0188e24 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -35,7 +35,7 @@ makedocs(
             "tutorial-double-integrator.md",
             "tutorial-initial-guess.md",
             "tutorial-continuation.md",
-            #"tutorial-nlp.md",
+            "tutorial-nlp.md",
             "tutorial-plot.md",
             "tutorial-lqr-basic.md",
             "tutorial-iss.md",
diff --git a/docs/src/tutorial-abstract.md b/docs/src/tutorial-abstract.md
index 9a1bfbe8..a5d5b9fd 100644
--- a/docs/src/tutorial-abstract.md
+++ b/docs/src/tutorial-abstract.md
@@ -19,7 +19,7 @@ A variable (only one is allowed) is a finite dimensional vector or reals that wi
 ```@example main
 using OptimalControl #hide
 @def ocp begin
-    v ∈ R^2, variable
+    v ∈ R², variable
 end
 ```
 
@@ -27,14 +27,14 @@ Aliases `v₁` and `v₂` are automatically defined and can be used in subsequen
 
 ```@example main
 @def ocp begin
-    v = (a, b) ∈ R^2, variable
+    v = (a, b) ∈ R², variable
 end
 ```
 
 ## Time declaration
 
 ```julia
-   :( $t ∈ [ $t0, $tf ], time        ) 
+   :( $t ∈ [$t0, $tf], time        ) 
 ```
 
 The independent variable or *time* is a scalar bound to a given interval. Its name is arbitrary.
@@ -50,7 +50,7 @@ One (or even the two bounds) can be variable, typically for minimum time problem
 
 ```@example main
 @def ocp begin
-    v = (T, λ) ∈ R^2, variable
+    v = (T, λ) ∈ R², variable
     t ∈ [0, T], time
 end
 ```
@@ -61,9 +61,21 @@ end
    :( $x ∈ R   , state               ) 
 ```
 
-As for the variable, there are automatic aliases (`x₁` for `x[1]`, *etc.*) and the used can define her own aliases (one per scalar component of the state):
+The state declaration defines the name and the dimension of the state:
 
-XXXXX
+```@example main
+@def ocp begin
+       x ∈ R⁴, state
+end
+```
+
+As for the variable, there are automatic aliases (`x₁` for `x[1]`, *etc.*) and the user can define her own aliases (one per scalar component of the state):
+
+```@example main
+@def ocp begin
+       x = (q₁, q₂, v₁, v₂) ∈ R⁴, state
+end
+```
 
 ## Control declaration
 ```julia
@@ -71,7 +83,21 @@ XXXXX
    :( $u ∈ R   , control             ) 
 ```
 
-As for the variable and the state, automatic or user-defined aliases are available for control components.
+The control declaration defines the name and the dimension of the control:
+
+```@example main
+@def ocp begin
+    u ∈ R², control
+end
+```
+
+As before, there are automatic aliases (`u₁` for `u[1]`, *etc.*) and the user can define her own aliases (one per scalar component of the state):
+
+```@example main
+@def ocp begin
+    u = (α, β) ∈ R², control
+end
+```
 
 ## Dynamics
 ```julia

From 54281a8889018c36d5874260a03fcab9e50e8f75 Mon Sep 17 00:00:00 2001
From: Jean-Baptiste Caillau <jean-baptiste.caillau@univ-cotedazur.fr>
Date: Sat, 3 Aug 2024 17:41:07 +0200
Subject: [PATCH 07/17] doc

---
 docs/Project.toml                  |   1 -
 docs/src/tutorial-abstract.md      | 121 ++++++++++++++++++++++-------
 docs/src/tutorial-basic-example.md |   2 +-
 docs/src/tutorial-continuation.md  |   2 +-
 docs/src/tutorial-goddard.md       |   2 +-
 docs/src/tutorial-initial-guess.md |   2 +-
 docs/src/tutorial-iss.md           |   2 +-
 docs/src/tutorial-lqr-basic.md     |   2 +-
 docs/src/tutorial-nlp.md           |   4 +-
 docs/src/tutorial-plot.md          |  20 ++---
 10 files changed, 109 insertions(+), 49 deletions(-)

diff --git a/docs/Project.toml b/docs/Project.toml
index f4abba17..94aa6a73 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -13,7 +13,6 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MINPACK = "4854310b-de5a-5eb6-a2a5-c1dee2bd17f9"
 NLPModelsIpopt = "f4238b75-b362-5c4c-b852-0801c9a21d71"
 NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
-OptimalControl = "5f98b655-cc9a-415a-b60e-744165666948"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 Percival = "01435c0c-c90d-11e9-3788-63660f8fbccc"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
diff --git a/docs/src/tutorial-abstract.md b/docs/src/tutorial-abstract.md
index a5d5b9fd..a6d04675 100644
--- a/docs/src/tutorial-abstract.md
+++ b/docs/src/tutorial-abstract.md
@@ -6,7 +6,7 @@ The full grammar of OptimalControl.jl DSL[^1] is given below. The idea is to use
 
 While the syntax will be transparent to those users familiar with Julia expressions, we provide examples for every case that should be widely understandable. We rely heavily on [MLStyle.jl](https://github.com/thautwarm/MLStyle.jl) and its pattern matching abilities for the semantic pass. Abstract definitions use the macro `@def`.
 
-## Variable declaration
+## Variable
 
 ```julia
    # variable                    
@@ -31,7 +31,15 @@ Aliases `v₁` and `v₂` are automatically defined and can be used in subsequen
 end
 ```
 
-## Time declaration
+A one dimensional can be declared according to
+
+```@example main
+@def ocp begin
+    v ∈ R, variable
+end
+```
+
+## Time
 
 ```julia
    :( $t ∈ [$t0, $tf], time        ) 
@@ -55,7 +63,8 @@ One (or even the two bounds) can be variable, typically for minimum time problem
 end
 ```
 
-## State declaration
+## State
+
 ```julia
    :( $x ∈ R^$n, state               ) 
    :( $x ∈ R   , state               ) 
@@ -77,7 +86,8 @@ As for the variable, there are automatic aliases (`x₁` for `x[1]`, *etc.*) and
 end
 ```
 
-## Control declaration
+## Control
+
 ```julia
    :( $u ∈ R^$m, control             ) 
    :( $u ∈ R   , control             ) 
@@ -100,80 +110,131 @@ end
 ```
 
 ## Dynamics
+
 ```julia
    :( ∂($x)($t) == $e1               ) 
-   :( ∂($x)($t) == $e1, $label       ) 
 ```
-- add x\dot (automatic alias)
-- no line to line / component declaration
-- labeled or not
+
+The dynamics is given in the standard vectorial ODE form:
+
+```math
+    \dot{x}(t) = f(x(t), u(t)) \quad \text{or} \quad f(t, x(t), u(t))
+```
+
+depending on whether it is autonomous or not (the parser will detect dependence time, which entails that time and state must be declared prior to dynamics - an error will be issued otherwise). The symbol `∂` or the dotted state name
+(`ẋ` can be used):
+
+```@example main
+@def ocp begin
+    t ∈ [0, 1], time
+    x ∈ R², state
+    u ∈ R, control
+    ∂(x)(t) == [x₂(t), u(t)]
+end
+```
+
+or
+
+```@example main
+@def ocp begin
+    t ∈ [0, 1], time
+    x ∈ R², state
+    u ∈ R, control
+    ẋ(t) == [x₂(t), u(t)]
+end
+```
+
+Any Julia code can be used, so the following is also OK: 
+
+```@example main
+\cdot  
+@def ocp begin
+    t ∈ [0, 1], time
+    x ∈ R², state
+    u ∈ R, control
+    ẋ(t) == F₀(x(t)) + u * F₁(x(t))
+end
+
+F₀(x) = [x[2], 0]
+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 [section ](@aliases) below):
+
+```@example main
+def double_integrator begin
+    tf ∈ R, variable
+    t ∈ [0, tf], time
+    x = (q, v) ∈ R², state
+    u ∈ R, control
+    q̇ = v
+    v̇ = u
+    ẋ(t) == [q̇, v̇]
+end
+```
 
 ## Constraints
-- boundary, control, state, mixed
-- ranges (linear) of general (*a priori* nonlinear)
-- equalities, one or two-sided inequalities
-- labeled or not (future use: retrieve associated multipliers on the discretised problem)
 
 ```julia
    :( $e1 == $e2                     ) 
-   :( $e1 == $e2, $label             ) 
    :( $e1 ≤  $e2 ≤  $e3              ) 
-   :( $e1 ≤  $e2 ≤  $e3, $label      ) 
    :(        $e2 ≤  $e3              ) 
-   :(        $e2 ≤  $e3, $label      ) 
    :( $e3 ≥  $e2 ≥  $e1              ) 
-   :( $e3 ≥  $e2 ≥  $e1, $label      ) 
    :( $e2 ≥  $e1                     ) 
-   :( $e2 ≥  $e1,        $label      ) 
 ```
 
-# Mayer cost
- ```julia                                      
+Constraints
+- boundary, control, state, mixed
+- ranges (linear) of general (*a priori* nonlinear)
+- equalities, one or two-sided inequalities
+- labeled or not (future use: retrieve associated multipliers on the discretised problem)
+
+## Mayer cost
+
+```julia                                      
    :( $e1                      → min ) 
    :( $e1                      → max ) 
 ```
 
-# Lagrange cost
+## Lagrange cost
+
 ```julia
    :(             ∫($e1)       → min ) 
    :(           - ∫($e1)       → min ) 
    :(       $e1 * ∫($e2)       → min ) 
-                                       
    :(             ∫($e1)       → max ) 
    :(           - ∫($e1)       → max ) 
    :(       $e1 * ∫($e2)       → max ) 
 ```
 - caveat:  ... + ... + ...
 
-# Bolza cost
+## Bolza cost
+
 ```julia
    :( $e1 +       ∫($e2)       → min ) 
    :( $e1 + $e2 * ∫($e3)       → min ) 
-                                       
    :( $e1 -       ∫($e2)       → min ) 
    :( $e1 - $e2 * ∫($e3)       → min ) 
-                                       
    :( $e1 +       ∫($e2)       → max ) 
    :( $e1 + $e2 * ∫($e3)       → max ) 
-                                       
    :( $e1 -       ∫($e2)       → max ) 
    :( $e1 - $e2 * ∫($e3)       → max ) 
-                                       
    :(             ∫($e2) + $e1 → min ) 
    :(       $e2 * ∫($e3) + $e1 → min ) 
-                                       
    :(             ∫($e2) - $e1 → min ) 
    :(       $e2 * ∫($e3) - $e1 → min ) 
-                                       
    :(             ∫($e2) + $e1 → max ) 
    :(       $e2 * ∫($e3) + $e1 → max ) 
-                                       
    :(             ∫($e2) - $e1 → max ) 
    :(       $e2 * ∫($e3) - $e1 → max ) 
  ```
 
-## Aliases
+## [Aliases](@id aliases)
 
+- labels for dyn and constraints (future use, retrieve multipliers...)
 - order: declaration first, then constraint and cost (no ordering for these two)
 - examples for most features
 - caveat's (check isssue) (case by base)
diff --git a/docs/src/tutorial-basic-example.md b/docs/src/tutorial-basic-example.md
index cf7f8d0d..ea940443 100644
--- a/docs/src/tutorial-basic-example.md
+++ b/docs/src/tutorial-basic-example.md
@@ -35,7 +35,7 @@ Then, we can define the problem
 
 ```@example main
 @def ocp begin
-    t ∈ [ 0, 1 ], time
+    t ∈ [0, 1], time
     x ∈ R², state
     u ∈ R, control
     x(0) == [ -1, 0 ]
diff --git a/docs/src/tutorial-continuation.md b/docs/src/tutorial-continuation.md
index a3b5efd6..230a0981 100644
--- a/docs/src/tutorial-continuation.md
+++ b/docs/src/tutorial-continuation.md
@@ -24,7 +24,7 @@ and write a function that returns the OCP for a given final time
 ```@example main
 function ocp_T(T)
     @def ocp begin
-        t ∈ [ 0, T ], time
+        t ∈ [0, T], time
         x ∈ R², state
         u ∈ R, control
         q = x₁
diff --git a/docs/src/tutorial-goddard.md b/docs/src/tutorial-goddard.md
index 61b05d19..6b65c165 100644
--- a/docs/src/tutorial-goddard.md
+++ b/docs/src/tutorial-goddard.md
@@ -65,7 +65,7 @@ mf = 0.6    # final mass to target
 @def ocp begin # definition of the optimal control problem
 
     tf ∈ R, variable
-    t ∈ [ t0, tf ], time
+    t ∈ [t0, tf], time
     x = (r, v, m) ∈ R³, state
     u ∈ R, control
 
diff --git a/docs/src/tutorial-initial-guess.md b/docs/src/tutorial-initial-guess.md
index 4344ac01..27c347a2 100644
--- a/docs/src/tutorial-initial-guess.md
+++ b/docs/src/tutorial-initial-guess.md
@@ -24,7 +24,7 @@ tf = 10
 α  = 5
 
 @def ocp begin
-    t ∈ [ t0, tf ], time
+    t ∈ [t0, tf], time
     v ∈ R, variable
     x ∈ R², state
     u ∈ R, control
diff --git a/docs/src/tutorial-iss.md b/docs/src/tutorial-iss.md
index 2111718b..1a3d0085 100644
--- a/docs/src/tutorial-iss.md
+++ b/docs/src/tutorial-iss.md
@@ -38,7 +38,7 @@ xf = 0
 α  = 1.5
 @def ocp begin
 
-    t ∈ [ t0, tf ], time
+    t ∈ [t0, tf], time
     x ∈ R, state
     u ∈ R, control
 
diff --git a/docs/src/tutorial-lqr-basic.md b/docs/src/tutorial-lqr-basic.md
index 10258d57..050c07b0 100644
--- a/docs/src/tutorial-lqr-basic.md
+++ b/docs/src/tutorial-lqr-basic.md
@@ -51,7 +51,7 @@ B  = [ 0
 function lqr(tf)
 
     @def ocp begin
-        t ∈ [ 0, tf ], time
+        t ∈ [0, tf], time
         x ∈ R², state
         u ∈ R, control
         x(0) == x0
diff --git a/docs/src/tutorial-nlp.md b/docs/src/tutorial-nlp.md
index 612f8df1..21e912d7 100644
--- a/docs/src/tutorial-nlp.md
+++ b/docs/src/tutorial-nlp.md
@@ -28,7 +28,7 @@ We define a test problem
 ```@example main
 @def ocp begin
 
-    t ∈ [ 0, 1 ], time
+    t ∈ [0, 1], time
     x ∈ R², state
     u ∈ R, control
 
@@ -97,4 +97,4 @@ using Percival
 
 output = percival(nlp)
 print(output)
-```
\ No newline at end of file
+```
diff --git a/docs/src/tutorial-plot.md b/docs/src/tutorial-plot.md
index 90e0a67c..a6caf6dc 100644
--- a/docs/src/tutorial-plot.md
+++ b/docs/src/tutorial-plot.md
@@ -14,14 +14,14 @@ Then, we define a simple optimal control problem and solve it.
 ```@example main
 @def ocp begin
 
-    t ∈ [ 0, 1 ], time
+    t ∈ [0, 1], time
     x ∈ R², state
     u ∈ R, control
 
-    x(0) == [ -1, 0 ]
-    x(1) == [ 0, 0 ]
+    x(0) == [-1, 0]
+    x(1) == [0, 0]
 
-    ẋ(t) == [ x₂(t), u(t) ]
+    ẋ(t) == [x₂(t), u(t)]
 
     ∫( 0.5u(t)^2 ) → min
 
@@ -112,14 +112,14 @@ You can plot the solution of a second optimal control problem on the same figure
 ```@example main
 @def ocp begin
 
-    t ∈ [ 0, 1 ], time
+    t ∈ [0, 1], time
     x ∈ R², state
     u ∈ R, control
 
-    x(0) == [ -0.5, -0.5 ]
-    x(1) == [ 0, 0 ]
+    x(0) == [-0.5, -0.5]
+    x(1) == [0, 0]
 
-    ẋ(t) == [ x₂(t), u(t) ]
+    ẋ(t) == [x₂(t), u(t)]
 
     ∫( 0.5u(t)^2 ) → min
 
@@ -187,7 +187,7 @@ B  = [ 0
 function lqr(tf)
 
     @def ocp begin
-        t ∈ [ 0, tf ], time
+        t ∈ [0, tf], time
         x ∈ R², state
         u ∈ R, control
         x(0) == x0
@@ -220,4 +220,4 @@ pu  = plot(plt[5]; legend=false, xlabel="s", ylabel="u")
 
 using Plots.PlotMeasures # for leftmargin, bottommargin
 plot(px1, px2, pu; layout=(1, 3), size=(800, 300), leftmargin=5mm, bottommargin=5mm)
-```
\ No newline at end of file
+```

From 2cba25f4b2c81622bc9cca08f50f5546ed99c55d Mon Sep 17 00:00:00 2001
From: Jean-Baptiste Caillau <jean-baptiste.caillau@univ-cotedazur.fr>
Date: Sat, 3 Aug 2024 18:56:06 +0200
Subject: [PATCH 08/17] doc

---
 docs/Project.toml                      |   1 +
 docs/src/tutorial-abstract.md          | 120 +++++++++++++++++++++++--
 docs/src/tutorial-double-integrator.md |   2 +-
 3 files changed, 113 insertions(+), 10 deletions(-)

diff --git a/docs/Project.toml b/docs/Project.toml
index 94aa6a73..f4abba17 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -13,6 +13,7 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MINPACK = "4854310b-de5a-5eb6-a2a5-c1dee2bd17f9"
 NLPModelsIpopt = "f4238b75-b362-5c4c-b852-0801c9a21d71"
 NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
+OptimalControl = "5f98b655-cc9a-415a-b60e-744165666948"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 Percival = "01435c0c-c90d-11e9-3788-63660f8fbccc"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
diff --git a/docs/src/tutorial-abstract.md b/docs/src/tutorial-abstract.md
index a6d04675..af0f1b0a 100644
--- a/docs/src/tutorial-abstract.md
+++ b/docs/src/tutorial-abstract.md
@@ -147,7 +147,6 @@ end
 Any Julia code can be used, so the following is also OK: 
 
 ```@example main
-\cdot  
 @def ocp begin
     t ∈ [0, 1], time
     x ∈ R², state
@@ -165,13 +164,13 @@ F₁(x) = [0, 1]
 Currently, it is not possible to declare the dynamics component after component, but a simple workaround is to use *aliases* (check the relevant [section ](@aliases) below):
 
 ```@example main
-def double_integrator begin
+@def damped_integrator begin
     tf ∈ R, variable
     t ∈ [0, tf], time
     x = (q, v) ∈ R², state
     u ∈ R, control
-    q̇ = v
-    v̇ = u
+    q̇ = v(t)
+    v̇ = u(t) - c(t)
     ẋ(t) == [q̇, v̇]
 end
 ```
@@ -186,11 +185,53 @@ end
    :( $e2 ≥  $e1                     ) 
 ```
 
-Constraints
-- boundary, control, state, mixed
-- ranges (linear) of general (*a priori* nonlinear)
-- equalities, one or two-sided inequalities
-- labeled or not (future use: retrieve associated multipliers on the discretised problem)
+Constraints 
+- can be of five types: boundary, control, state, mixed, variable,
+- be linear (ranges) or nonlinear (not ranges),
+- equalities or (one or two-sided) inequalities.
+
+Boundary conditions are detected when the expression contains evaluations of the state at initial and / or final time bounds (*e.g.*, `x(0)`), and may not involve the control. Conversely control, state or mixed constraints will involve control, state or both evaluated at the declared time (*e.g.*, `x(t) + u(t)`). 
+Other combinations should be detected as incorrect by the parser 🤞🏾. The variable may be involved in any of the four previous constraints. Constraints involving the variable only are variable constraints, either linear or nonlinear.
+In the example below, there are
+- two linear boundary constraints,
+- one linear variable constraint,
+- one linear state constraint,
+- one (two-sided) nonlinear control constraint.
+
+```@example main
+@def ocp begin
+    tf ∈ R, variable
+    t ∈ [0, tf], time
+    x ∈ R², state
+    u ∈ R, control
+    x(0) == [-1, 0]
+    x(tf) == [0, 0]
+    ẋ(t) == [x₂(t), u(t)]
+    tf ≥ 0 
+    x₂(t) ≤ 1
+    u(t)^2 ≤ 1
+   end
+```
+
+!!! caveat
+    Write either `u(t)^2` or `(u^2)(t)`, not `u^2(t)` since in Julia the latter is means `u^(2t)`. Moreover,
+    in the case of equalities or of one-sided inequalities, the control and / or the state must belong the *left-hand side*. The following will error:
+
+```julia
+julia> @def ocp begin
+    t ∈ [0, 2], time
+    x ∈ R², state
+    u ∈ R, control
+    x(0) == [-1, 0]
+    x(2) == [0, 0]
+    ẋ(t) == [x₂(t), u(t)]
+    1 ≤ x₂(t)
+    -1 ≤ u(t) ≤ 1
+end
+ERROR: ParsingError: 
+Line 7: 1 ≤ x₂(t)
+bad constraint declaration
+```
 
 ## Mayer cost
 
@@ -199,6 +240,28 @@ Constraints
    :( $e1                      → max ) 
 ```
 
+Mayer costs are defined in a similar way to boundary conditions and follow the same rules. The symbol `→` is used
+to denote minimisation of maximisation, the latter being treated by minimising the opposite cost.
+
+```@example main
+@def ocp begin
+    tf ∈ R, variable
+    t ∈ [0, tf], time
+    x = (q, v) ∈ R², state
+    u ∈ R, control
+    tf ≥ 0
+    -1 ≤ u(t) ≤ 1
+    q(0) == 1
+    v(0) == 2
+    q(tf) == 0
+    v(tf) == 0
+    0 ≤ q(t) ≤ 5
+   -2 ≤ v(t) ≤ 3
+    ẋ(t) == [v(t), u(t)]
+    tf → min
+end
+```
+
 ## Lagrange cost
 
 ```julia
@@ -235,6 +298,43 @@ Constraints
 ## [Aliases](@id aliases)
 
 - labels for dyn and constraints (future use, retrieve multipliers...)
+
+!!! hint
+    You can use a trace mode for the macro `@def` to look at your code after expansions of the aliases adding `true` after your `begin ... end` block:
+
+```@example main
+@def damped_integrator begin
+    tf ∈ R, variable
+    t ∈ [0, tf], time
+    x = (q, v) ∈ R², state
+    u ∈ R, control
+    q̇ = v(t)
+    v̇ = u(t) - c(t)
+    ẋ(t) == [q̇, v̇]
+end true
+```
+
+!!! caveat
+    The dynamics of an OCP is indeed a particular constraint, be careful to use `==` and not a single `=` that would try define an alias:
+
+```julia
+julia> @def double_integrator begin
+       tf ∈ R, variable
+       t ∈ [0, tf], time
+       x = (q, v) ∈ R², state
+       u ∈ R, control
+       q̇ = v
+       v̇ = u
+       ẋ(t) = [q̇, v̇]
+       end
+ERROR: ParsingError: 
+Line 7: ẋ(t) = begin
+        #= REPL[35]:8 =#
+        [q̇, v̇]
+    end
+forbidden alias name: (∂(x))(t)
+```
+
 - order: declaration first, then constraint and cost (no ordering for these two)
 - examples for most features
 - caveat's (check isssue) (case by base)
@@ -242,6 +342,8 @@ Constraints
 - expressions should evaluate at run
 - aliases (vars, and in general)
 - link to example + API for functional syntax
+- parsing error should be explicit
 - point towards examples for further use
 
+
 [^1]: Domain Specific Language
\ No newline at end of file
diff --git a/docs/src/tutorial-double-integrator.md b/docs/src/tutorial-double-integrator.md
index e7d74995..8e49755c 100644
--- a/docs/src/tutorial-double-integrator.md
+++ b/docs/src/tutorial-double-integrator.md
@@ -68,7 +68,7 @@ nothing # hide
 Solve it
 
 ```@example main
-sol = solve(ocp; grid_size=200)
+sol = solve(ocp; grid_size=100)
 nothing # hide
 ```
 

From ba90ccfef642ed0c25248393cfb7508a4f6c08e6 Mon Sep 17 00:00:00 2001
From: Jean-Baptiste Caillau <jean-baptiste.caillau@univ-cotedazur.fr>
Date: Sun, 4 Aug 2024 16:14:29 +0200
Subject: [PATCH 09/17] doc

---
 docs/src/tutorial-abstract.md | 127 ++++++++++++++++++++++++++++------
 1 file changed, 105 insertions(+), 22 deletions(-)

diff --git a/docs/src/tutorial-abstract.md b/docs/src/tutorial-abstract.md
index af0f1b0a..8b45adfb 100644
--- a/docs/src/tutorial-abstract.md
+++ b/docs/src/tutorial-abstract.md
@@ -1,15 +1,14 @@
 # [Abstract syntax](@id abstract)
 
-The full grammar of OptimalControl.jl DSL[^1] is given below. The idea is to use a syntax that is
+The full grammar of OptimalControl.jl small Domain Specific Language is given below. The idea is to use a syntax that is
 - pure Julia (and, as such, effortlessly analysed by the standard Julia parser),
 - as close as possible to the mathematical description of an optimal control problem. 
 
-While the syntax will be transparent to those users familiar with Julia expressions, we provide examples for every case that should be widely understandable. We rely heavily on [MLStyle.jl](https://github.com/thautwarm/MLStyle.jl) and its pattern matching abilities for the semantic pass. Abstract definitions use the macro `@def`.
+While the syntax will be transparent to those users familiar with Julia expressions (`Expr`'s), we provide examples for every case that should be widely understandable. We rely heavily on [MLStyle.jl](https://github.com/thautwarm/MLStyle.jl) and its pattern matching abilities for the semantic pass. Abstract definitions use the macro `@def`.
 
-## Variable
+## [Variable](@id variable)
 
 ```julia
-   # variable                    
    :( $v ∈ R^$q, variable            ) 
    :( $v ∈ R   , variable            ) 
 ```
@@ -63,7 +62,7 @@ One (or even the two bounds) can be variable, typically for minimum time problem
 end
 ```
 
-## State
+## [State](@id state)
 
 ```julia
    :( $x ∈ R^$n, state               ) 
@@ -86,7 +85,7 @@ As for the variable, there are automatic aliases (`x₁` for `x[1]`, *etc.*) and
 end
 ```
 
-## Control
+## [Control](@id control)
 
 ```julia
    :( $u ∈ R^$m, control             ) 
@@ -109,7 +108,7 @@ As before, there are automatic aliases (`u₁` for `u[1]`, *etc.*) and the user
 end
 ```
 
-## Dynamics
+## [Dynamics](@id dynamics)
 
 ```julia
    :( ∂($x)($t) == $e1               ) 
@@ -151,7 +150,7 @@ Any Julia code can be used, so the following is also OK:
     t ∈ [0, 1], time
     x ∈ R², state
     u ∈ R, control
-    ẋ(t) == F₀(x(t)) + u * F₁(x(t))
+    ẋ(t) == F₀(x(t)) + u(t) * F₁(x(t))
 end
 
 F₀(x) = [x[2], 0]
@@ -161,7 +160,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 [section ](@aliases) below):
+Currently, it is not possible to declare the dynamics component after component, but a simple workaround is to use *aliases* (check the relevant [section ](#aliases) below):
 
 ```@example main
 @def damped_integrator begin
@@ -272,7 +271,24 @@ end
    :(           - ∫($e1)       → max ) 
    :(       $e1 * ∫($e2)       → max ) 
 ```
-- caveat:  ... + ... + ...
+
+Lagrange (integral) costs are defined used the symbol `∫`, *with parenthesis:
+
+```@example main
+@def ocp begin
+    t ∈ [0, 1], time
+    x = (q, v) ∈ R², state
+    u ∈ R, control
+    0.5∫( q(t) + u(t)^2 ) → min
+end
+```
+
+The integration range is implicitly equal to the time range, so the cost above is to be understood as
+```math
+\int_0^1 q(t) + u^2(t)\,\mathrm{d}t \to \min.
+```
+
+As for the dynamics, the parser will detect whether the integrand depends or not on time (autonomous / non-autonomous case).
 
 ## Bolza cost
 
@@ -295,9 +311,66 @@ end
    :(       $e2 * ∫($e3) - $e1 → max ) 
  ```
 
+Quite readily, Mayer and Lagrange costs can be combined into genral Bolza costs. For instance as follows:
+
+```@example main
+@def ocp begin
+    p = (t0, tf) ∈ R², variable
+    t ∈ [t0, tf], time
+    x = (q, v) ∈ R², state
+    u ∈ R², control
+    (tf - t0) + 0.5∫( c(t) * u(t)^2 ) → min
+end
+```
+
+!!! caveat
+    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
+julia> @def ocp begin
+           p = (t0, tf) ∈ R², variable
+           t ∈ [t0, tf], time
+           x = (q, v) ∈ R², state
+           u ∈ R², control
+           (tf - t0) + q(tf) + 0.5∫( c(t) * u(t)^2 ) → min
+       end
+ERROR: ParsingError: 
+Line 5: (tf - t0) + q(tf) + 0.5 * ∫(c(t) * u(t) ^ 2) → min
+bad objective declaration resulting in a Mayer term with trailing ∫
+```
+
+The correct syntax is
+```@example main
+@def ocp begin
+    p = (t0, tf) ∈ R², variable
+    t ∈ [t0, tf], time
+    x = (q, v) ∈ R², state
+    u ∈ R², control
+    ( (tf - t0) + q(tf) ) + 0.5∫( c(t) * u(t)^2 ) → min
+end
+```
+
 ## [Aliases](@id aliases)
 
-- labels for dyn and constraints (future use, retrieve multipliers...)
+```julia
+        :( $a = $e1 )
+```
+
+The single `=` symbol is used to define not a constraint but an alias, that is a purely syntactic replacement. There are some automatic aliases, *e.g.* `x₁` for `x[1]` if `x` is the state, and we have also seen that the user can define her own aliases when declaring the [variable](#variable), [state](#state) and [control](#control). Arbitrary aliases can be further defined, as below (compare with previous examples in the [dynamics](#dynamics) section):
+
+```@example main
+@def ocp begin
+    t ∈ [0, 1], time
+    x ∈ R², state
+    u ∈ R, control
+    F₀ = [x₂(t), 0]
+    F₁ = [0, 1]
+    ẋ(t) == F₀ + u(t) * F₁
+end
+
+```
+
+!!! caveat 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
     You can use a trace mode for the macro `@def` to look at your code after expansions of the aliases adding `true` after your `begin ... end` block:
@@ -311,7 +384,7 @@ end
     q̇ = v(t)
     v̇ = u(t) - c(t)
     ẋ(t) == [q̇, v̇]
-end true
+end
 ```
 
 !!! caveat
@@ -335,15 +408,25 @@ Line 7: ẋ(t) = begin
 forbidden alias name: (∂(x))(t)
 ```
 
-- order: declaration first, then constraint and cost (no ordering for these two)
-- examples for most features
-- caveat's (check isssue) (case by base)
-- error should be OK (give an example)
-- expressions should evaluate at run
-- aliases (vars, and in general)
-- link to example + API for functional syntax
-- parsing error should be explicit
-- point towards examples for further use
+## Misc
+
+- Declarations (variable - in any -, time, state and control) must be done first. Then, dynamics, constraints and cost can be introduced in an arbitrary order.
+- It is possible to provide numbers / labels (as math equations) for the constraints to improve readability (this is mostly for future use, typically to retrieve the Lagrange multiplier associated with the discretisation of a given constraint):
 
+```@example main
+@def damped_integrator begin
+    tf ∈ R, variable
+    t ∈ [0, tf], time
+    x = (q, v) ∈ R², state
+    u ∈ R, control
+    tf ≥ 0, (1)
+    q(0) == 2, (♡)
+    q̇ = v(t)
+    v̇ = u(t) - c(t)
+    ẋ(t) == [q̇, v̇]
+    x(t).^2  ≤ [1, 2], (state_con) 
+end
+```
 
-[^1]: Domain Specific Language
\ No newline at end of file
+- Parsing errors should be explicit enough (with line number in the `@def` `begin ... end` block indicated) 🤞🏾
+- Check tutorials and applications in the documentation for further use.
\ No newline at end of file

From 13d94e3b6971c8a438b71cead299205c9c044866 Mon Sep 17 00:00:00 2001
From: Jean-Baptiste Caillau <jean-baptiste.caillau@univ-cotedazur.fr>
Date: Sun, 4 Aug 2024 16:46:45 +0200
Subject: [PATCH 10/17] doc

---
 docs/Project.toml             |  1 -
 docs/src/tutorial-abstract.md | 72 ++++++++++++++++++-----------------
 2 files changed, 37 insertions(+), 36 deletions(-)

diff --git a/docs/Project.toml b/docs/Project.toml
index f4abba17..94aa6a73 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -13,7 +13,6 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MINPACK = "4854310b-de5a-5eb6-a2a5-c1dee2bd17f9"
 NLPModelsIpopt = "f4238b75-b362-5c4c-b852-0801c9a21d71"
 NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
-OptimalControl = "5f98b655-cc9a-415a-b60e-744165666948"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 Percival = "01435c0c-c90d-11e9-3788-63660f8fbccc"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
diff --git a/docs/src/tutorial-abstract.md b/docs/src/tutorial-abstract.md
index 8b45adfb..75e82bc3 100644
--- a/docs/src/tutorial-abstract.md
+++ b/docs/src/tutorial-abstract.md
@@ -1,19 +1,19 @@
 # [Abstract syntax](@id abstract)
 
-The full grammar of OptimalControl.jl small Domain Specific Language is given below. The idea is to use a syntax that is
+The full grammar of OptimalControl.jl small *Domain Specific Language* is given below. The idea is to use a syntax that is
 - pure Julia (and, as such, effortlessly analysed by the standard Julia parser),
 - as close as possible to the mathematical description of an optimal control problem. 
 
-While the syntax will be transparent to those users familiar with Julia expressions (`Expr`'s), we provide examples for every case that should be widely understandable. We rely heavily on [MLStyle.jl](https://github.com/thautwarm/MLStyle.jl) and its pattern matching abilities for the semantic pass. Abstract definitions use the macro `@def`.
+While the syntax will be transparent to those users familiar with Julia expressions (`Expr`'s), we provide examples for every case that should be widely understandable. We rely heavily on [MLStyle.jl](https://github.com/thautwarm/MLStyle.jl) and its pattern matching abilities 👍🏽 for the semantic pass. Abstract definitions use the macro `@def`.
 
 ## [Variable](@id variable)
 
 ```julia
-   :( $v ∈ R^$q, variable            ) 
-   :( $v ∈ R   , variable            ) 
+   :( $v ∈ R^$q, variable ) 
+   :( $v ∈ R   , variable ) 
 ```
 
-A variable (only one is allowed) is a finite dimensional vector or reals that will be optimised along with state and control values. To define an (almost empty!) optimal control problem, named `ocp`, having a dimension two variable named `v`, do the following:
+A variable (only one is allowed) is a finite dimensional vector or reals that will be *optimised* along with state and control values. To define an (almost empty!) optimal control problem, named `ocp`, having a dimension two variable named `v`, do the following:
 
 ```@example main
 using OptimalControl #hide
@@ -22,7 +22,7 @@ using OptimalControl #hide
 end
 ```
 
-Aliases `v₁` and `v₂` are automatically defined and can be used in subsequent expressions instead of `v[1]` and `v[2]`. The user can also define her own aliases for the components (one per dimension):
+Aliases `v₁` and `v₂` are automatically defined and can be used in subsequent expressions instead of `v[1]` and `v[2]`. The user can also define her own aliases for the components (one alias per dimension):
 
 ```@example main
 @def ocp begin
@@ -30,7 +30,7 @@ Aliases `v₁` and `v₂` are automatically defined and can be used in subsequen
 end
 ```
 
-A one dimensional can be declared according to
+A one dimensional variable can be declared according to
 
 ```@example main
 @def ocp begin
@@ -41,19 +41,20 @@ end
 ## Time
 
 ```julia
-   :( $t ∈ [$t0, $tf], time        ) 
+   :( $t ∈ [$t0, $tf], time ) 
 ```
 
 The independent variable or *time* is a scalar bound to a given interval. Its name is arbitrary.
 
 ```@example main
+t0 = 1
 tf = 5
 @def ocp begin
-    t ∈ [0, tf], time
+    t ∈ [t0, tf], time
 end
 ```
 
-One (or even the two bounds) can be variable, typically for minimum time problems:
+One (or even the two bounds) can be variable, typically for minimum time problems (see [Mayer cost](#mayer) section):
 
 ```@example main
 @def ocp begin
@@ -65,8 +66,8 @@ end
 ## [State](@id state)
 
 ```julia
-   :( $x ∈ R^$n, state               ) 
-   :( $x ∈ R   , state               ) 
+   :( $x ∈ R^$n, state ) 
+   :( $x ∈ R   , state ) 
 ```
 
 The state declaration defines the name and the dimension of the state:
@@ -88,8 +89,8 @@ end
 ## [Control](@id control)
 
 ```julia
-   :( $u ∈ R^$m, control             ) 
-   :( $u ∈ R   , control             ) 
+   :( $u ∈ R^$m, control ) 
+   :( $u ∈ R   , control ) 
 ```
 
 The control declaration defines the name and the dimension of the control:
@@ -111,13 +112,13 @@ end
 ## [Dynamics](@id dynamics)
 
 ```julia
-   :( ∂($x)($t) == $e1               ) 
+   :( ∂($x)($t) == $e1 ) 
 ```
 
 The dynamics is given in the standard vectorial ODE form:
 
 ```math
-    \dot{x}(t) = f(x(t), u(t)) \quad \text{or} \quad f(t, x(t), u(t))
+    \dot{x}(t) = f(x(t), u(t)) \quad \text{or} \quad \dot{x}(t) = f(t, x(t), u(t))
 ```
 
 depending on whether it is autonomous or not (the parser will detect dependence time, which entails that time and state must be declared prior to dynamics - an error will be issued otherwise). The symbol `∂` or the dotted state name
@@ -155,12 +156,13 @@ end
 
 F₀(x) = [x[2], 0]
 F₁(x) = [0, 1]
+nothing # hide
 ```
 
 !!! 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 [section ](#aliases) below):
+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 aliases) section below):
 
 ```@example main
 @def damped_integrator begin
@@ -177,11 +179,11 @@ end
 ## Constraints
 
 ```julia
-   :( $e1 == $e2                     ) 
-   :( $e1 ≤  $e2 ≤  $e3              ) 
-   :(        $e2 ≤  $e3              ) 
-   :( $e3 ≥  $e2 ≥  $e1              ) 
-   :( $e2 ≥  $e1                     ) 
+   :( $e1 == $e2        ) 
+   :( $e1 ≤  $e2 ≤  $e3 ) 
+   :(        $e2 ≤  $e3 ) 
+   :( $e3 ≥  $e2 ≥  $e1 ) 
+   :( $e2 ≥  $e1        ) 
 ```
 
 Constraints 
@@ -232,11 +234,11 @@ Line 7: 1 ≤ x₂(t)
 bad constraint declaration
 ```
 
-## Mayer cost
+## [Mayer cost](@id mayer)
 
 ```julia                                      
-   :( $e1                      → min ) 
-   :( $e1                      → max ) 
+   :( $e1 → min ) 
+   :( $e1 → max ) 
 ```
 
 Mayer costs are defined in a similar way to boundary conditions and follow the same rules. The symbol `→` is used
@@ -264,12 +266,12 @@ end
 ## Lagrange cost
 
 ```julia
-   :(             ∫($e1)       → min ) 
-   :(           - ∫($e1)       → min ) 
-   :(       $e1 * ∫($e2)       → min ) 
-   :(             ∫($e1)       → max ) 
-   :(           - ∫($e1)       → max ) 
-   :(       $e1 * ∫($e2)       → max ) 
+   :(       ∫($e1) → min ) 
+   :(     - ∫($e1) → min ) 
+   :( $e1 * ∫($e2) → min ) 
+   :(       ∫($e1) → max ) 
+   :(     - ∫($e1) → max ) 
+   :( $e1 * ∫($e2) → max ) 
 ```
 
 Lagrange (integral) costs are defined used the symbol `∫`, *with parenthesis:
@@ -279,13 +281,13 @@ Lagrange (integral) costs are defined used the symbol `∫`, *with parenthesis:
     t ∈ [0, 1], time
     x = (q, v) ∈ R², state
     u ∈ R, control
-    0.5∫( q(t) + u(t)^2 ) → min
+    0.5∫(q(t) + u(t)^2) → min
 end
 ```
 
 The integration range is implicitly equal to the time range, so the cost above is to be understood as
 ```math
-\int_0^1 q(t) + u^2(t)\,\mathrm{d}t \to \min.
+\int_0^1 (q(t) + u^2(t))\,\mathrm{d}t \to \min.
 ```
 
 As for the dynamics, the parser will detect whether the integrand depends or not on time (autonomous / non-autonomous case).
@@ -309,7 +311,7 @@ As for the dynamics, the parser will detect whether the integrand depends or not
    :(       $e2 * ∫($e3) + $e1 → max ) 
    :(             ∫($e2) - $e1 → max ) 
    :(       $e2 * ∫($e3) - $e1 → max ) 
- ```
+```
 
 Quite readily, Mayer and Lagrange costs can be combined into genral Bolza costs. For instance as follows:
 
@@ -319,7 +321,7 @@ Quite readily, Mayer and Lagrange costs can be combined into genral Bolza costs.
     t ∈ [t0, tf], time
     x = (q, v) ∈ R², state
     u ∈ R², control
-    (tf - t0) + 0.5∫( c(t) * u(t)^2 ) → min
+    (tf - t0) + 0.5∫(c(t) * u(t)^2) → min
 end
 ```
 

From 9e14dee53c2a3243d34c2cb2a1688b40434d2437 Mon Sep 17 00:00:00 2001
From: Jean-Baptiste Caillau <jean-baptiste.caillau@univ-cotedazur.fr>
Date: Sun, 4 Aug 2024 16:52:37 +0200
Subject: [PATCH 11/17] doc

---
 docs/Project.toml             | 1 +
 docs/src/tutorial-abstract.md | 9 +++++----
 2 files changed, 6 insertions(+), 4 deletions(-)

diff --git a/docs/Project.toml b/docs/Project.toml
index 94aa6a73..f4abba17 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -13,6 +13,7 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MINPACK = "4854310b-de5a-5eb6-a2a5-c1dee2bd17f9"
 NLPModelsIpopt = "f4238b75-b362-5c4c-b852-0801c9a21d71"
 NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
+OptimalControl = "5f98b655-cc9a-415a-b60e-744165666948"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 Percival = "01435c0c-c90d-11e9-3788-63660f8fbccc"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
diff --git a/docs/src/tutorial-abstract.md b/docs/src/tutorial-abstract.md
index 75e82bc3..7b04937b 100644
--- a/docs/src/tutorial-abstract.md
+++ b/docs/src/tutorial-abstract.md
@@ -162,7 +162,7 @@ nothing # hide
 !!! 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 aliases) section below):
+Currently, it is not possible to declare the dynamics component after component, but a simple workaround is to use *aliases* (check the relevant [aliases](#aliases) section below):
 
 ```@example main
 @def damped_integrator begin
@@ -372,7 +372,8 @@ end
 
 ```
 
-!!! caveat 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.
+!!! caveat
+    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
     You can use a trace mode for the macro `@def` to look at your code after expansions of the aliases adding `true` after your `begin ... end` block:
@@ -386,7 +387,7 @@ end
     q̇ = v(t)
     v̇ = u(t) - c(t)
     ẋ(t) == [q̇, v̇]
-end
+end true
 ```
 
 !!! caveat
@@ -412,7 +413,7 @@ forbidden alias name: (∂(x))(t)
 
 ## Misc
 
-- Declarations (variable - in any -, time, state and control) must be done first. Then, dynamics, constraints and cost can be introduced in an arbitrary order.
+- Declarations (variable - if any -, time, state and control) must be done first. Then, dynamics, constraints and cost can be introduced in an arbitrary order.
 - It is possible to provide numbers / labels (as math equations) for the constraints to improve readability (this is mostly for future use, typically to retrieve the Lagrange multiplier associated with the discretisation of a given constraint):
 
 ```@example main

From ac5a212a319fa7d14d4a4b71c8daf3fa54a5dbf7 Mon Sep 17 00:00:00 2001
From: Jean-Baptiste Caillau <jean-baptiste.caillau@univ-cotedazur.fr>
Date: Sun, 4 Aug 2024 17:01:48 +0200
Subject: [PATCH 12/17] doc

---
 docs/Project.toml             |  1 -
 docs/src/tutorial-abstract.md | 32 ++++++++++++++++++++------------
 2 files changed, 20 insertions(+), 13 deletions(-)

diff --git a/docs/Project.toml b/docs/Project.toml
index f4abba17..94aa6a73 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -13,7 +13,6 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MINPACK = "4854310b-de5a-5eb6-a2a5-c1dee2bd17f9"
 NLPModelsIpopt = "f4238b75-b362-5c4c-b852-0801c9a21d71"
 NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
-OptimalControl = "5f98b655-cc9a-415a-b60e-744165666948"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 Percival = "01435c0c-c90d-11e9-3788-63660f8fbccc"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
diff --git a/docs/src/tutorial-abstract.md b/docs/src/tutorial-abstract.md
index 7b04937b..5b9d412e 100644
--- a/docs/src/tutorial-abstract.md
+++ b/docs/src/tutorial-abstract.md
@@ -378,16 +378,24 @@ end
 !!! hint
     You can use a trace mode for the macro `@def` to look at your code after expansions of the aliases adding `true` after your `begin ... end` block:
 
-```@example main
-@def damped_integrator begin
-    tf ∈ R, variable
-    t ∈ [0, tf], time
-    x = (q, v) ∈ R², state
-    u ∈ R, control
-    q̇ = v(t)
-    v̇ = u(t) - c(t)
-    ẋ(t) == [q̇, v̇]
-end true
+```julia
+julia> @def damped_integrator begin
+           tf ∈ R, variable
+           t ∈ [0, tf], time
+           x = (q, v) ∈ R², state
+           u ∈ R, control
+           q̇ = v(t)
+           v̇ = u(t) - c(t)
+           ẋ(t) == [q̇, v̇]
+       end true
+
+variable: tf, dim: 1
+time: t, initial time: 0, final time: tf
+state: x, dim: 2
+control: u, dim: 1
+alias: q̇ = (x[2])(t)
+alias: v̇ = u(t) - c(t)
+dynamics: ẋ(t) == [(x[2])(t), u(t) - c(t)]
 ```
 
 !!! caveat
@@ -413,8 +421,8 @@ forbidden alias name: (∂(x))(t)
 
 ## Misc
 
-- Declarations (variable - if any -, time, state and control) must be done first. Then, dynamics, constraints and cost can be introduced in an arbitrary order.
-- It is possible to provide numbers / labels (as math equations) for the constraints to improve readability (this is mostly for future use, typically to retrieve the Lagrange multiplier associated with the discretisation of a given constraint):
+- Declarations (of variable - if any -, time, state and control) must be done first. Then, dynamics, constraints and cost can be introduced in an arbitrary order.
+- It is possible to provide numbers / labels (as in math equations) for the constraints to improve readability (this is mostly for future use, typically to retrieve the Lagrange multiplier associated with the discretisation of a given constraint):
 
 ```@example main
 @def damped_integrator begin

From 31da6e3a1e281c587fb969ed5ff6e4619809143d Mon Sep 17 00:00:00 2001
From: Jean-Baptiste Caillau <jean-baptiste.caillau@univ-cotedazur.fr>
Date: Sun, 4 Aug 2024 23:01:53 +0200
Subject: [PATCH 13/17] Update tutorial-basic-example.md

---
 docs/src/tutorial-basic-example.md | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/docs/src/tutorial-basic-example.md b/docs/src/tutorial-basic-example.md
index ea940443..2e76a30a 100644
--- a/docs/src/tutorial-basic-example.md
+++ b/docs/src/tutorial-basic-example.md
@@ -59,6 +59,8 @@ And plot the solution
 plot(sol)
 ```
 
+For a comprehensive introduction to the syntax used above to describe the optimal control problem, check [this tutorial(@ref abstract). 
+
 We can save the solution in a julia `.jld2` data file and reload it later, and also export a discretised version of the solution in a more portable [JSON](https://en.wikipedia.org/wiki/JSON) format.
 
 ```@example main

From 5aa1801904c8792f49732d6877e94082ac3bea84 Mon Sep 17 00:00:00 2001
From: Jean-Baptiste Caillau <jean-baptiste.caillau@univ-cotedazur.fr>
Date: Sun, 4 Aug 2024 23:02:57 +0200
Subject: [PATCH 14/17] Update tutorial-basic-example.md

---
 docs/src/tutorial-basic-example.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/docs/src/tutorial-basic-example.md b/docs/src/tutorial-basic-example.md
index 2e76a30a..8d1b14ef 100644
--- a/docs/src/tutorial-basic-example.md
+++ b/docs/src/tutorial-basic-example.md
@@ -59,7 +59,7 @@ And plot the solution
 plot(sol)
 ```
 
-For a comprehensive introduction to the syntax used above to describe the optimal control problem, check [this tutorial(@ref abstract). 
+For a comprehensive introduction to the syntax used above to describe the optimal control problem, check [this tutorial](@ref abstract). 
 
 We can save the solution in a julia `.jld2` data file and reload it later, and also export a discretised version of the solution in a more portable [JSON](https://en.wikipedia.org/wiki/JSON) format.
 

From c2a1db88a817dc5c16dcefe5b31f00e90cf38f8f Mon Sep 17 00:00:00 2001
From: Olivier Cots <olivier.cots@irit.fr>
Date: Sat, 3 Aug 2024 17:29:45 +0200
Subject: [PATCH 15/17] add MINPACK

---
 docs/Project.toml            |   1 +
 docs/src/tutorial-goddard.md | 126 +++++++++++++++++++++++++++++------
 docs/src/tutorial-iss.md     |  83 +++++++++++++++++++----
 src/solve.jl                 |   4 +-
 4 files changed, 178 insertions(+), 36 deletions(-)

diff --git a/docs/Project.toml b/docs/Project.toml
index 94aa6a73..2fb63658 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -3,6 +3,7 @@ CTBase = "54762871-cc72-4466-b8e8-f6c8b58076cd"
 CTDirect = "790bbbee-bee9-49ee-8912-a9de031322d5"
 CTFlows = "1c39547c-7794-42f7-af83-d98194f657c2"
 CommonSolve = "38540f10-b2f7-11e9-35d8-d573e4eb0ff2"
+DifferentiationInterface = "a0c0ee7d-e4b9-4e03-894e-1c5f64a51d63"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 DocumenterMermaid = "a078cd44-4d9c-4618-b545-3ab9d77f9177"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
diff --git a/docs/src/tutorial-goddard.md b/docs/src/tutorial-goddard.md
index 6b65c165..30bcdce7 100644
--- a/docs/src/tutorial-goddard.md
+++ b/docs/src/tutorial-goddard.md
@@ -33,25 +33,30 @@ $v(t) \leq v_{\max}$. The initial state is fixed while only the final mass is pr
     The Hamiltonian is affine with respect to the control, so singular arcs may occur,
     as well as constrained arcs due to the path constraint on the velocity (see below).
 
-## Direct method
-
 ```@setup main
 using Suppressor # to suppress warnings
 ```
 
-We import the `OptimalControl.jl` package to define the optimal control problem and 
-`NLPModelsIpopt.jl` to solve it. We import the `Plots.jl` package to plot the solution. The 
-`OrdinaryDiffEq.jl` package is used to define the shooting function for the indirect method and 
-the `NonlinearSolve.jl` package permits solve the shooting equation.
+We import the `OptimalControl.jl` package to define the optimal control problem and
+[`NLPModelsIpopt.jl`](https://github.com/JuliaSmoothOptimizers/NLPModelsIpopt.jl) to solve it. 
+We import the [`Plots.jl`](https://github.com/JuliaPlots/Plots.jl) package to plot the solution. 
+The [`OrdinaryDiffEq.jl`](https://github.com/SciML/OrdinaryDiffEq.jl) package is used to 
+define the shooting function for the indirect method and the 
+[`NonlinearSolve.jl`](https://github.com/SciML/NonlinearSolve.jl) and 
+[`MINPACK.jl`](https://github.com/sglyon/MINPACK.jl) packages permit to solve the shooting 
+equation.
 
 ```@example main
-using OptimalControl
-using NLPModelsIpopt
-using OrdinaryDiffEq # to get the Flow function
-using NonlinearSolve # NLE solver
-using Plots
+using OptimalControl  # to define the optimal control problem and more
+using NLPModelsIpopt  # to solve the problem via a direct method
+using OrdinaryDiffEq  # to get the Flow function from OptimalControl
+using NonlinearSolve  # interface to NLE solvers
+using MINPACK         # NLE solver: use to solve the shooting equation
+using Plots           # to plot the solution
 ```
 
+## Optimal control problem
+
 We define the problem
 
 ```@example main
@@ -100,6 +105,8 @@ end
 nothing # hide
 ```
 
+## Direct method
+
 We then solve it
 
 ```@example main
@@ -113,7 +120,7 @@ and plot the solution
 plt = plot(direct_sol, solution_label="(direct)", size=(800, 800))
 ```
 
-## Indirect method
+## Structure of the solution
 
 We first determine visually the structure of the optimal solution which is composed of a
 bang arc with maximal control, followed by a singular arc, then by a boundary arc and the final
@@ -216,6 +223,8 @@ fb = Flow(ocp, (x, p, tf) -> ub(x), (x, u, tf) -> g(x), (x, p, tf) -> μ(x, p))
 nothing # hide
 ```
 
+## Shooting function
+
 Then, we define the shooting function according to the optimal structure we have determined,
 that is a concatenation of four arcs.
 
@@ -240,6 +249,8 @@ end
 nothing # hide
 ```
 
+## Initial guess
+
 To solve the problem by an indirect shooting method, we then need a good initial guess,
 that is a good approximation of the initial costate, the three switching times and the
 final time.
@@ -266,22 +277,33 @@ s = similar(p0, 7)
 @suppress_err begin # hide
 shoot!(s, p0, t1, t2, t3, tf)
 end # hide
-println("Norm of the shooting function: ‖s‖ = ", norm(s), "\n")
+println("\nNorm of the shooting function: ‖s‖ = ", norm(s), "\n")
+```
+
+## Indirect shooting
+
+We aggregate the data to define the initial guess vector.
+
+```@example main
+ξ = [ p0 ; t1 ; t2 ; t3 ; tf ] # initial guess
 ```
 
+### NonlinearSolve.jl
+
 Finally, we can solve the shooting equations thanks to the [NonlinearSolve.jl](https://github.com/SciML/NonlinearSolve.jl).
 
 ```@example main
-nle = (s, ξ, λ) -> shoot!(s, ξ[1:3], ξ[4], ξ[5], ξ[6], ξ[7])   # auxiliary function
-                                                               # with aggregated inputs
-ξ = [ p0 ; t1 ; t2 ; t3 ; tf ]                                 # initial guess
+# auxiliary function with aggregated inputs
+nle! = (s, ξ, λ) -> shoot!(s, ξ[1:3], ξ[4], ξ[5], ξ[6], ξ[7])
+
+# NLE problem
+prob = NonlinearProblem(nle!, ξ)
 
-prob = NonlinearProblem(nle, ξ)
 global indirect_sol =      # hide
 @suppress_err begin # hide
 NonlinearSolve.solve(prob)      # hide
-# resolution of S(p0) = 0
-indirect_sol = NonlinearSolve.solve(prob; abstol=1e-8, reltol=1e-8, show_trace=Val(true))
+# resolution of S(ξ) = 0
+indirect_sol = solve(prob; abstol=1e-8, reltol=1e-8, show_trace=Val(true))
 end                 # hide
 
 # we retrieve the costate solution together with the times
@@ -291,6 +313,7 @@ t2 = indirect_sol.u[5]
 t3 = indirect_sol.u[6]
 tf = indirect_sol.u[7]
 
+println("")
 println("p0 = ", p0)
 println("t1 = ", t1)
 println("t2 = ", t2)
@@ -302,10 +325,71 @@ s = similar(p0, 7)
 @suppress_err begin # hide
 shoot!(s, p0, t1, t2, t3, tf)
 end # hide
-println("Norm of the shooting function: ‖s‖ = ", norm(s), "\n")
+println("\nNorm of the shooting function: ‖s‖ = ", norm(s), "\n")
 ```
 
-We plot the solution of the indirect solution (in red) over the solution of the direct method (in blue).
+### MINPACK.jl
+
+Instead of the [`NonlinearSolve.jl`](https://github.com/SciML/NonlinearSolve.jl) package we can use the [`MINPACK.jl`](https://github.com/sglyon/MINPACK.jl) package to solve 
+the shooting equation. To compute the Jacobian of the shooting function we use the 
+[`DifferentiationInterface.jl`](https://gdalle.github.io/DifferentiationInterface.jl/DifferentiationInterface) package with 
+[`ForwardDiff`](https://github.com/JuliaDiff/ForwardDiff.jl) backend.
+
+```@example main
+using DifferentiationInterface
+import ForwardDiff
+backend = AutoForwardDiff()
+nothing # hide
+```
+
+We are now in position to compute the Jacobian of the shooting function and use the `hybrj` solver
+from `MINPACK` through the `fsolve` function, providing the Jacobian.
+
+```@example main
+# auxiliary function with aggregated inputs
+nle! = (s, ξ) -> shoot!(s, ξ[1:3], ξ[4], ξ[5], ξ[6], ξ[7])
+
+# Jacobian of the shooting function
+jnle!(js, ξ) = jacobian!(nle!, js, backend, ξ)
+
+global indirect_sol =      # hide
+@suppress_err begin # hide
+fsolve(nle!, ξ, show_trace=false)      # hide
+# resolution of S(ξ) = 0
+try # hide
+indirect_sol = fsolve(nle!, jnle!, ξ, show_trace=true)
+catch e # hide
+println("hybrj not supported. Replaced by hybrd even if it is not visible on the doc.") # hide
+indirect_sol = fsolve(nle!, ξ, show_trace=true) # hide
+end # hide
+end                 # hide
+
+# we retrieve the costate solution together with the times
+p0 = indirect_sol.x[1:3]
+t1 = indirect_sol.x[4]
+t2 = indirect_sol.x[5]
+t3 = indirect_sol.x[6]
+tf = indirect_sol.x[7]
+
+println("")
+println("p0 = ", p0)
+println("t1 = ", t1)
+println("t2 = ", t2)
+println("t3 = ", t3)
+println("tf = ", tf)
+
+# Norm of the shooting function at solution
+s = similar(p0, 7)
+@suppress_err begin # hide
+shoot!(s, p0, t1, t2, t3, tf)
+end # hide
+println("\nNorm of the shooting function: ‖s‖ = ", norm(s), "\n")
+```
+
+## Plot of the solution
+
+We plot the solution of the indirect solution (in red) over the solution of the direct method 
+(in blue).
 
 ```@example main
 f = f1 * (t1, fs) * (t2, fb) * (t3, f0) # concatenation of the flows
diff --git a/docs/src/tutorial-iss.md b/docs/src/tutorial-iss.md
index 1a3d0085..83ded7a9 100644
--- a/docs/src/tutorial-iss.md
+++ b/docs/src/tutorial-iss.md
@@ -9,12 +9,15 @@ using Suppressor # to suppress warnings
 Let us start by importing the necessary packages.
 
 ```@example main
-using OptimalControl
-using OrdinaryDiffEq        # to get the Flow function from OptimalControl
-using NonlinearSolve        # NLE solver: we get the fsolve function
-using Plots
+using OptimalControl  # to define the optimal control problem and its flow
+using OrdinaryDiffEq  # to get the Flow function from OptimalControl
+using NonlinearSolve  # interface to NLE solvers
+using MINPACK         # NLE solver: use to solve the shooting equation
+using Plots           # to plot the solution
 ```
 
+## Optimal control problem
+
 Let us consider the following optimal control problem:
 
 ```math
@@ -53,6 +56,8 @@ end;
 nothing # hide
 ```
 
+## Boundary value problem
+
 The **pseudo-Hamiltonian** of this problem is
 
 ```math
@@ -84,6 +89,8 @@ where $[t]~=  (x(t),p(t),u(x(t), p(t)))$.
 
     Our goal is to solve this (BVP). Solving (BVP) consists in solving the Pontryagin Maximum Principle which provides necessary conditions of optimality.
 
+## Shooting function
+
 To achive our goal, let us first introduce the pseudo-Hamiltonian vector field
 
 ```math
@@ -135,31 +142,81 @@ Now, to solve the (BVP) we introduce the **shooting function**.
     \end{array}
 ```
 
-At the end, solving (BVP) is equivalent to solve $S(p_0) = 0$.
+```@example main
+S(p0) = π( φ(t0, x0, p0, tf) ) - xf # shooting function
+nothing # hide
+```
+
+## Resolution of the shooting equation
 
-This is what we call the **indirect simple shooting method**.
+At the end, solving (BVP) is equivalent to solve $S(p_0) = 0$. This is what we call the 
+**indirect simple shooting method**. We define an initial guess.
 
 ```@example main
-S(p0) = π( φ(t0, x0, p0, tf) ) - xf;                # shooting function
+ξ = [ 0.0 ] # initial guess
+nothing # hide
+```
 
-nle = (s, ξ, λ) -> s[1] = S(ξ[1])                   # auxiliary function
-ξ = [ 0.0 ]                                         # initial guess
+### NonlinearSolve.jl
 
-prob = NonlinearProblem(nle, ξ)
+```@example main
+nle! = (s, ξ, λ) -> s[1] = S(ξ[1])                  # auxiliary function
+prob = NonlinearProblem(nle!, ξ)                    # NLE problem
 global indirect_sol =      # hide
 @suppress_err begin # hide
 NonlinearSolve.solve(prob)      # hide
-indirect_sol = NonlinearSolve.solve(prob)           # resolution of S(p0) = 0
+indirect_sol = solve(prob, show_trace=Val(true))    # resolution of S(p0) = 0  
 end                 # hide
-
 p0_sol = indirect_sol.u[1]                          # costate solution
-println("costate:    p0 = ", p0_sol)
+println("\ncostate:    p0 = ", p0_sol)
 @suppress_err begin # hide
 println("shoot: |S(p0)| = ", abs(S(p0_sol)), "\n")
 end # hide
 nothing # hide
 ```
 
+### MINPACK.jl
+
+Instead of the [`NonlinearSolve.jl`](https://github.com/SciML/NonlinearSolve.jl) package we can use the [`MINPACK.jl`](https://github.com/sglyon/MINPACK.jl) package to solve 
+the shooting equation. To compute the Jacobian of the shooting function we use the 
+[`DifferentiationInterface.jl`](https://gdalle.github.io/DifferentiationInterface.jl/DifferentiationInterface) package with 
+[`ForwardDiff`](https://github.com/JuliaDiff/ForwardDiff.jl) backend.
+
+```@example main
+using DifferentiationInterface
+import ForwardDiff
+backend = AutoForwardDiff()
+nothing # hide
+```
+
+We are now in position to compute the Jacobian of the shooting function and use the `hybrj` solver
+from `MINPACK` through the `fsolve` function, providing the Jacobian.
+
+```@example main
+nle! = (s, ξ) -> s[1] = S(ξ[1])                        # auxiliary function
+jnle!(js, ξ) = jacobian!(nle!, js, backend, ξ)         # Jacobian of the shooting function
+global indirect_sol =      # hide
+@suppress_err begin # hide
+fsolve(nle!, ξ)      # hide
+try # hide
+indirect_sol = fsolve(nle!, jnle!, ξ, show_trace=true) # resolution of S(p0) = 0
+catch e # hide
+println("hybrj not supported. Replaced by hybrd even if it is not visible on the doc.") # hide
+indirect_sol = fsolve(nle!, ξ, show_trace=true) # hide
+end # hide
+end                 # hide
+
+# costate solution
+p0_sol = indirect_sol.x[1]
+println("\ncostate:    p0 = ", p0_sol)
+@suppress_err begin # hide
+println("shoot: |S(p0)| = ", abs(S(p0_sol)), "\n")
+end # hide
+nothing # hide
+```
+
+## Plot of the solution
+
 The solution can be plot calling first the flow.
 
 ```@example main
diff --git a/src/solve.jl b/src/solve.jl
index f84faffd..c7a39e88 100644
--- a/src/solve.jl
+++ b/src/solve.jl
@@ -1,6 +1,6 @@
 # --------------------------------------------------------------------------------------------------
 # Resolution
-import CommonSolve: solve
+import CommonSolve: solve, CommonSolve
 
 # by order of preference
 algorithms = ()
@@ -66,7 +66,7 @@ julia> sol = solve(ocp, init=(state=t->[-1+t, t*(t-1)], control=t->6-12*t))
 ```
 
 """
-function solve(ocp::OptimalControlModel, description::Symbol...;
+function CommonSolve.solve(ocp::OptimalControlModel, description::Symbol...;
     init=__ocp_init(),
     grid_size::Integer=CTDirect.__grid_size_direct(),
     display::Bool=CTDirect.__display(),

From 5faadddfdfddfd915cdeebb0773707e38f5b5ab7 Mon Sep 17 00:00:00 2001
From: Olivier Cots <66357348+ocots@users.noreply.github.com>
Date: Sun, 4 Aug 2024 13:34:06 +0000
Subject: [PATCH 16/17] foo

---
 docs/Project.toml            |   1 +
 docs/src/tutorial-goddard.md |  96 ++++++++++++++++++++++--------
 docs/src/tutorial-iss.md     | 111 +++++++++++++++++++++++++----------
 3 files changed, 155 insertions(+), 53 deletions(-)

diff --git a/docs/Project.toml b/docs/Project.toml
index 2fb63658..4533f990 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -1,4 +1,5 @@
 [deps]
+BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 CTBase = "54762871-cc72-4466-b8e8-f6c8b58076cd"
 CTDirect = "790bbbee-bee9-49ee-8912-a9de031322d5"
 CTFlows = "1c39547c-7794-42f7-af83-d98194f657c2"
diff --git a/docs/src/tutorial-goddard.md b/docs/src/tutorial-goddard.md
index 30bcdce7..a73f933d 100644
--- a/docs/src/tutorial-goddard.md
+++ b/docs/src/tutorial-goddard.md
@@ -290,21 +290,41 @@ We aggregate the data to define the initial guess vector.
 
 ### NonlinearSolve.jl
 
-Finally, we can solve the shooting equations thanks to the [NonlinearSolve.jl](https://github.com/SciML/NonlinearSolve.jl).
+We first use the [NonlinearSolve.jl](https://github.com/SciML/NonlinearSolve.jl) package to solve the shooting
+equation. Let us define the problem.
 
 ```@example main
 # auxiliary function with aggregated inputs
 nle! = (s, ξ, λ) -> shoot!(s, ξ[1:3], ξ[4], ξ[5], ξ[6], ξ[7])
 
-# NLE problem
+# NLE problem with initial guess
 prob = NonlinearProblem(nle!, ξ)
+nothing # hide
+```
 
-global indirect_sol =      # hide
-@suppress_err begin # hide
-NonlinearSolve.solve(prob)      # hide
+Let us do some benchmarking.
+
+```@example main
+using BenchmarkTools
+@benchmark solve(prob; abstol=1e-8, reltol=1e-8, show_trace=Val(false))
+```
+
+For small nonlinear systems, it could be faster to use the 
+[`SimpleNewtonRaphson()` descent algorithm](https://docs.sciml.ai/NonlinearSolve/stable/tutorials/code_optimization/).
+
+```@example main
+@benchmark solve(prob, SimpleNewtonRaphson(); abstol=1e-8, reltol=1e-8, show_trace=Val(false))
+```
+
+Now, let us solve the problem and retrieve the initial costate solution.
+
+```@example main
+global indirect_sol =               # hide
+@suppress_err begin                 # hide
+solve(prob; show_trace=Val(false))  # hide
 # resolution of S(ξ) = 0
 indirect_sol = solve(prob; abstol=1e-8, reltol=1e-8, show_trace=Val(true))
-end                 # hide
+end                                 # hide
 
 # we retrieve the costate solution together with the times
 p0 = indirect_sol.u[1:3]
@@ -330,7 +350,22 @@ println("\nNorm of the shooting function: ‖s‖ = ", norm(s), "\n")
 
 ### MINPACK.jl
 
-Instead of the [`NonlinearSolve.jl`](https://github.com/SciML/NonlinearSolve.jl) package we can use the [`MINPACK.jl`](https://github.com/sglyon/MINPACK.jl) package to solve 
+```@setup main
+using MINPACK
+function fsolve(f, j, x; kwargs...)
+    try
+        MINPACK.fsolve(f, j, x; kwargs...)
+    catch e
+        println("Erreur using MINPACK")
+        println(e)
+        println("hybrj not supported. Replaced by hybrd even if it is not visible on the doc.")
+        MINPACK.fsolve(f, x; kwargs...)
+    end
+end
+```
+
+Instead of the [`NonlinearSolve.jl`](https://github.com/SciML/NonlinearSolve.jl) package we can use the 
+[`MINPACK.jl`](https://github.com/sglyon/MINPACK.jl) package to solve 
 the shooting equation. To compute the Jacobian of the shooting function we use the 
 [`DifferentiationInterface.jl`](https://gdalle.github.io/DifferentiationInterface.jl/DifferentiationInterface) package with 
 [`ForwardDiff`](https://github.com/JuliaDiff/ForwardDiff.jl) backend.
@@ -342,27 +377,42 @@ backend = AutoForwardDiff()
 nothing # hide
 ```
 
-We are now in position to compute the Jacobian of the shooting function and use the `hybrj` solver
-from `MINPACK` through the `fsolve` function, providing the Jacobian.
+Let us define the problem to solve.
 
 ```@example main
 # auxiliary function with aggregated inputs
-nle! = (s, ξ) -> shoot!(s, ξ[1:3], ξ[4], ξ[5], ξ[6], ξ[7])
+nle!  = ( s, ξ) -> shoot!(s, ξ[1:3], ξ[4], ξ[5], ξ[6], ξ[7])
 
-# Jacobian of the shooting function
-jnle!(js, ξ) = jacobian!(nle!, js, backend, ξ)
+ # Jacobian of the (auxiliary) shooting function
+jnle! = (js, ξ) -> jacobian!(nle!, similar(ξ), js, backend, ξ)
+nothing # hide
+```
 
-global indirect_sol =      # hide
-@suppress_err begin # hide
-fsolve(nle!, ξ, show_trace=false)      # hide
+We are now in position to solve the problem with the `hybrj` solver from `MINPACK` through the `fsolve` 
+function, providing the Jacobian. Let us do some benchmarking.
+
+```@example main
+@benchmark fsolve(nle!, jnle!, ξ; show_trace=false) # initial guess given to the solver
+```
+
+We can also use the [preparation step](https://gdalle.github.io/DifferentiationInterface.jl/DifferentiationInterface/stable/tutorial1/#Preparing-for-multiple-gradients) of `DifferentiationInterface.jl`.
+
+```@example main
+extras = prepare_jacobian(nle!, similar(ξ), backend, ξ)
+jnle_prepared!(js, ξ) = jacobian!(nle!, similar(ξ), js, backend, ξ, extras)
+@benchmark fsolve(nle!, jnle_prepared!, ξ; show_trace=false)
+```
+
+Now, let us solve the problem and retrieve the initial costate solution.
+
+```@example main
+
+global indirect_sol =                                   # hide
+@suppress_err begin                                     # hide
+fsolve(nle!, jnle!, ξ; show_trace=false)                # hide
 # resolution of S(ξ) = 0
-try # hide
 indirect_sol = fsolve(nle!, jnle!, ξ, show_trace=true)
-catch e # hide
-println("hybrj not supported. Replaced by hybrd even if it is not visible on the doc.") # hide
-indirect_sol = fsolve(nle!, ξ, show_trace=true) # hide
-end # hide
-end                 # hide
+end                                                     # hide
 
 # we retrieve the costate solution together with the times
 p0 = indirect_sol.x[1:3]
@@ -380,9 +430,9 @@ println("tf = ", tf)
 
 # Norm of the shooting function at solution
 s = similar(p0, 7)
-@suppress_err begin # hide
+@suppress_err begin                                     # hide
 shoot!(s, p0, t1, t2, t3, tf)
-end # hide
+end                                                     # hide
 println("\nNorm of the shooting function: ‖s‖ = ", norm(s), "\n")
 ```
 
diff --git a/docs/src/tutorial-iss.md b/docs/src/tutorial-iss.md
index 83ded7a9..97789fe2 100644
--- a/docs/src/tutorial-iss.md
+++ b/docs/src/tutorial-iss.md
@@ -159,25 +159,63 @@ nothing # hide
 
 ### NonlinearSolve.jl
 
+We first use the [NonlinearSolve.jl](https://github.com/SciML/NonlinearSolve.jl) package to solve the shooting
+equation. Let us define the problem.
+
 ```@example main
-nle! = (s, ξ, λ) -> s[1] = S(ξ[1])                  # auxiliary function
-prob = NonlinearProblem(nle!, ξ)                    # NLE problem
-global indirect_sol =      # hide
-@suppress_err begin # hide
-NonlinearSolve.solve(prob)      # hide
-indirect_sol = solve(prob, show_trace=Val(true))    # resolution of S(p0) = 0  
-end                 # hide
+nle! = (s, ξ, λ) -> s[1] = S(ξ[1]) # auxiliary function
+prob = NonlinearProblem(nle!, ξ)   # NLE problem with initial guess
+nothing # hide
+```
+
+Let us do some benchmarking.
+
+```@example main
+using BenchmarkTools
+@benchmark solve(prob; show_trace=Val(false))
+```
+
+For small nonlinear systems, it could be faster to use the 
+[`SimpleNewtonRaphson()` descent algorithm](https://docs.sciml.ai/NonlinearSolve/stable/tutorials/code_optimization/).
+
+```@example main
+@benchmark solve(prob, SimpleNewtonRaphson(); show_trace=Val(false))
+```
+
+Now, let us solve the problem and retrieve the initial costate solution.
+
+```@example main
+global indirect_sol =                               # hide
+@suppress_err begin                                 # hide
+solve(prob; show_trace=Val(false))                  # hide
+indirect_sol = solve(prob; show_trace=Val(true))    # resolution of S(p0) = 0  
+end                                                 # hide
 p0_sol = indirect_sol.u[1]                          # costate solution
 println("\ncostate:    p0 = ", p0_sol)
-@suppress_err begin # hide
+@suppress_err begin                                 # hide
 println("shoot: |S(p0)| = ", abs(S(p0_sol)), "\n")
-end # hide
-nothing # hide
+end                                                 # hide
+nothing                                             # hide
 ```
 
 ### MINPACK.jl
 
-Instead of the [`NonlinearSolve.jl`](https://github.com/SciML/NonlinearSolve.jl) package we can use the [`MINPACK.jl`](https://github.com/sglyon/MINPACK.jl) package to solve 
+```@setup main
+using MINPACK
+function fsolve(f, j, x; kwargs...)
+    try
+        MINPACK.fsolve(f, j, x; kwargs...)
+    catch e
+        println("Erreur using MINPACK")
+        println(e)
+        println("hybrj not supported. Replaced by hybrd even if it is not visible on the doc.")
+        MINPACK.fsolve(f, x; kwargs...)
+    end
+end
+```
+
+Instead of the [`NonlinearSolve.jl`](https://github.com/SciML/NonlinearSolve.jl) package we can use the 
+[`MINPACK.jl`](https://github.com/sglyon/MINPACK.jl) package to solve 
 the shooting equation. To compute the Jacobian of the shooting function we use the 
 [`DifferentiationInterface.jl`](https://gdalle.github.io/DifferentiationInterface.jl/DifferentiationInterface) package with 
 [`ForwardDiff`](https://github.com/JuliaDiff/ForwardDiff.jl) backend.
@@ -189,30 +227,43 @@ backend = AutoForwardDiff()
 nothing # hide
 ```
 
-We are now in position to compute the Jacobian of the shooting function and use the `hybrj` solver
-from `MINPACK` through the `fsolve` function, providing the Jacobian.
+Let us define the problem to solve.
 
 ```@example main
-nle! = (s, ξ) -> s[1] = S(ξ[1])                        # auxiliary function
-jnle!(js, ξ) = jacobian!(nle!, js, backend, ξ)         # Jacobian of the shooting function
-global indirect_sol =      # hide
-@suppress_err begin # hide
-fsolve(nle!, ξ)      # hide
-try # hide
-indirect_sol = fsolve(nle!, jnle!, ξ, show_trace=true) # resolution of S(p0) = 0
-catch e # hide
-println("hybrj not supported. Replaced by hybrd even if it is not visible on the doc.") # hide
-indirect_sol = fsolve(nle!, ξ, show_trace=true) # hide
-end # hide
-end                 # hide
+nle!  = ( s, ξ) -> s[1] = S(ξ[1])                              # auxiliary function
+jnle! = (js, ξ) -> jacobian!(nle!, similar(ξ), js, backend, ξ) # Jacobian of nle
+nothing # hide
+```
+
+We are now in position to solve the problem with the `hybrj` solver from `MINPACK` through the `fsolve` 
+function, providing the Jacobian. Let us do some benchmarking.
+
+```@example main
+@benchmark fsolve(nle!, jnle!, ξ; show_trace=false) # initial guess given to the solver
+```
 
-# costate solution
-p0_sol = indirect_sol.x[1]
+We can also use the [preparation step](https://gdalle.github.io/DifferentiationInterface.jl/DifferentiationInterface/stable/tutorial1/#Preparing-for-multiple-gradients) of `DifferentiationInterface.jl`.
+
+```@example main
+extras = prepare_jacobian(nle!, similar(ξ), backend, ξ)
+jnle_prepared!(js, ξ) = jacobian!(nle!, similar(ξ), js, backend, ξ, extras)
+@benchmark fsolve(nle!, jnle_prepared!, ξ; show_trace=false)
+```
+
+Now, let us solve the problem and retrieve the initial costate solution.
+
+```@example main
+global indirect_sol =                                   # hide
+@suppress_err begin                                     # hide
+fsolve(nle!, jnle!, ξ; show_trace=false)                # hide
+indirect_sol = fsolve(nle!, jnle!, ξ; show_trace=true)  # resolution of S(p0) = 0
+end                                                     # hide
+p0_sol = indirect_sol.x[1]                              # costate solution
 println("\ncostate:    p0 = ", p0_sol)
-@suppress_err begin # hide
+@suppress_err begin                                     # hide
 println("shoot: |S(p0)| = ", abs(S(p0_sol)), "\n")
-end # hide
-nothing # hide
+end                                                     # hide
+nothing                                                 # hide
 ```
 
 ## Plot of the solution

From a47228b8b2aee3e56033a1df55cd038126950515 Mon Sep 17 00:00:00 2001
From: Olivier Cots <olivier.cots@irit.fr>
Date: Sun, 4 Aug 2024 23:07:54 +0200
Subject: [PATCH 17/17] restructuration

---
 docs/make.jl                         | 24 ++++---
 docs/src/tutorial-abstract.md        | 94 ++++++++++++++--------------
 docs/src/tutorial-basic-example-f.md |  2 +-
 docs/src/tutorial-basic-example.md   |  2 +-
 docs/src/tutorial-goddard.md         |  2 +-
 docs/src/tutorial-initial-guess.md   |  2 +-
 docs/src/tutorial-lqr-basic.md       |  2 +-
 docs/src/tutorial-plot.md            |  2 +-
 8 files changed, 67 insertions(+), 63 deletions(-)

diff --git a/docs/make.jl b/docs/make.jl
index d0188e24..4b9f172e 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -29,18 +29,22 @@ makedocs(
     ),
     pages = [
         "Introduction" => "index.md",
+        "Basic example" => [
+            "Energy min (abstract syntax)"   => "tutorial-basic-example.md", 
+            "Energy min (functional syntax)" => "tutorial-basic-example-f.md",
+            "Time minimisation"              => "tutorial-double-integrator.md",
+        ],
+        "Manual" => [
+            "Abstract syntax" => "tutorial-abstract.md",
+            "Initial guess"   => "tutorial-initial-guess.md",
+            "Plot a solution" => "tutorial-plot.md",
+        ],
         "Tutorials" => [
-            "tutorial-basic-example.md", 
-            "tutorial-basic-example-f.md", 
-            "tutorial-double-integrator.md",
-            "tutorial-initial-guess.md",
             "tutorial-continuation.md",
-            "tutorial-nlp.md",
-            "tutorial-plot.md",
-            "tutorial-lqr-basic.md",
+            "Goddard: direct, indirect"  => "tutorial-goddard.md",
             "tutorial-iss.md",
-            "tutorial-goddard.md",
-            "tutorial-abstract.md",
+            "Linear–quadratic regulator" => "tutorial-lqr-basic.md",
+            "tutorial-nlp.md",
         ],
         "API" => [
             "api-optimalcontrol.md",
@@ -51,7 +55,7 @@ makedocs(
             ],
         ],
         "Developers" => [
-            #"OptimalControl.jl" => "dev-optimalcontrol.md",
+            "OptimalControl.jl" => "dev-optimalcontrol.md",
             "Subpackages" => [
                 "CTBase.jl" => "dev-ctbase.md",
                 "CTDirect.jl" => "dev-ctdirect.md",
diff --git a/docs/src/tutorial-abstract.md b/docs/src/tutorial-abstract.md
index 5b9d412e..9cf401d4 100644
--- a/docs/src/tutorial-abstract.md
+++ b/docs/src/tutorial-abstract.md
@@ -1,4 +1,4 @@
-# [Abstract syntax](@id abstract)
+# [The abstract syntax to define an optimal control problem](@id abstract)
 
 The full grammar of OptimalControl.jl small *Domain Specific Language* is given below. The idea is to use a syntax that is
 - pure Julia (and, as such, effortlessly analysed by the standard Julia parser),
@@ -9,8 +9,8 @@ While the syntax will be transparent to those users familiar with Julia expressi
 ## [Variable](@id variable)
 
 ```julia
-   :( $v ∈ R^$q, variable ) 
-   :( $v ∈ R   , variable ) 
+:( $v ∈ R^$q, variable ) 
+:( $v ∈ R   , variable ) 
 ```
 
 A variable (only one is allowed) is a finite dimensional vector or reals that will be *optimised* along with state and control values. To define an (almost empty!) optimal control problem, named `ocp`, having a dimension two variable named `v`, do the following:
@@ -41,7 +41,7 @@ end
 ## Time
 
 ```julia
-   :( $t ∈ [$t0, $tf], time ) 
+:( $t ∈ [$t0, $tf], time ) 
 ```
 
 The independent variable or *time* is a scalar bound to a given interval. Its name is arbitrary.
@@ -66,15 +66,15 @@ end
 ## [State](@id state)
 
 ```julia
-   :( $x ∈ R^$n, state ) 
-   :( $x ∈ R   , state ) 
+:( $x ∈ R^$n, state ) 
+:( $x ∈ R   , state ) 
 ```
 
 The state declaration defines the name and the dimension of the state:
 
 ```@example main
 @def ocp begin
-       x ∈ R⁴, state
+    x ∈ R⁴, state
 end
 ```
 
@@ -82,15 +82,15 @@ As for the variable, there are automatic aliases (`x₁` for `x[1]`, *etc.*) and
 
 ```@example main
 @def ocp begin
-       x = (q₁, q₂, v₁, v₂) ∈ R⁴, state
+    x = (q₁, q₂, v₁, v₂) ∈ R⁴, state
 end
 ```
 
 ## [Control](@id control)
 
 ```julia
-   :( $u ∈ R^$m, control ) 
-   :( $u ∈ R   , control ) 
+:( $u ∈ R^$m, control ) 
+:( $u ∈ R   , control ) 
 ```
 
 The control declaration defines the name and the dimension of the control:
@@ -112,7 +112,7 @@ end
 ## [Dynamics](@id dynamics)
 
 ```julia
-   :( ∂($x)($t) == $e1 ) 
+:( ∂($x)($t) == $e1 ) 
 ```
 
 The dynamics is given in the standard vectorial ODE form:
@@ -179,11 +179,11 @@ end
 ## Constraints
 
 ```julia
-   :( $e1 == $e2        ) 
-   :( $e1 ≤  $e2 ≤  $e3 ) 
-   :(        $e2 ≤  $e3 ) 
-   :( $e3 ≥  $e2 ≥  $e1 ) 
-   :( $e2 ≥  $e1        ) 
+:( $e1 == $e2        ) 
+:( $e1 ≤  $e2 ≤  $e3 ) 
+:(        $e2 ≤  $e3 ) 
+:( $e3 ≥  $e2 ≥  $e1 ) 
+:( $e2 ≥  $e1        ) 
 ```
 
 Constraints 
@@ -211,7 +211,7 @@ In the example below, there are
     tf ≥ 0 
     x₂(t) ≤ 1
     u(t)^2 ≤ 1
-   end
+end
 ```
 
 !!! caveat
@@ -237,8 +237,8 @@ bad constraint declaration
 ## [Mayer cost](@id mayer)
 
 ```julia                                      
-   :( $e1 → min ) 
-   :( $e1 → max ) 
+:( $e1 → min ) 
+:( $e1 → max ) 
 ```
 
 Mayer costs are defined in a similar way to boundary conditions and follow the same rules. The symbol `→` is used
@@ -266,12 +266,12 @@ end
 ## Lagrange cost
 
 ```julia
-   :(       ∫($e1) → min ) 
-   :(     - ∫($e1) → min ) 
-   :( $e1 * ∫($e2) → min ) 
-   :(       ∫($e1) → max ) 
-   :(     - ∫($e1) → max ) 
-   :( $e1 * ∫($e2) → max ) 
+:(       ∫($e1) → min ) 
+:(     - ∫($e1) → min ) 
+:( $e1 * ∫($e2) → min ) 
+:(       ∫($e1) → max ) 
+:(     - ∫($e1) → max ) 
+:( $e1 * ∫($e2) → max ) 
 ```
 
 Lagrange (integral) costs are defined used the symbol `∫`, *with parenthesis:
@@ -287,7 +287,7 @@ end
 
 The integration range is implicitly equal to the time range, so the cost above is to be understood as
 ```math
-\int_0^1 (q(t) + u^2(t))\,\mathrm{d}t \to \min.
+\int_0^1 \left( q(t) + u^2(t) \right) \mathrm{d}t \to \min.
 ```
 
 As for the dynamics, the parser will detect whether the integrand depends or not on time (autonomous / non-autonomous case).
@@ -295,25 +295,25 @@ As for the dynamics, the parser will detect whether the integrand depends or not
 ## Bolza cost
 
 ```julia
-   :( $e1 +       ∫($e2)       → min ) 
-   :( $e1 + $e2 * ∫($e3)       → min ) 
-   :( $e1 -       ∫($e2)       → min ) 
-   :( $e1 - $e2 * ∫($e3)       → min ) 
-   :( $e1 +       ∫($e2)       → max ) 
-   :( $e1 + $e2 * ∫($e3)       → max ) 
-   :( $e1 -       ∫($e2)       → max ) 
-   :( $e1 - $e2 * ∫($e3)       → max ) 
-   :(             ∫($e2) + $e1 → min ) 
-   :(       $e2 * ∫($e3) + $e1 → min ) 
-   :(             ∫($e2) - $e1 → min ) 
-   :(       $e2 * ∫($e3) - $e1 → min ) 
-   :(             ∫($e2) + $e1 → max ) 
-   :(       $e2 * ∫($e3) + $e1 → max ) 
-   :(             ∫($e2) - $e1 → max ) 
-   :(       $e2 * ∫($e3) - $e1 → max ) 
-```
-
-Quite readily, Mayer and Lagrange costs can be combined into genral Bolza costs. For instance as follows:
+:( $e1 +       ∫($e2)       → min ) 
+:( $e1 + $e2 * ∫($e3)       → min ) 
+:( $e1 -       ∫($e2)       → min ) 
+:( $e1 - $e2 * ∫($e3)       → min ) 
+:( $e1 +       ∫($e2)       → max ) 
+:( $e1 + $e2 * ∫($e3)       → max ) 
+:( $e1 -       ∫($e2)       → max ) 
+:( $e1 - $e2 * ∫($e3)       → max ) 
+:(             ∫($e2) + $e1 → min ) 
+:(       $e2 * ∫($e3) + $e1 → min ) 
+:(             ∫($e2) - $e1 → min ) 
+:(       $e2 * ∫($e3) - $e1 → min ) 
+:(             ∫($e2) + $e1 → max ) 
+:(       $e2 * ∫($e3) + $e1 → max ) 
+:(             ∫($e2) - $e1 → max ) 
+:(       $e2 * ∫($e3) - $e1 → max ) 
+```
+
+Quite readily, Mayer and Lagrange costs can be combined into general Bolza costs. For instance as follows:
 
 ```@example main
 @def ocp begin
@@ -355,7 +355,7 @@ end
 ## [Aliases](@id aliases)
 
 ```julia
-        :( $a = $e1 )
+:( $a = $e1 )
 ```
 
 The single `=` symbol is used to define not a constraint but an alias, that is a purely syntactic replacement. There are some automatic aliases, *e.g.* `x₁` for `x[1]` if `x` is the state, and we have also seen that the user can define her own aliases when declaring the [variable](#variable), [state](#state) and [control](#control). Arbitrary aliases can be further defined, as below (compare with previous examples in the [dynamics](#dynamics) section):
@@ -399,7 +399,7 @@ dynamics: ẋ(t) == [(x[2])(t), u(t) - c(t)]
 ```
 
 !!! caveat
-    The dynamics of an OCP is indeed a particular constraint, be careful to use `==` and not a single `=` that would try define an alias:
+    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:
 
 ```julia
 julia> @def double_integrator begin
diff --git a/docs/src/tutorial-basic-example-f.md b/docs/src/tutorial-basic-example-f.md
index 43ed0a3c..794daa3a 100644
--- a/docs/src/tutorial-basic-example-f.md
+++ b/docs/src/tutorial-basic-example-f.md
@@ -1,4 +1,4 @@
-# [Basic example (functional version)](@id basic-f)
+# [Double integrator: energy min (functional syntax)](@id basic-f)
 
 Let us consider a wagon moving along a rail, whom acceleration can be controlled by a force $u$.
 We denote by $x = (x_1, x_2)$ the state of the wagon, that is its position $x_1$ and its velocity $x_2$.
diff --git a/docs/src/tutorial-basic-example.md b/docs/src/tutorial-basic-example.md
index 8d1b14ef..48db765e 100644
--- a/docs/src/tutorial-basic-example.md
+++ b/docs/src/tutorial-basic-example.md
@@ -1,4 +1,4 @@
-# [Basic example](@id basic)
+# [Double integrator: energy min (abstract syntax)](@id basic)
 
 Let us consider a wagon moving along a rail, whom acceleration can be controlled by a force $u$.
 We denote by $x = (x_1, x_2)$ the state of the wagon, that is its position $x_1$ and its velocity $x_2$.
diff --git a/docs/src/tutorial-goddard.md b/docs/src/tutorial-goddard.md
index a73f933d..b829438d 100644
--- a/docs/src/tutorial-goddard.md
+++ b/docs/src/tutorial-goddard.md
@@ -1,4 +1,4 @@
-# [Goddard problem](@id goddard)
+# [Direct and indirect methods for the Goddard problem](@id goddard)
 
 ## Introduction
 
diff --git a/docs/src/tutorial-initial-guess.md b/docs/src/tutorial-initial-guess.md
index 27c347a2..e6f473cb 100644
--- a/docs/src/tutorial-initial-guess.md
+++ b/docs/src/tutorial-initial-guess.md
@@ -1,4 +1,4 @@
-# Initial guess options
+# Initial guess (or iterate) for the resolution
 
 ```@meta
 CurrentModule =  OptimalControl
diff --git a/docs/src/tutorial-lqr-basic.md b/docs/src/tutorial-lqr-basic.md
index 050c07b0..292368fa 100644
--- a/docs/src/tutorial-lqr-basic.md
+++ b/docs/src/tutorial-lqr-basic.md
@@ -1,4 +1,4 @@
-# LQR example
+# A simple Linear–quadratic regulator example
 
 We consider the following Linear Quadratic Regulator (LQR) problem which consists in minimising
 
diff --git a/docs/src/tutorial-plot.md b/docs/src/tutorial-plot.md
index a6caf6dc..85527eca 100644
--- a/docs/src/tutorial-plot.md
+++ b/docs/src/tutorial-plot.md
@@ -1,4 +1,4 @@
-# Plot a solution
+# How to plot a solution
 
 In this tutorial we explain the different ways to plot a solution of an optimal control problem.