717 dev add some integration tests for the new solve

test/problems/TestProblems.jl

@@ -4,7 +4,9 @@ module TestProblems
 
     include("beam.jl")
     include("goddard.jl")
+    include("quadrotor.jl")
+    include("transfer.jl")
 
-    export Beam, Goddard
+    export Beam, Goddard, Quadrotor, Transfer
 
-end
+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/goddard.jl

@@ -58,7 +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
 
     return (ocp=goddard, obj=1.01257, name="goddard", init=init)
 end

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

test/suite/solve/test_canonical.jl

@@ -9,7 +9,7 @@
 module TestCanonical
 
 import Test
-import OptimalControl
+import OptimalControl: OptimalControl, GPU
 
 # Import du module d'affichage (DIP - dépend de l'abstraction)
 include(joinpath(@__DIR__, "..", "..", "helpers", "print_utils.jl"))
@@ -68,8 +68,10 @@ function test_canonical()
         ]
 
         problems = [
+            #("Transfer", Transfer()), # debug: add later (currently issue with :exa)
             ("Beam",    Beam()),
             ("Goddard", Goddard()),
+            ("Quadrotor", Quadrotor()),
         ]
 
         # ----------------------------------------------------------------
@@ -138,8 +140,8 @@ function test_canonical()
         # GPU tests (only if CUDA is available)
         # ----------------------------------------------------------------
         if is_cuda_on()
-            gpu_modeler  = ("Exa/GPU", OptimalControl.Exa(backend=CUDA.CUDABackend()))
-            gpu_solver   = ("MadNLP/GPU",    OptimalControl.MadNLP(print_level=MadNLP.ERROR, linear_solver=MadNLPGPU.CUDSSSolver))
+            gpu_modeler  = ("Exa/GPU", OptimalControl.Exa{GPU}())
+            gpu_solver   = ("MadNLP/GPU", OptimalControl.MadNLP{GPU}(print_level=MadNLP.ERROR))
 
             for (pname, pb) in problems
                 Test.@testset "GPU / $pname" begin