diff --git a/test/problems/TestProblems.jl b/test/problems/TestProblems.jl index 2ec09eea..694bcb91 100644 --- a/test/problems/TestProblems.jl +++ b/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 \ No newline at end of file +end diff --git a/test/problems/beam.jl b/test/problems/beam.jl index 542d7500..5927cbee 100644 --- a/test/problems/beam.jl +++ b/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 diff --git a/test/problems/goddard.jl b/test/problems/goddard.jl index 310adcdb..03588a30 100644 --- a/test/problems/goddard.jl +++ b/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 diff --git a/test/problems/quadrotor.jl b/test/problems/quadrotor.jl new file mode 100644 index 00000000..b9379844 --- /dev/null +++ b/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 diff --git a/test/problems/transfer.jl b/test/problems/transfer.jl new file mode 100644 index 00000000..072b7339 --- /dev/null +++ b/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 diff --git a/test/suite/solve/test_canonical.jl b/test/suite/solve/test_canonical.jl index b883a52f..a5a2c602 100644 --- a/test/suite/solve/test_canonical.jl +++ b/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