Merge pull request #292 from control-toolbox/253-document-algorithms

Add tuto solve

docs/make.jl

@@ -37,9 +37,10 @@ makedocs(
         ],
         "Manual" => [
             "Abstract syntax" => "tutorial-abstract.md",
-            "Flow"            => "tutorial-flow.md",
+            "Solve"           => "tutorial-solve.md",
             "Initial guess"   => "tutorial-initial-guess.md",
             "Plot a solution" => "tutorial-plot.md",
+            "Flow"            => "tutorial-flow.md",
         ],
         "Tutorials" => [
             "tutorial-continuation.md",

docs/src/tutorial-flow.md

@@ -1,6 +1,6 @@
 # [How to compute flows](@id manual-flow)
 
-In this tutorial, we explain the `Flow` function from `OptimalControl` package.
+In this tutorial, we explain the `Flow` function from `OptimalControl.jl` package.
 
 ## Basic usage
 

docs/src/tutorial-initial-guess.md

@@ -1,4 +1,4 @@
-# Initial guess (or iterate) for the resolution
+# [Initial guess (or iterate) for the resolution](@id tutorial-init)
 
 ```@meta
 CurrentModule =  OptimalControl

docs/src/tutorial-iss.md

@@ -5,11 +5,11 @@ In this tutorial we present the indirect simple shooting method on a simple exam
 Let us start by importing the necessary packages.
 
 ```@example main
-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
+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
@@ -129,7 +129,7 @@ nothing # hide
     ```
     where $\mathbf{H}(z) = H(z, u(z))$ and $\vec{\mathbf{H}} = (\nabla_p \mathbf{H}, -\nabla_x \mathbf{H})$. This is what is actually computed by `Flow`.
 
-Now, to solve the (BVP) we introduce the **shooting function**.
+Now, to solve the (BVP) we introduce the **shooting function**:
 
 ```math
     \begin{array}{rlll}
@@ -139,7 +139,7 @@ Now, to solve the (BVP) we introduce the **shooting function**.
 ```
 
 ```@example main
-S(p0) = π( φ(t0, x0, p0, tf) ) - xf # shooting function
+S(p0) = π( φ(t0, x0, p0, tf) ) - xf    # shooting function
 nothing # hide
 ```
 
@@ -149,7 +149,7 @@ At the end, solving (BVP) is equivalent to solve $S(p_0) = 0$. This is what we c
 **indirect simple shooting method**. We define an initial guess.
 
 ```@example main
-ξ = [ 0.0 ] # initial guess
+ξ = [ 0.1 ]    # initial guess
 nothing # hide
 ```
 
@@ -159,8 +159,8 @@ We first use the [NonlinearSolve.jl](https://github.com/SciML/NonlinearSolve.jl)
 equation. Let us define the problem.
 
 ```@example main
-nle! = (s, ξ, λ) -> s[1] = S(ξ[1]) # auxiliary function
-prob = NonlinearProblem(nle!, ξ)   # NLE problem with initial guess
+nle! = (s, ξ, λ) -> s[1] = S(ξ[1])    # auxiliary function
+prob = NonlinearProblem(nle!, ξ)      # NLE problem with initial guess
 nothing # hide
 ```
 
@@ -181,8 +181,8 @@ For small nonlinear systems, it could be faster to use the
 Now, let us solve the problem and retrieve the initial costate solution.
 
 ```@example main
-indirect_sol = solve(prob; show_trace=Val(true))    # resolution of S(p0) = 0  
-p0_sol = indirect_sol.u[1]                          # costate solution
+indirect_sol = solve(prob; show_trace=Val(true))      # resolution of S(p0) = 0  
+p0_sol = indirect_sol.u[1]                            # costate solution
 println("\ncostate:    p0 = ", p0_sol)
 println("shoot: |S(p0)| = ", abs(S(p0_sol)), "\n")
 ```
@@ -219,16 +219,16 @@ nothing # hide
 Let us define the problem to solve.
 
 ```@example main
-nle!  = ( s, ξ) -> s[1] = S(ξ[1])                              # auxiliary function
-jnle! = (js, ξ) -> jacobian!(nle!, similar(ξ), js, backend, ξ) # Jacobian of nle
+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
+@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`.
@@ -242,8 +242,8 @@ jnle_prepared!(js, ξ) = jacobian!(nle!, similar(ξ), js, backend, ξ, extras)
 Now, let us solve the problem and retrieve the initial costate solution.
 
 ```@example main
-indirect_sol = fsolve(nle!, jnle!, ξ; show_trace=true)  # resolution of S(p0) = 0
-p0_sol = indirect_sol.x[1]                              # costate solution
+indirect_sol = fsolve(nle!, jnle!, ξ; show_trace=true)    # resolution of S(p0) = 0
+p0_sol = indirect_sol.x[1]                                # costate solution
 println("\ncostate:    p0 = ", p0_sol)
 println("shoot: |S(p0)| = ", abs(S(p0_sol)), "\n")
 ```
@@ -253,8 +253,8 @@ println("shoot: |S(p0)| = ", abs(S(p0_sol)), "\n")
 The solution can be plot calling first the flow.
 
 ```@example main
-sol = φ( (t0, tf), x0, p0_sol )
-plot(sol, size=(800, 600))
+sol = φ((t0, tf), x0, p0_sol)
+plot(sol)
 ```
 
 In the indirect shooting method, the research of the optimal control is replaced by the computation

docs/src/tutorial-nlp.md

@@ -61,14 +61,15 @@ For a first example we use the `ipopt` solver from [`NLPModelsIpopt.jl`](https:/
 
 ```@example main
 using NLPModelsIpopt
+
 nlp_sol = ipopt(nlp; print_level=5, mu_strategy="adaptive", tol=1e-8, sb="yes")
 nothing # hide
 ```
 
 Then we can rebuild and plot an optimal control problem solution (note that the multipliers are optional, but the OCP costate will not be retrieved if the multipliers are not provided).
 
 ```@example main
-sol = OptimalControlSolution(docp, primal=nlp_sol.solution, dual=nlp_sol.multipliers)
+sol = OptimalControlSolution(docp; primal=nlp_sol.solution, dual=nlp_sol.multipliers)
 plot(sol)
 ```
 ## Change the NLP solver
@@ -86,15 +87,15 @@ Another  possible NLP solver is [`Percival.jl`](https://jso.dev/Percival.jl).
 ```@example main
 using Percival
 
-nlp_sol = percival(nlp, verbose=1)
+nlp_sol = percival(nlp; verbose=1)
 ```
 
 ## Initial guess
 
 An initial guess, including warm start, can be passed to [`direct_transcription`](@ref) the same way as for `solve`.
 
 ```@example main
-docp, nlp = direct_transcription(ocp, init=sol)
+docp, nlp = direct_transcription(ocp; init=sol)
 nothing # hide
 ```
 

docs/src/tutorial-solve.md

@@ -0,0 +1,187 @@
+# [The solve function](@id manual-solve)
+
+In this tutorial, we explain the [`solve`](@ref) function from `OptimalControl.jl` package.
+
+## Basic usage
+
+Les us define a basic optimal control problem.
+
+```@example main
+using OptimalControl
+
+t0 = 0
+tf = 1
+x0 = [-1, 0]
+
+ocp = @def begin
+
+    t ∈ [ t0, tf ], time
+    x = (q, v) ∈ R², state
+    u ∈ R, control
+
+    x(t0) == x0
+    x(tf) == [ 0, 0 ]
+
+    ẋ(t)  == [ v(t), u(t) ]
+
+    ∫( 0.5u(t)^2 ) → min
+
+end
+nothing # hide
+```
+
+Let us try to solve the problem:
+
+```@setup main_repl
+using OptimalControl
+t0 = 0
+tf = 1
+x0 = [-1, 0]
+ocp = @def begin
+    t ∈ [ t0, tf ], time
+    x = (q, v) ∈ R², state
+    u ∈ R, control
+    x(t0) == x0
+    x(tf) == [ 0, 0 ]
+    ẋ(t)  == [ v(t), u(t) ]
+    ∫( 0.5u(t)^2 ) → min
+end
+```
+
+```@repl main_repl
+solve(ocp)
+```
+
+As you can see, an error occured since we need the package [`NLPModelsIpopt.jl`](https://jso.dev/NLPModelsIpopt.jl).
+
+Actually, the default solving method is what we call a 
+[direct method](https://en.wikipedia.org/wiki/Optimal_control#Numerical_methods_for_optimal_control). 
+In a direct method, the optimal control problem is transcribed to a nonlinear optimization problem (NLP) of the form
+
+```math
+\text{minimize}\quad F(y), \quad\text{subject to the constraints}\quad g(y)=0, \quad h(y)\le 0. 
+```
+
+The `OptimalControl.jl` package makes the transcription but it needs a package to modelise the NLP problem and 
+another one to solve it. The `NLPModelsIpopt.jl` package provides an interface to the well-known solver 
+[`Ipopt`](https://coin-or.github.io/Ipopt/) that can be used to solve general nonlinear programming problems.
+
+```@example main
+using NLPModelsIpopt
+
+solve(ocp)
+nothing # hide
+```
+
+## Solvers
+
+`OptimalControl.jl` offers a list of methods to solve your optimal control problem. To get the list of methods,
+simply call `available_methods`.
+
+```@example main
+available_methods()
+```
+
+Each line is a method. The priority is given from top to bottom. This means that 
+
+```julia
+solve(ocp)
+```
+
+is equivalent to 
+
+```julia
+solve(ocp, :direct, :adnlp, :ipopt)
+```
+
+Let us detail the three symbols. As you can see, there are only direct methods. The symbol `:adnlp` is for the 
+choice of modeler. As said before, the NLP problem needs to be modelised in `Julia` code. We use 
+[`ADNLPModels.jl`](https://jso.dev/ADNLPModels.jl) which provides automatic differentiation (AD)-based model implementations that conform to the [`NLPModels.jl`](https://github.com/JuliaSmoothOptimizers/ADNLPModels.jl) API.
+
+The last symbol is what distinguish the two available methods. It corresponds to the NLP solver. By default, we
+use the `NLPModelsIpopt.jl` package but you can choose [`MadNLP.jl`](https://madnlp.github.io/MadNLP.jl) which
+is an open-source nonlinear programming solver, purely implemented in Julia and which implements a filter 
+line-search interior-point algorithm, as used in Ipopt.
+
+Note that you are not compelled to provide the full description of the method to use it. A partial description,
+if not ambiguous, will work. Just remember that if several full descriptions contain the partial one, then, the 
+priority is given to first one in the list. Hence, these calls are all equivalent:
+
+```julia
+solve(ocp)
+solve(ocp, :direct                )
+solve(ocp,          :adnlp        )
+solve(ocp,                  :ipopt)
+solve(ocp, :direct, :adnlp        )
+solve(ocp, :direct,         :ipopt)
+solve(ocp, :direct, :adnlp, :ipopt)
+```
+
+Let us try `MadNLP.jl`.
+
+```@example main
+using MadNLP
+
+solve(ocp, :direct, :adnlp, :madnlp; display=true)
+nothing # hide
+```
+
+Again, there are several equivalent manners to use `MadNLP.jl`.
+
+```julia
+solve(ocp,                  :madnlp)
+solve(ocp, :direct,         :madnlp)
+solve(ocp,          :adnlp, :madnlp)
+solve(ocp, :direct, :adnlp, :madnlp)
+```
+
+## Options
+
+### Direct method
+
+The options common to all the direct methods may be seen directly in the 
+[code](https://github.com/control-toolbox/CTDirect.jl/blob/23c76374fb91f65f8d10c03b7d8ff5c75a5f4076/src/solve.jl#L68).
+There are `init`, `grid_size` and `time_grid`. 
+
+- The `init` option can be used to set an initial guess for the solver. See this [tutorial](@ref tutorial-init). 
+- The `grid_size` option corresponds to the size of the (uniform) time discretization grid. 
+More precisely, it is the number of steps, that is if `N = grid_size` and if the initial and final times are 
+denoted respectively `t0` and `tf`, then we have:
+```julia
+Δt = (tf - t0) / N
+```
+- The `time_grid` option is the grid of times: `t0, t1, ..., tf`. If the initial and/or the final times are free, then you can provide a normalised grid between 0 and 1. Note that you can set either `grid_size` or `time_grid` but not both.
+
+```@example main
+sol = solve(ocp; grid_size=10, display=false)
+sol.times
+```
+
+Or with `MadNLP`:
+
+```@example main
+sol = solve(ocp, :madnlp; grid_size=10, display=false)
+sol.times
+```
+
+### NLPModelsIpopt
+
+You can provide any option of `NLPModelsIpopt` or `Ipopt` with a pair `keyword=value`. Please check
+the list of [`Ipopt` options](https://coin-or.github.io/Ipopt/OPTIONS.html) and the 
+[`NLPModelsIpopt.jl` documentation](https://jso.dev/NLPModelsIpopt.jl).
+
+```@example main
+solve(ocp, :ipopt; max_iter=0)
+nothing # hide
+```
+
+### MadNLP
+
+If you use the `MadNLP` solver, then you can provide any option of it. Please check the
+[`MadNLP.jl` documentation](https://madnlp.github.io/MadNLP.jl) and the list of 
+[`MadNLP.jl`options](https://madnlp.github.io/MadNLP.jl/stable/options/).
+
+```@example main
+solve(ocp, :madnlp; max_iter=1, display=true)
+nothing # hide
+```