up tuto disc

Author: ocots

docs/Project.toml

@@ -6,21 +6,23 @@ ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MINPACK = "4854310b-de5a-5eb6-a2a5-c1dee2bd17f9"
 MadNLP = "2621e9c9-9eb4-46b1-8089-e8c72242dfb6"
+NLPModels = "a4795742-8479-5a88-8948-cc11e1c8c1a6"
 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"
-Tutorials = "cb10daa6-a5e5-4c25-a171-ae181b8ea3c9"
 
 [compat]
+DataFrames = "1.7"
 DifferentiationInterface = "0.7"
 Documenter = "1.8"
 ForwardDiff = "0.10, 1"
 LinearAlgebra = "1"
 MINPACK = "1.3"
 MadNLP = "0.8"
+NLPModels = "0.21"
 NLPModelsIpopt = "0.10"
 NonlinearSolve = "4.6"
 OptimalControl = "1.0"

docs/extras/disc.jl

@@ -0,0 +1,131 @@
+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
+    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
+Cd = 310
+Tmax = 3.5
+β = 500
+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
+
+# Solve the problem with different discretization methods
+solutions = []
+schemes = [:trapeze, :midpoint, :euler, :euler_implicit, :gauss_legendre_2, :gauss_legendre_3]
+for scheme in schemes
+    sol = solve(ocp; disc_method=scheme, tol=1e-8, display=false)
+    push!(solutions, (scheme, sol))
+end
+
+# Plot the results
+x_style = (legend=:none,)
+p_style = (legend=:none,)
+styles = (state_style=x_style, costate_style=p_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...)
+end
+plot(plt; size=(800, 800))
+
+# Benchmark the problem with different AD backends
+using NLPModels
+using DataFrames
+using BenchmarkTools
+
+# DataFrame to store the results
+data = DataFrame(
+    Backend=Symbol[],
+    NNZO=Int[],
+    NNZJ=Int[],
+    NNZH=Int[],
+    PrepTime=Float64[],
+    ExecTime=Float64[],
+    Objective=Float64[], 
+    Iterations=Int[]
+)
+
+# The different AD backends to test
+backends = [:optimized, :manual]
+
+# Loop over the backends
+for adnlp_backend in backends
+
+    println("Testing backend: ", adnlp_backend)
+
+    # Discretize the problem with a large grid size and Gauss-Legendre method
+    bt = @btimed direct_transcription($ocp; 
+        disc_method=:euler, 
+        grid_size=1000, 
+        adnlp_backend=$adnlp_backend,
+    )
+    docp, nlp = bt.value
+    prepa_time = bt.time
+
+    # 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
+    exec_time = bt.time
+
+    # Store the results in the DataFrame
+    push!(data, 
+        (
+            Backend=adnlp_backend, 
+            NNZO=nnzo, 
+            NNZJ=nnzj, 
+            NNZH=nnzh, 
+            PrepTime=prepa_time, 
+            ExecTime=exec_time, 
+            Objective=-nlp_sol.objective, 
+            Iterations=nlp_sol.iter
+        )
+    )
+end
+println(data)

docs/src/assets/Manifest.toml

@@ -2,7 +2,7 @@
 
 julia_version = "1.11.5"
 manifest_format = "2.0"
-project_hash = "32b2794783d26475326d2bd2122ce29a5a993624"
+project_hash = "a0f7a38986e41d62db81da09b6b12da3c238b679"
 
 [[deps.ADNLPModels]]
 deps = ["ADTypes", "ForwardDiff", "LinearAlgebra", "NLPModels", "Requires", "ReverseDiff", "SparseArrays", "SparseConnectivityTracer", "SparseMatrixColorings"]

docs/src/assets/Project.toml

@@ -6,6 +6,7 @@ ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MINPACK = "4854310b-de5a-5eb6-a2a5-c1dee2bd17f9"
 MadNLP = "2621e9c9-9eb4-46b1-8089-e8c72242dfb6"
+NLPModels = "a4795742-8479-5a88-8948-cc11e1c8c1a6"
 NLPModelsIpopt = "f4238b75-b362-5c4c-b852-0801c9a21d71"
 NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
 OptimalControl = "5f98b655-cc9a-415a-b60e-744165666948"
@@ -15,12 +16,14 @@ Suppressor = "fd094767-a336-5f1f-9728-57cf17d0bbfb"
 Tutorials = "cb10daa6-a5e5-4c25-a171-ae181b8ea3c9"
 
 [compat]
+DataFrames = "1.7"
 DifferentiationInterface = "0.7"
 Documenter = "1.8"
 ForwardDiff = "0.10, 1"
 LinearAlgebra = "1"
 MINPACK = "1.3"
 MadNLP = "0.8"
+NLPModels = "0.21"
 NLPModelsIpopt = "0.10"
 NonlinearSolve = "4.6"
 OptimalControl = "1.0"

docs/src/tutorial-continuation.md

@@ -1,10 +1,5 @@
 # Discrete continuation
 
-```@meta
-CurrentModule =  OptimalControl
-Draft = false
-```
-
 By using the warm start option, it is easy to implement a basic discrete continuation method, in which a sequence of problems is solved by using each solution as the initial guess for the next problem.
 This approach typically leads to faster and more reliable convergence than solving each problem with the same initial guess and is particularly useful for problems that require a good initial guess to converge.