Refine discretisation tutorial

- Change disc_method parameter to scheme
- Remove const from parameters and wrap OCP in goddard_problem() function
- Add comments to dynamics functions explaining physics
- Add successful(sol) checks after solve calls
- Update AD backend section to use explicit solve mode with strategy instances
- Add documentation links using @extref syntax
- Simplify plotting code with enumerate and color indexing
- Add :external_cross_references to warnonly in make.jl

Author: ocots

docs/make.jl

@@ -55,12 +55,14 @@ draft = true  # Draft mode: if true, @example blocks in markdown are not execute
 # ═══════════════════════════════════════════════════════════════════════════════
 makedocs(;
     draft=draft,
-    warnonly=[:cross_references, :autodocs_block],
+    warnonly=[:cross_references, :autodocs_block, :external_cross_references],
     sitename="Tutorials",
     format=Documenter.HTML(;
         repolink="https://" * repo_url,
         prettyurls=false,
-        size_threshold_ignore=["tutorial-discretisation.md"],
+        size_threshold_ignore=[
+            "tutorial-discretisation.md",
+        ],
         assets=[
             asset("https://control-toolbox.org/assets/css/documentation.css"),
             asset("https://control-toolbox.org/assets/js/documentation.js"),

docs/src/tutorial-continuation.md

@@ -17,7 +17,7 @@ using Printf
 using Plots
 ```
 
-The `init` parameter of the `solve` function allows providing an initial guess. See the [initial guess documentation](https://control-toolbox.org/OptimalControl.jl/stable/manual-initial-guess.html) for more details.
+The `init` parameter of the `solve` function allows providing an initial guess. See the [initial guess documentation](@extref OptimalControl manual-initial-guess) for more details.
 
 We write a function that returns the OCP for a given final time:
 

docs/src/tutorial-discretisation.md

@@ -13,14 +13,16 @@ using NLPModels       # to retrieve information about the NLP problem
 using NLPModelsIpopt  # to solve the problem via a direct method
 using Plots           # to plot the solution
 
-const t0 = 0      # initial time
-const r0 = 1      # initial altitude
-const v0 = 0      # initial speed
-const m0 = 1      # initial mass
-const vmax = 0.1  # maximal authorized speed
-const mf = 0.6    # final mass to target
+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
+# Goddard problem function
+function goddard_problem()
+    ocp = @def begin # definition of the optimal control problem
 
     tf ∈ R, variable
     t ∈ [t0, tf], time
@@ -37,39 +39,46 @@ ocp = @def begin # definition of the optimal control problem
 
     r(tf) → max
 
-end;
+    end
+    return ocp
+end
+
+ocp = goddard_problem()
 
 # Dynamics
-const Cd = 310
-const Tmax = 3.5
-const β = 500
-const b = 2
+Cd = 310
+Tmax = 3.5
+β = 500
+b = 2
 
 F0(x) = begin
+    # Uncontrolled dynamics: gravity and drag
     r, v, m = x
-    D = Cd * v^2 * exp(-β*(r - 1)) # Drag force
-    return [ v, -D/m - 1/r^2, 0 ]
+    D = Cd * v^2 * exp(-β*(r - 1)) # Aerodynamic drag
+    return [ v, -D/m - 1/r^2, 0 ]   # [dr/dt, dv/dt, dm/dt]
 end
 
 F1(x) = begin
+    # Control dynamics: thrust contribution
     r, v, m = x
-    return [ 0, Tmax/m, -b*Tmax ]
+    return [ 0, Tmax/m, -b*Tmax ]   # [dr/dt, dv/dt, dm/dt] due to thrust
 end
 nothing # hide
 ```
 
 ## Discretisation schemes
 
-When calling `solve`, the option `disc_method=...` can be used to specify the discretization scheme. In addition to the default implicit `:midpoint` method, other available options include the implicit `:trapeze` method (also known as Crank–Nicolson), and Gauss–Legendre collocation schemes with 2 and 3 stages: `:gauss_legendre_2` and `:gauss_legendre_3`, of order 4 and 6 respectively. Note that higher-order methods typically result in larger NLP problems for the same number of time steps, and their accuracy also depends on the smoothness of the problem.
+When calling `solve`, the option `scheme=...` can be used to specify the discretization scheme. In addition to the default implicit `:midpoint` method, other available options include the implicit `:trapeze` method (also known as Crank–Nicolson), and Gauss–Legendre collocation schemes with 2 and 3 stages: `:gauss_legendre_2` and `:gauss_legendre_3`, of order 4 and 6 respectively. Note that higher-order methods typically result in larger NLP problems for the same number of time steps, and their accuracy also depends on the smoothness of the problem.
 
 Let us first solve the problem with the default `:midpoint` method and display the solution.
 
 ```@example main-disc
-sol = solve(ocp; disc_method=:midpoint, display=false)
+sol = solve(ocp; scheme=:midpoint, display=false)
+@assert successful(sol) "Solution failed"
 plot(sol; size=(800, 800))
 ```
 
-Let us now compare different discretization schemes to evaluate their accuracy and performance.
+Let us now compare different discretization schemes to evaluate their accuracy and performance. See the [Strategy options](@id manual-solve-strategy-options) for information on default values for parameters like `tol`.
 
 ```@example main-disc
 # Solve the problem with different discretization methods
@@ -89,8 +98,9 @@ schemes = [
     :gauss_legendre_3
 ]
 for scheme in schemes
-    bt = @btimed solve($ocp; disc_method=$scheme, tol=1e-8, display=false)
+    bt = @btimed solve($ocp; scheme=$scheme, tol=1e-8, display=false)
     local sol = bt.value
+    #@assert successful(sol) "Solution failed for scheme=$scheme"
     push!(solutions, (scheme, sol))
     push!(data, (
         Scheme=scheme, 
@@ -109,24 +119,28 @@ p_style = (legend=:none,)
 u_style = (legend=:topright,)
 styles = (state_style=x_style, costate_style=p_style, control_style=u_style)
 
-scheme, sol = solutions[1]
-plt = plot(sol; label=string(scheme), styles...)
-for (scheme, sol) in solutions[2:end]
-    plt = plot!(sol; label=string(scheme), styles...)
+plt = plot()
+for (i, (scheme, sol)) in enumerate(solutions)
+    plt = plot!(sol; label=string(scheme), color=i, styles...)
 end
 plot(plt; size=(800, 800))
 ```
 
 ## Large problems and AD backend
 
-For some large problems, you may notice that the solver takes a long time before the iterations actually begin. This is due to the computation of sparse derivatives — specifically, the Jacobian of the constraints and the Hessian of the Lagrangian — which can be quite costly. One possible alternative is to set the option `adnlp_backend=:manual`, which uses simpler sparsity patterns. The resulting matrices are faster to compute but are also less sparse, so this represents a trade-off between automatic differentiation (AD) preparation time and the efficiency of the optimization itself.
+For some large problems, you may notice that the solver takes a long time before the iterations actually begin. This is due to the computation of sparse derivatives — specifically, the Jacobian of the constraints and the Hessian of the Lagrangian — which can be quite costly. One possible alternative is to set the option `backend=:manual` in the modeler, which uses simpler sparsity patterns. The resulting matrices are faster to compute but are also less sparse, so this represents a trade-off between automatic differentiation (AD) preparation time and the efficiency of the optimization itself.
 
 ```@example main-disc
-solve(ocp; disc_method=:gauss_legendre_3, grid_size=1000, adnlp_backend=:manual)
+disc = OptimalControl.Collocation(grid_size=1000, scheme=:gauss_legendre_3)
+mod = OptimalControl.ADNLP(backend=:manual)
+sol = OptimalControl.Ipopt()
+solve(ocp; discretizer=disc, modeler=mod, solver=sol)
 nothing # hide
 ```
 
-Let us now compare the performance of the two backends. We will use the `@btimed` macro from the `BenchmarkTools` package to measure the time taken for both the preparation of the NLP problem and the execution of the solver. Besides, we will also collect the number of non-zero elements in the Jacobian and Hessian matrices, which can be useful to understand the sparsity of the problem, thanks to the functions `get_nnzo`, `get_nnzj`, and `get_nnzj` from the `NLPModels` package. The problem is first discretized with the `direct_transcription` method and then solved with the `ipopt` solver, see the [tutorial on direct transcription](@ref tutorial-nlp) for more details.
+This explicit approach with strategy instances is less high-level than the descriptive method used earlier. It provides more control over the solving process and is detailed in the [explicit mode documentation](@extref OptimalControl manual-solve-explicit).
+
+Let us now compare the performance of the two backends. We will use the `@btimed` macro from the `BenchmarkTools` package to measure the time taken for both the preparation of the NLP problem and the execution of the solver. Separating these measurements allows us to understand whether the bottleneck is in the discretization/modeling phase (which depends on the AD backend) or in the solver phase (which depends on the sparsity of the matrices). We use a lower-level approach with explicit strategy instances to separately measure the discretization/modeling time and the solver time.
 
 ```@example main-disc
 # DataFrame to store the results
@@ -147,30 +161,32 @@ backends = [:optimized, :manual]
 # Loop over the backends
 for adnlp_backend in backends
 
-    # Discretize the problem with a large grid size and Gauss-Legendre method
-    bt = @btimed direct_transcription($ocp; 
-        disc_method=:gauss_legendre_3, 
-        grid_size=1000, 
-        adnlp_backend=$adnlp_backend,
-    )
-    docp = bt.value
-    nlp = nlp_model(docp)
+    # Create strategy instances
+    discretizer = OptimalControl.Collocation(grid_size=1000, scheme=:gauss_legendre_3)
+    modeler = OptimalControl.ADNLP(backend=adnlp_backend)
+    solver = OptimalControl.Ipopt(print_level=0)
+
+    # Initial guess
+    init = OptimalControl.build_initial_guess(ocp, nothing)
+
+    # Discretize and build NLP model
+    # This lower-level approach with discretize/nlp_model/solve is detailed in the [tutorial on direct transcription](@ref tutorial-nlp)
+    bt = @btimed begin
+        docp = discretize($ocp, $discretizer)
+        nlp_model(docp, $init, $modeler)
+    end
     prepa_time = bt.time
+    nlp = bt.value
 
     # Get the number of non-zero elements
     nnzo = get_nnzo(nlp) # Gradient of the Objective
     nnzj = get_nnzj(nlp) # Jacobian of the constraints
     nnzh = get_nnzh(nlp) # Hessian of the Lagrangian
 
     # Solve the problem
-    bt = @btimed ipopt($nlp; 
-        print_level=0, 
-        mu_strategy="adaptive", 
-        tol=1e-8,
-        sb="yes",
-    )
-    nlp_sol = bt.value
+    bt = @btimed solve($nlp, $solver; display=false)
     exec_time = bt.time
+    nlp_sol = bt.value
 
     # Store the results in the DataFrame
     push!(data,