add tutos

Author: ocots

docs/Project.toml

@@ -1,5 +1,28 @@
 [deps]
+DifferentiationInterface = "a0c0ee7d-e4b9-4e03-894e-1c5f64a51d63"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+MINPACK = "4854310b-de5a-5eb6-a2a5-c1dee2bd17f9"
+MadNLP = "2621e9c9-9eb4-46b1-8089-e8c72242dfb6"
+NLPModelsIpopt = "f4238b75-b362-5c4c-b852-0801c9a21d71"
+NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
+OptimalControl = "5f98b655-cc9a-415a-b60e-744165666948"
+OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+Suppressor = "fd094767-a336-5f1f-9728-57cf17d0bbfb"
 
 [compat]
-Documenter = "1"
+DifferentiationInterface = "0.6"
+Documenter = "1.8"
+ForwardDiff = "0.10"
+LinearAlgebra = "1"
+MINPACK = "1.3"
+MadNLP = "0.8"
+NLPModelsIpopt = "0.10"
+NonlinearSolve = "4.6"
+OptimalControl = "1.0"
+OrdinaryDiffEq = "6.93"
+Plots = "1.40"
+Suppressor = "0.2"
+julia = "1"

docs/make.jl

@@ -1,21 +1,34 @@
 using Documenter
+using OptimalControl
 
 repo_url = "github.com/control-toolbox/Tutorials.jl"
 
 makedocs(;
-    remotes=nothing,
-    warnonly=:cross_references,
-    sitename="Tutorials",
+    warnonly=[:cross_references, :autodocs_block],
+    sitename="Tutorials.jl",
     format=Documenter.HTML(;
         repolink="https://" * repo_url,
-        prettyurls=false,
-        size_threshold_ignore=["index.md"],
+        prettyurls=true,
+        size_threshold_ignore=[
+            ""
+        ],
         assets=[
             asset("https://control-toolbox.org/assets/css/documentation.css"),
             asset("https://control-toolbox.org/assets/js/documentation.js"),
         ],
     ),
-    pages=["Introduction" => "index.md"],
+    pages=[
+        "Getting Started" => "index.md",
+        "Tutorials and Advanced Features" => [
+            "Discrete continuation" => "tutorial-continuation.md",
+            "Discretisation methods" => "tutorial-discretisation.md",
+            "NLP manipulations" => "tutorial-nlp.md",
+            "Indirect simple shooting" => "tutorial-iss.md",
+            "Goddard: direct, indirect" => "tutorial-goddard.md",
+            "Linear–quadratic regulator" => "tutorial-lqr-basic.md",
+            "Minimal action" => "tutorial-mam.md",
+        ],
+    ],
 )
 
-deploydocs(; repo=repo_url * ".git", devbranch="main")
+deploydocs(; repo=repo_url * ".git", devbranch="main")

docs/src/assets/Goddard_and_Rocket.jpg

@@ -1,12 +1,69 @@
 # Tutorials
 
-Documentation for [Tutorials](https://github.com/control-toolbox/Tutorials.jl).
+This collection of tutorials is part of the [control-toolbox ecosystem](https://github.com/control-toolbox). The control-toolbox ecosystem gathers Julia packages for mathematical control and applications. It aims to provide tools to model and solve optimal control problems with ordinary differential equations by direct and indirect methods. If you want to define an optimal control problem and solve it, please check the [documentation](https://control-toolbox.org/OptimalControl.jl).
 
-## Dependencies
+From this page, you can find a list of tutorials to solve optimal control problems with OptimalControl.
 
-All the numerical simulations to generate this documentation are performed with the following packages.
+## Reproducibility
+
+```@raw html
+<details><summary>The documentation of this package was built using these direct dependencies,</summary>
+```
+
+```@example
+using Pkg # hide
+Pkg.status() # hide
+```
+
+```@raw html
+</details>
+```
+
+```@raw html
+<details><summary>and using this machine and Julia version.</summary>
+```
 
 ```@example
-using Pkg
-Pkg.status()
+using InteractiveUtils # hide
+versioninfo() # hide
 ```
+
+```@raw html
+</details>
+```
+
+```@raw html
+<details><summary>A more complete overview of all dependencies and their versions is also provided.</summary>
+```
+
+```@example
+using Pkg # hide
+Pkg.status(; mode = PKGMODE_MANIFEST) # hide
+```
+
+```@raw html
+</details>
+```
+
+```@eval
+using TOML
+using Markdown
+version = TOML.parse(read("../../Project.toml", String))["version"]
+name = TOML.parse(read("../../Project.toml", String))["name"]
+link_manifest = "https://github.com/SciML/" *
+                name *
+                ".jl/tree/gh-pages/v" *
+                version *
+                "/assets/Manifest.toml"
+link_project = "https://github.com/SciML/" *
+               name *
+               ".jl/tree/gh-pages/v" *
+               version *
+               "/assets/Project.toml"
+Markdown.parse("""You can also download the
+[manifest]($link_manifest)
+file and the
+[project]($link_project)
+file.
+""")
+```

docs/src/index.md

@@ -0,0 +1,141 @@
+# Discrete continuation
+
+```@meta
+CurrentModule =  OptimalControl
+```
+
+Using the warm start option, it is easy to implement a basic discrete continuation method, where a sequence of problems is solved using each solution as initial guess for the next problem.
+This usually gives better and faster convergence than solving each problem with the same initial guess, and is a way to handle problems that require a good initial guess.
+
+
+## Continuation on parametric OCP
+
+The most compact syntax to perform a discrete continuation is to use a function that returns the OCP for a given value of the continuation parameter, and solve a sequence of these problems. We illustrate this on a very basic double integrator with increasing fixed final time.
+
+First we load the required packages
+
+```@example main
+using OptimalControl
+using NLPModelsIpopt
+using Printf
+using Plots
+```
+
+and write a function that returns the OCP for a given final time
+
+```@example main
+function ocp_T(T)
+    ocp = @def begin
+        t ∈ [0, T], time
+        x ∈ R², state
+        u ∈ R, control
+        q = x₁
+        v = x₂
+        q(0) == 0
+        v(0) == 0
+        q(T) == 1
+        v(T) == 0
+        ẋ(t) == [ v(t), u(t) ]
+        ∫(u(t)^2) → min
+    end
+    return ocp
+end
+nothing # hide
+```
+
+Then we perform the continuation with a simple *for* loop, using each solution to initialize the next problem.
+
+```@example main
+init1 = ()
+for T=1:5
+    ocp1 = ocp_T(T) 
+    sol1 = solve(ocp1; display=false, init=init1)
+    global init1 = sol1
+    @printf("T %.2f objective %9.6f iterations %d\n", T, objective(sol1), iterations(sol1))
+end
+```
+
+## Continuation on global variable
+
+As a second example, we show how to avoid redefining a new OCP each time, and modify the original one instead.
+More precisely we now solve a Goddard problem for a decreasing maximal thrust. If we store the value for *Tmax* in a global variable, we can simply modify this variable and keep the same OCP problem during the continuation.
+
+Let us first define the Goddard problem (note that the formulation below illustrates all the possible constraints types, and the problem could be defined in a more compact way).
+
+```@example main
+Cd = 310
+Tmax = 3.5
+β = 500
+b = 2
+function F0(x)
+    r, v, m = x
+    D = Cd * v^2 * exp(-β*(r - 1))
+    return [ v, -D/m - 1/r^2, 0 ]
+end
+function F1(x)
+    r, v, m = x
+    return [ 0, Tmax/m, -b*Tmax ]
+end
+r0 = 1
+v0 = 0
+m0 = 1
+mf = 0.6
+x0 = [r0, v0, m0]
+vmax = 0.1
+
+@def ocp begin
+    tf ∈ R, variable
+    t ∈ [0, tf], time
+    x ∈ R^3, state
+    u ∈ R, control
+    0.01 ≤ tf ≤ Inf
+    r = x[1]
+    v = x[2]
+    m = x[3]
+    x(0) == x0
+    m(tf) == mf
+    r0 ≤ r(t) ≤ r0 + 0.1
+    v0 ≤ v(t) ≤ vmax
+    mf ≤ m(t) ≤ m0
+    0 ≤ u(t) ≤ 1
+    ẋ(t) == F0(x(t)) + u(t) * F1(x(t))
+    r(tf) → max
+end
+
+sol0 = solve(ocp; display=false)
+@printf("Objective for reference solution %.6f\n", objective(sol0))
+```
+
+Then we perform the continuation on the maximal thrust.
+
+```@example main
+sol       = sol0
+Tmax_list = []
+obj_list  = []
+for Tmax_local=3.5:-0.5:1
+    global Tmax = Tmax_local  
+    global sol = solve(ocp; display=false, init=sol)
+    @printf("Tmax %.2f objective %.6f iterations %d\n", Tmax, objective(sol), iterations(sol))
+    push!(Tmax_list, Tmax)
+    push!(obj_list, objective(sol))
+end 
+```
+
+We plot now the objective w.r.t the maximal thrust, as well as both solutions for *Tmax*=3.5 and *Tmax*=1.
+
+```@example main
+using Plots.PlotMeasures # for leftmargin
+
+plt_obj = plot(Tmax_list, obj_list;
+    seriestype=:scatter,
+    title="Goddard problem",
+    label="r(tf)", 
+    xlabel="Maximal thrust (Tmax)",
+    ylabel="Maximal altitude r(tf)")
+
+plt_sol = plot(sol0; solution_label="(Tmax = "*string(Tmax_list[1])*")")
+plot!(plt_sol, sol;  solution_label="(Tmax = "*string(Tmax_list[end])*")")
+
+layout = grid(2, 1, heights=[0.2, 0.8])
+plot(plt_obj, plt_sol; layout=layout, size=(800, 1000), leftmargin=5mm)
+```

docs/src/tutorial-continuation.md

@@ -0,0 +1,98 @@
+# [Discretisation methods](@id tutorial-discretisation-methods)
+
+## Discretisation formulas
+When calling `solve`, the option `disc_method=...` can be used to set the discretisation scheme.
+In addition to the default implicit `:trapeze` method (aka Crank-Nicolson), other choices are available, namely implicit `:midpoint` and the Gauss-Legendre collocations with 2 and  stages, `:gauss_legendre_2` and `:gauss_legendre_3`, of order 4 and 6 respectively. 
+Note that higher order methods will typically lead to larger NLP problems for the same number of time steps, and that accuracy will also depend on the smoothness of the problem.
+
+As an example we will use the [Goddard problem](@ref tutorial-goddard)
+```@example main
+using OptimalControl  # to define the optimal control problem and more
+using NLPModelsIpopt  # to solve the problem via a direct method
+using Plots           # to plot the solution
+
+t0 = 0      # initial time
+r0 = 1      # initial altitude
+v0 = 0      # initial speed
+m0 = 1      # initial mass
+vmax = 0.1  # maximal authorized speed
+mf = 0.6    # final mass to target
+
+ocp = @def begin # definition of the optimal control problem
+
+    tf ∈ R, variable
+    t ∈ [t0, tf], time
+    x = (r, v, m) ∈ R³, state
+    u ∈ R, control
+
+    x(t0) == [ r0, v0, m0 ]
+    m(tf) == mf,         (1)
+    0 ≤ u(t) ≤ 1
+    r(t) ≥ r0
+    0 ≤ v(t) ≤ vmax
+
+    ẋ(t) == F0(x(t)) + u(t) * F1(x(t))
+
+    r(tf) → max
+
+end;
+
+# Dynamics
+const Cd = 310
+const Tmax = 3.5
+const β = 500
+const b = 2
+
+F0(x) = begin
+    r, v, m = x
+    D = Cd * v^2 * exp(-β*(r - 1)) # Drag force
+    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
+```
+Now let us compare different discretisations
+```@example main
+sol_trapeze = solve(ocp; tol=1e-8)
+plot(sol_trapeze)
+
+sol_midpoint = solve(ocp, disc_method=:midpoint; tol=1e-8)
+plot!(sol_midpoint)
+
+sol_euler = solve(ocp, disc_method=:euler; tol=1e-8)
+plot!(sol_euler)
+
+sol_euler_imp = solve(ocp, disc_method=:euler_implicit; tol=1e-8)
+plot!(sol_euler_imp)
+
+sol_gl2 = solve(ocp, disc_method=:gauss_legendre_2; tol=1e-8)
+plot!(sol_gl2)
+
+sol_gl3 = solve(ocp, disc_method=:gauss_legendre_3; tol=1e-8)
+plot!(sol_gl3)
+```
+
+## Large problems and AD backend
+For some large problems, you may notice that solving spends a long time before the iterations actually begin.
+This is due to the computing of the sparse derivatives, namely the Jacobian of the constraints and the Hessian of the Lagrangian, that can become quite costly.
+A possible alternative is to set the option `adnlp_backend=:manual`, which will use more basic sparsity patterns.
+The resulting matrices are faster to compute but are also less sparse, so this is a trade-off bewteen the AD preparation and the optimization itself.
+
+```@example main
+solve(ocp, disc_method=:gauss_legendre_3, grid_size=1000, adnlp_backend=:manual)
+nothing # hide
+```
+
+## Explicit time grid
+The option `time_grid=...` allows to pass the complete time grid vector `t0, t1, ..., tf`, which is typically useful if one wants a non uniform grid. 
+In the case of a free initial and/or final time, provide a normalised grid between 0 and 1. 
+Note that `time_grid` will override `grid_size` if both are present.
+
+```@example main
+sol = solve(ocp, time_grid=[0, 0.1, 0.5, 0.9, 1], display=false)
+println(time_grid(sol))
+```