Refactor canonical tests with separated modelers/solvers

test/helpers/print_utils.jl

@@ -110,7 +110,7 @@ function print_test_header(show_memory::Bool=false)
     print(" │ ")
     printstyled(rpad("Type", 4); bold=true)
     print(" │ ")
-    printstyled(rpad("Problem", 8); bold=true)
+    printstyled(rpad("Problem", 14); bold=true)
     print(" │ ")
     printstyled(rpad("Disc", 8); bold=true)
     print(" │ ")
@@ -140,7 +140,7 @@ function print_test_header(show_memory::Bool=false)
     print("─┼─")
     print("────")
     print("─┼─")
-    print("────────")
+    print("──────────────")
     print("─┼─")
     print("────────")
     print("─┼─")
@@ -179,7 +179,7 @@ end
         success::Bool,
         time::Real,
         obj::Real,
-        ref_obj::Real,
+        ref_obj::Union{Real, Nothing},
         iterations::Union{Int, Nothing} = nothing,
         memory_bytes::Union{Int, Nothing} = nothing,
         show_memory::Bool = false
@@ -225,13 +225,13 @@ function print_test_line(
     success::Bool,
     time::Real,
     obj::Real,
-    ref_obj::Real,
+    ref_obj::Union{Real,Nothing},
     iterations::Union{Int,Nothing}=nothing,
     memory_bytes::Union{Int,Nothing}=nothing,
     show_memory::Bool=false,
 )
     # Relative error calculation
-    rel_error = abs(obj - ref_obj) / abs(ref_obj) * 100
+    rel_error = ref_obj === nothing ? missing : abs(obj - ref_obj) / abs(ref_obj) * 100
 
     # Colored status (✓ green or ✗ red)
     if success
@@ -248,7 +248,7 @@ function print_test_line(
 
     # Fixed columns with rpad/lpad (like CTBenchmarks)
     # Problem: 8 characters
-    printstyled(rpad(problem, 8); color=:cyan, bold=true)
+    printstyled(rpad(problem, 14); color=:cyan, bold=true)
     print(" │ ")
 
     # Discretizer: 8 characters
@@ -282,16 +282,16 @@ function print_test_line(
     print(" │ ")
 
     # Reference: scientific format with phantom sign
-    ref_str = @sprintf("%.6e", ref_obj)
-    if ref_obj >= 0
+    ref_str = ref_obj === nothing ? "N/A" : @sprintf("%.6e", ref_obj)
+    if ref_obj !== nothing && ref_obj >= 0
         ref_str = " " * ref_str  # Phantom sign
     end
     print(lpad(ref_str, 14))
     print(" │ ")
 
     # Error: scientific notation with 2 decimal places
-    err_str = @sprintf("%.2e", rel_error / 100)  # Convert to fraction then scientific format
-    err_color = rel_error < 1 ? :green : (rel_error < 5 ? :yellow : :red)
+    err_str = rel_error === missing ? "N/A" : @sprintf("%.2e", rel_error / 100)  # Convert to fraction then scientific format
+    err_color = rel_error === missing ? :white : (rel_error < 1 ? :green : (rel_error < 5 ? :yellow : :red))
     printstyled(lpad(err_str, 7); color=err_color)
 
     # Memory: optional, right-aligned, 10 characters

test/problems/TestProblems.jl

@@ -5,8 +5,11 @@ using OptimalControl
 include("beam.jl")
 include("goddard.jl")
 include("double_integrator.jl")
+include("quadrotor.jl")
+include("transfer.jl")
 
 export Beam, Goddard
 export DoubleIntegratorTime, DoubleIntegratorEnergy, DoubleIntegratorEnergyConstrained
+export Quadrotor, Transfer
 
 end

test/problems/beam.jl

@@ -22,7 +22,10 @@ function Beam()
         ∫(u(t)^2) → min
     end
 
-    init = (state=[0.05, 0.1], control=0.1)
+    init = @init ocp begin
+        x(t) := [0.05, 0.1]
+        u(t) := 0.1
+    end
 
     return (ocp=ocp, obj=8.898598, name="beam", init=init)
 end

test/problems/double_integrator.jl

@@ -56,6 +56,7 @@ function DoubleIntegratorTime()
         ocp=ocp,
         obj=2.0,
         name="double_integrator_time",
+        init=nothing,
         x0=x0,
         xf=xf,
         t0=t0,
@@ -107,7 +108,7 @@ function DoubleIntegratorEnergy()
         0.5∫(u(t)^2) → min
     end
 
-    return (ocp=ocp, obj=6.0, name="double_integrator_energy", x0=x0, xf=xf, t0=t0, tf=tf)
+    return (ocp=ocp, obj=6.0, name="double_integrator_energy", init=nothing, x0=x0, xf=xf, t0=t0, tf=tf)
 end
 
 """
@@ -159,8 +160,9 @@ function DoubleIntegratorEnergyConstrained()
 
     return (
         ocp=ocp,
-        obj=nothing,
+        obj=7.680030,
         name="double_integrator_energy_constrained",
+        init=nothing,
         x0=x0,
         xf=xf,
         t0=t0,

test/problems/goddard.jl

@@ -58,8 +58,11 @@ function Goddard(; vmax=0.1, Tmax=3.5)
     end
 
     # Components for a reasonable initial guess around a feasible trajectory.
-    init = (state=[1.01, 0.05, 0.8], control=0.5, variable=[0.1])
-
+    init = @init goddard begin
+        x(t) := [1.01, 0.05, 0.8]
+        u(t) := 0.5
+        tf := 0.1
+    end
     # Dynamics functions for indirect methods
     function F0(x)
         r, v, m = x

test/problems/quadrotor.jl

@@ -0,0 +1,54 @@
+# Quadrotor optimal control problem definition used by tests and examples.
+#
+# Returns a NamedTuple with fields:
+#   - ocp  :: the CTParser-defined optimal control problem
+#   - obj  :: reference optimal objective value (Ipopt / MadNLP, Collocation)
+#   - name :: a short problem name
+#   - init :: NamedTuple of components for CTSolvers.initial_guess
+function Quadrotor(; T=1, g=9.8, r=0.1)
+
+    ocp = @def begin
+        t ∈ [0, T], time
+        x ∈ R⁹, state
+        u ∈ R⁴, control
+
+        x(0) == zeros(9)
+
+        ∂(x₁)(t) == x₂(t)
+        ∂(x₂)(t) ==
+        u₁(t) * cos(x₇(t)) * sin(x₈(t)) * cos(x₉(t)) + u₁(t) * sin(x₇(t)) * sin(x₉(t))
+        ∂(x₃)(t) == x₄(t)
+        ∂(x₄)(t) ==
+        u₁(t) * cos(x₇(t)) * sin(x₈(t)) * sin(x₉(t)) - u₁(t) * sin(x₇(t)) * cos(x₉(t))
+        ∂(x₅)(t) == x₆(t)
+        ∂(x₆)(t) == u₁(t) * cos(x₇(t)) * cos(x₈(t)) - g
+        ∂(x₇)(t) == u₂(t) * cos(x₇(t)) / cos(x₈(t)) + u₃(t) * sin(x₇(t)) / cos(x₈(t))
+        ∂(x₈)(t) == -u₂(t) * sin(x₇(t)) + u₃(t) * cos(x₇(t))
+        ∂(x₉)(t) ==
+        u₂(t) * cos(x₇(t)) * tan(x₈(t)) + u₃(t) * sin(x₇(t)) * tan(x₈(t)) + u₄(t)
+
+        dt1 = sin(2π * t / T)
+        df1 = 0
+        dt3 = 2sin(4π * t / T)
+        df3 = 0
+        dt5 = 2t / T
+        df5 = 2
+
+        0.5∫(
+            (x₁(t) - dt1)^2 +
+            (x₃(t) - dt3)^2 +
+            (x₅(t) - dt5)^2 +
+            x₇(t)^2 +
+            x₈(t)^2 +
+            x₉(t)^2 +
+            r * (u₁(t)^2 + u₂(t)^2 + u₃(t)^2 + u₄(t)^2),
+        ) → min
+    end
+
+    init = @init ocp begin
+        x(t) := 0.1 * ones(9) 
+        u(t) := 0.1 * ones(4) 
+    end
+
+    return (ocp=ocp, obj=4.2679623758, name="quadrotor", init=init)
+end

test/problems/transfer.jl

@@ -0,0 +1,92 @@
+# Transfer optimal control problem definition used by tests and examples.
+#
+# Returns a NamedTuple with fields:
+#   - ocp  :: the CTParser-defined optimal control problem
+#   - obj  :: reference optimal objective value (Ipopt / MadNLP, Collocation)
+#   - name :: a short problem name
+#   - init :: NamedTuple of components for CTSolvers.initial_guess
+
+asqrt(x; ε=1e-9) = sqrt(sqrt(x^2+ε^2))     # Avoid issues with AD
+
+const μ = 5165.8620912                     # Earth gravitation constant
+
+function F0(x)
+    P, ex, ey, hx, hy, L = x
+    pdm = asqrt(P / μ)
+    cl = cos(L)
+    sl = sin(L)
+    w = 1 + ex * cl + ey * sl
+    F = [0, 0, 0, 0, 0, w^2 / (P * pdm)]
+    return F
+end
+
+function F1(x)
+    P, ex, ey, hx, hy, L = x
+    pdm = asqrt(P / μ)
+    cl = cos(L)
+    sl = sin(L)
+    F = pdm * [0, sl, -cl, 0, 0, 0]
+    return F
+end
+
+function F2(x)
+    P, ex, ey, hx, hy, L = x
+    pdm = asqrt(P / μ)
+    cl = cos(L)
+    sl = sin(L)
+    w = 1 + ex * cl + ey * sl
+    F = pdm * [2 * P / w, cl + (ex + cl) / w, sl + (ey + sl) / w, 0, 0, 0]
+    return F
+end
+
+function F3(x)
+    P, ex, ey, hx, hy, L = x
+    pdm = asqrt(P / μ)
+    cl = cos(L)
+    sl = sin(L)
+    w = 1 + ex * cl + ey * sl
+    pdmw = pdm / w
+    zz = hx * sl - hy * cl
+    uh = (1 + hx^2 + hy^2) / 2
+    F = pdmw * [0, -zz * ey, zz * ex, uh * cl, uh * sl, zz]
+    return F
+end
+
+function Transfer(; Tmax=60)
+
+    cTmax = 3600^2 / 1e6; T = Tmax * cTmax     # Conversion from Newtons to kg x Mm / h²
+    mass0 = 1500                               # Initial mass of the spacecraft
+    β = 1.42e-02                               # Engine specific impulsion
+    P0 = 11.625                                # Initial semilatus rectum
+    ex0, ey0 = 0.75, 0                         # Initial eccentricity
+    hx0, hy0 = 6.12e-2, 0                      # Initial ascending node and inclination
+    L0 = π                                     # Initial longitude
+    Pf = 42.165                                # Final semilatus rectum
+    exf, eyf = 0, 0                            # Final eccentricity
+    hxf, hyf = 0, 0                            # Final ascending node and inclination
+    Lf = 3π                                    # Estimation of final longitude
+    x0 = [P0, ex0, ey0, hx0, hy0, L0]          # Initial state
+    xf = [Pf, exf, eyf, hxf, hyf, Lf]          # Final state
+
+    ocp = @def begin
+        tf ∈ R, variable
+        t ∈ [0, tf], time
+        x = (P, ex, ey, hx, hy, L) ∈ R⁶, state
+        u ∈ R³, control
+        x(0) == x0
+        x[1:5](tf) == xf[1:5]
+        mass = mass0 - β * T * t
+        ẋ(t) == F0(x(t)) + T / mass * (u₁(t) * F1(x(t)) + u₂(t) * F2(x(t)) + u₃(t) * F3(x(t)))
+        u₁(t)^2 + u₂(t)^2 + u₃(t)^2 ≤ 1
+        tf → min
+    end
+
+    init = @init ocp begin 
+        tf_i = 15
+        x(t) := x0 + (xf - x0) * t / tf_i       # Linear interpolation
+        u(t) := [0.1, 0.5, 0.]                  # Initial guess for the control
+        tf := tf_i                              # Initial guess for final time
+    end
+
+    return (ocp=ocp, obj=14.79643132, name="transfer", init=init)
+end