Merge pull request #141 from control-toolbox/update-doc

docs and CTDirect update

Author: ocots

Project.toml

@@ -11,6 +11,6 @@ DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 
 [compat]
 CTBase = "0.7"
-CTDirect = "0.3"
+CTDirect = "0.4"
 CTFlows = "0.3"
-julia = "1.8"
+julia = "1.9"

docs/Project.toml

@@ -1,2 +1,3 @@
 [deps]
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+OptimalControl = "5f98b655-cc9a-415a-b60e-744165666948"

docs/make.jl

@@ -6,7 +6,9 @@ makedocs(
     format = Documenter.HTML(prettyurls = false),
     pages = [
         "Introduction" => "index.md",
-        "API" => "api.md"
+        "Tutorials" => ["basic-example.md", "goddard.md"],
+        "API" => "api.md",
+        "Developpers" => "dev-api.md"
     ]
 )
 

docs/src/basic-example.md

@@ -0,0 +1,60 @@
+# Basic example
+
+Consider we want to minimise the cost functional
+
+```math
+    \frac{1}{2}\int_{0}^{1} u^2(t) \, \mathrm{d}t
+```
+
+subject to the dynamical constraints for $t \in [0, 1]$
+
+```math
+    \dot x_1(t) = x_2(t), \quad \dot x_2(t) = u(t) \in \mathbb{R},
+```
+
+and the limit conditions
+
+```math
+    x(0) = (-1, 0), \quad x(1) = (0, 0).
+```
+
+First, we need to import the `OptimalControl.jl` package:
+
+```@example main
+using OptimalControl
+```
+
+Then, we can define the problem
+
+```@example main
+t0 = 0
+tf = 1
+A = [ 0 1
+      0 0 ]
+B = [ 0
+      1 ]
+
+@def ocp_di begin
+    t ∈ [ t0, tf ], time                         # time interval
+    x ∈ R², state                                # state
+    u ∈ R, control                               # control
+    x(t0) == [-1, 0],    (initial_con)           # initial condition
+    x(tf) == [0, 0],    (final_con)              # final condition
+    ẋ(t) == A * x(t) + B * u(t)                  # dynamics
+    ∫( 0.5u(t)^2 ) → min                         # objective
+end
+nothing # hide
+```
+
+Solve it
+
+```@example main
+sol_di = solve(ocp_di)
+nothing # hide
+```
+
+and plot the solution
+
+```@example main
+plot(sol_di, size=(700, 700))
+```

docs/src/dev-api.md

@@ -0,0 +1,11 @@
+# Internal functions
+
+```@meta
+CurrentModule = OptimalControl 
+```
+
+```@autodocs
+Modules = [OptimalControl]
+Order = [:module, :type, :function, :macro]
+Public = false
+```

docs/src/goddard.md

@@ -0,0 +1,105 @@
+# Advanced example
+
+This well-known problem[^1] [^2] models the ascent of a rocket through the atmosphere, and we restrict here ourselves to vertical (one dimensional) trajectories. The state variables are the altitude $r$, speed $v$ and mass $m$ of the rocket during the flight, for a total dimension of 3. The rocket is subject to gravity $g$, thrust $u$ and drag force $D$ (function of speed and altitude). The final time $T$ is free, and the objective is to reach a maximal altitude with a bounded fuel consumption.
+
+We thus want to solve the optimal control problem in Mayer form
+
+```math
+    \max\, r(T)
+```
+
+subject to the control dynamics
+
+```math
+    \dot{r} = v, \quad
+    \dot{v} = \frac{T_{\max}\,u - D(r,v)}{m} - g, \quad
+    \dot{m} = -u,
+```
+
+and subject to the control constraint $u(t) \in [0,1]$ and the state constraint
+$v(t) \leq v_{\max}$. The initial state is fixed while only the final mass is prescribed.
+
+!!! note
+
+    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).
+
+We import the `OptimalControl.jl` package:
+
+```@example main
+using OptimalControl
+```
+
+We define the problem
+
+```@example main
+# Parameters
+const Cd = 310
+const Tmax = 3.5
+const β = 500
+const b = 2
+
+t0 = 0
+r0 = 1
+v0 = 0
+vmax = 0.1
+m0 = 1
+mf = 0.6
+
+# Initial state
+x0 = [ r0, v0, m0 ]
+
+# Abstract model
+@def ocp_goddard begin
+
+    tf, variable
+    t ∈ [ t0, tf ], time
+    x ∈ R³, state
+    u ∈ R, control
+    
+    r = x₁
+    v = x₂
+    m = x₃
+   
+    x(t0) == [ r0, v0, m0 ]
+    0  ≤ u(t) ≤ 1
+         r(t) ≥ r0,     (1)
+    0  ≤ v(t) ≤ vmax,   (2)
+    mf ≤ m(t) ≤ m0,     (3)
+
+    ẋ(t) == F0(x(t)) + u(t) * F1(x(t))
+ 
+    r(tf) → max
+    
+end;
+
+F0(x) = begin
+    r, v, m = x
+    D = Cd * v^2 * exp(-β*(r - 1))
+    return [ v, -D/m - 1/r^2, 0 ]
+end
+
+F1(x) = begin
+    r, v, m = x
+    return [ 0, Tmax/m, -b*Tmax ]
+end
+nothing # hide
+```
+
+Solve it
+
+```@example main
+N = 50
+direct_sol_goddard = solve(ocp_goddard, grid_size=N)
+nothing # hide
+```
+
+and plot the solution
+
+```@example main
+plot(direct_sol_goddard, size=(700, 700))
+```
+
+[^1]: R.H. Goddard. A Method of Reaching Extreme Altitudes, volume 71(2) of Smithsonian Miscellaneous Collections. Smithsonian institution, City of Washington, 1919.
+
+[^2]: H. Seywald and E.M. Cliff. Goddard problem in presence of a dynamic pressure limit. Journal of Guidance, Control, and Dynamics, 16(4):776–781, 1993.

docs/src/index.md

@@ -1,6 +1,17 @@
-# Introduction to the OptimalControl.jl package
+# OptimalControl.jl
 
-The `OptimalControl.jl` package is part of the [control-toolbox ecosystem](https://github.com/control-toolbox). It aims to provide tools to solve optimal control problems by direct and indirect methods. An optimal control problem can be described as minimising the cost functional
+```@meta
+CurrentModule =  OptimalControl
+```
+
+The `OptimalControl.jl` package is part of the [control-toolbox ecosystem](https://github.com/control-toolbox).
+
+!!! note "Install"
+
+    To install a package from the control-toolbox ecosystem, 
+    please visit the [installation page](https://github.com/control-toolbox#installation).
+
+This package aims to provide tools to solve optimal control problems by direct and indirect methods. An optimal control problem can be described as minimising the cost functional
 
 ```math
 g(t_0, x(t_0), t_f, x(t_f)) + \int_{t_0}^{t_f} f^{0}(t, x(t), u(t))~\mathrm{d}t
@@ -23,75 +34,3 @@ and other constraints such as
 \phi_l &\le& \phi(t_0, x(t_0), t_f, x(t_f)) &\le& \phi_u.
 \end{array}
 ```
-
-**Contents.**
-
-```@contents
-Pages = ["index.md", "api.md"]
-Depth = 2
-```
-
-## Installation
-
-To install a package from the control-toolbox ecosystem, please visit the [installation page](https://github.com/control-toolbox#installation).
-
-## Basic usage
-
-Consider we want to minimise the cost functional
-
-```math
-    \frac{1}{2}\int_{0}^{1} u^2(t) \, \mathrm{d}t
-```
-
-subject to the dynamical constraints for $t \in [0, 1]$
-
-```math
-    \dot x_1(t) = x_2(t), \quad \dot x_2(t) = u(t) \in \mathbb{R},
-```
-
-and the limit conditions
-
-```math
-    x(0) = (-1, 0), \quad x(1) = (0, 0).
-```
-
-First, we need to import the `OptimalControl.jl` package:
-
-```@example main
-using OptimalControl
-```
-
-Then, we can define the problem
-
-```@example main
-ocp = Model()                                    # the model for the problem definition
-
-state!(ocp, 2)                                   # dimension of the state
-control!(ocp, 1)                                 # dimension of the control
-time!(ocp, [0, 1])                               # time interval
-
-objective!(ocp, :lagrange, (x, u) -> 0.5u^2)     # objective
-
-A = [ 0 1
-      0 0 ]
-B = [ 0
-      1 ]
-dynamics!(ocp, (x, u) -> A*x + B*u) # dynamics
-
-constraint!(ocp, :initial, [-1, 0])              # initial condition
-constraint!(ocp, :final,   [0, 0])               # final condition
-nothing # hide
-```
-
-Solve it
-
-```@example main
-sol = solve(ocp)
-nothing # hide
-```
-
-and plot the solution
-
-```@example main
-plot(sol, size=(700, 700))
-```

src/solve.jl

@@ -22,15 +22,15 @@ function solve(ocp::OptimalControlModel, description::Symbol...;
     method = getFullDescription(description, algorithmes)
 
     # todo: OptimalControlInit must be in CTBase
-    #=
+    
     if isnothing(init)
         init =  OptimalControlInit()
     elseif init isa CTBase.OptimalControlSolution
         init = OptimalControlInit(init)
     else
-        OptimalControlInit(init...)
+        init = OptimalControlInit(x_init=init[rg(1,ocp.state_dimension)],u_init=init[rg(ocp.state_dimension+1,ocp.state_dimension+ocp.control_dimension)],v_init=init[rg(ocp.state_dimension+ocp.control_dimension+1,lastindex(init))])
     end
-    =#
+    
 
     # print chosen method
     display ? println("\nMethod = ", method) : nothing
@@ -41,6 +41,8 @@ function solve(ocp::OptimalControlModel, description::Symbol...;
     end
 end
 
+rg(i,j) = i == j ? i : i:j
+
 function clean(d::Description)
     return d\(:direct, )
 end

test/test_goddard_indirect.jl

@@ -39,9 +39,9 @@ u1 = 1
 # singular control
 H0 = Lift(F0)
 H1 = Lift(F1)
-H01  = @Poisson {H0, H1}
-H001 = @Poisson {H0, H01}
-H101 = @Poisson {H1, H01}
+H01  = @Lie {H0, H1}
+H001 = @Lie {H0, H01}
+H101 = @Lie {H1, H01}
 us(x, p) = -H001(x, p) / H101(x, p)
 
 # boundary control