Skip to content

Commit 0073ea1

Browse files
authored
Merge pull request #737 from control-toolbox/add-tests
Refactor canonical tests with separated modelers/solvers
2 parents 30c046a + 28a2aa2 commit 0073ea1

File tree

8 files changed

+310
-202
lines changed

8 files changed

+310
-202
lines changed

test/helpers/print_utils.jl

Lines changed: 10 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -110,7 +110,7 @@ function print_test_header(show_memory::Bool=false)
110110
print("")
111111
printstyled(rpad("Type", 4); bold=true)
112112
print("")
113-
printstyled(rpad("Problem", 8); bold=true)
113+
printstyled(rpad("Problem", 14); bold=true)
114114
print("")
115115
printstyled(rpad("Disc", 8); bold=true)
116116
print("")
@@ -140,7 +140,7 @@ function print_test_header(show_memory::Bool=false)
140140
print("─┼─")
141141
print("────")
142142
print("─┼─")
143-
print("────────")
143+
print("──────────────")
144144
print("─┼─")
145145
print("────────")
146146
print("─┼─")
@@ -179,7 +179,7 @@ end
179179
success::Bool,
180180
time::Real,
181181
obj::Real,
182-
ref_obj::Real,
182+
ref_obj::Union{Real, Nothing},
183183
iterations::Union{Int, Nothing} = nothing,
184184
memory_bytes::Union{Int, Nothing} = nothing,
185185
show_memory::Bool = false
@@ -225,13 +225,13 @@ function print_test_line(
225225
success::Bool,
226226
time::Real,
227227
obj::Real,
228-
ref_obj::Real,
228+
ref_obj::Union{Real,Nothing},
229229
iterations::Union{Int,Nothing}=nothing,
230230
memory_bytes::Union{Int,Nothing}=nothing,
231231
show_memory::Bool=false,
232232
)
233233
# Relative error calculation
234-
rel_error = abs(obj - ref_obj) / abs(ref_obj) * 100
234+
rel_error = ref_obj === nothing ? missing : abs(obj - ref_obj) / abs(ref_obj) * 100
235235

236236
# Colored status (✓ green or ✗ red)
237237
if success
@@ -248,7 +248,7 @@ function print_test_line(
248248

249249
# Fixed columns with rpad/lpad (like CTBenchmarks)
250250
# Problem: 8 characters
251-
printstyled(rpad(problem, 8); color=:cyan, bold=true)
251+
printstyled(rpad(problem, 14); color=:cyan, bold=true)
252252
print("")
253253

254254
# Discretizer: 8 characters
@@ -282,16 +282,16 @@ function print_test_line(
282282
print("")
283283

284284
# Reference: scientific format with phantom sign
285-
ref_str = @sprintf("%.6e", ref_obj)
286-
if ref_obj >= 0
285+
ref_str = ref_obj === nothing ? "N/A" : @sprintf("%.6e", ref_obj)
286+
if ref_obj !== nothing && ref_obj >= 0
287287
ref_str = " " * ref_str # Phantom sign
288288
end
289289
print(lpad(ref_str, 14))
290290
print("")
291291

292292
# Error: scientific notation with 2 decimal places
293-
err_str = @sprintf("%.2e", rel_error / 100) # Convert to fraction then scientific format
294-
err_color = rel_error < 1 ? :green : (rel_error < 5 ? :yellow : :red)
293+
err_str = rel_error === missing ? "N/A" : @sprintf("%.2e", rel_error / 100) # Convert to fraction then scientific format
294+
err_color = rel_error === missing ? :white : (rel_error < 1 ? :green : (rel_error < 5 ? :yellow : :red))
295295
printstyled(lpad(err_str, 7); color=err_color)
296296

297297
# Memory: optional, right-aligned, 10 characters

test/problems/TestProblems.jl

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -5,8 +5,11 @@ using OptimalControl
55
include("beam.jl")
66
include("goddard.jl")
77
include("double_integrator.jl")
8+
include("quadrotor.jl")
9+
include("transfer.jl")
810

911
export Beam, Goddard
1012
export DoubleIntegratorTime, DoubleIntegratorEnergy, DoubleIntegratorEnergyConstrained
13+
export Quadrotor, Transfer
1114

1215
end

test/problems/beam.jl

Lines changed: 4 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -22,7 +22,10 @@ function Beam()
2222
(u(t)^2) min
2323
end
2424

25-
init = (state=[0.05, 0.1], control=0.1)
25+
init = @init ocp begin
26+
x(t) := [0.05, 0.1]
27+
u(t) := 0.1
28+
end
2629

2730
return (ocp=ocp, obj=8.898598, name="beam", init=init)
2831
end

test/problems/double_integrator.jl

Lines changed: 4 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -56,6 +56,7 @@ function DoubleIntegratorTime()
5656
ocp=ocp,
5757
obj=2.0,
5858
name="double_integrator_time",
59+
init=nothing,
5960
x0=x0,
6061
xf=xf,
6162
t0=t0,
@@ -107,7 +108,7 @@ function DoubleIntegratorEnergy()
107108
0.5(u(t)^2) min
108109
end
109110

110-
return (ocp=ocp, obj=6.0, name="double_integrator_energy", x0=x0, xf=xf, t0=t0, tf=tf)
111+
return (ocp=ocp, obj=6.0, name="double_integrator_energy", init=nothing, x0=x0, xf=xf, t0=t0, tf=tf)
111112
end
112113

113114
"""
@@ -159,8 +160,9 @@ function DoubleIntegratorEnergyConstrained()
159160

160161
return (
161162
ocp=ocp,
162-
obj=nothing,
163+
obj=7.680030,
163164
name="double_integrator_energy_constrained",
165+
init=nothing,
164166
x0=x0,
165167
xf=xf,
166168
t0=t0,

test/problems/goddard.jl

Lines changed: 5 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -58,8 +58,11 @@ function Goddard(; vmax=0.1, Tmax=3.5)
5858
end
5959

6060
# Components for a reasonable initial guess around a feasible trajectory.
61-
init = (state=[1.01, 0.05, 0.8], control=0.5, variable=[0.1])
62-
61+
init = @init goddard begin
62+
x(t) := [1.01, 0.05, 0.8]
63+
u(t) := 0.5
64+
tf := 0.1
65+
end
6366
# Dynamics functions for indirect methods
6467
function F0(x)
6568
r, v, m = x

test/problems/quadrotor.jl

Lines changed: 54 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,54 @@
1+
# Quadrotor optimal control problem definition used by tests and examples.
2+
#
3+
# Returns a NamedTuple with fields:
4+
# - ocp :: the CTParser-defined optimal control problem
5+
# - obj :: reference optimal objective value (Ipopt / MadNLP, Collocation)
6+
# - name :: a short problem name
7+
# - init :: NamedTuple of components for CTSolvers.initial_guess
8+
function Quadrotor(; T=1, g=9.8, r=0.1)
9+
10+
ocp = @def begin
11+
t [0, T], time
12+
x R⁹, state
13+
u R⁴, control
14+
15+
x(0) == zeros(9)
16+
17+
(x₁)(t) == x₂(t)
18+
(x₂)(t) ==
19+
u₁(t) * cos(x₇(t)) * sin(x₈(t)) * cos(x₉(t)) + u₁(t) * sin(x₇(t)) * sin(x₉(t))
20+
(x₃)(t) == x₄(t)
21+
(x₄)(t) ==
22+
u₁(t) * cos(x₇(t)) * sin(x₈(t)) * sin(x₉(t)) - u₁(t) * sin(x₇(t)) * cos(x₉(t))
23+
(x₅)(t) == x₆(t)
24+
(x₆)(t) == u₁(t) * cos(x₇(t)) * cos(x₈(t)) - g
25+
(x₇)(t) == u₂(t) * cos(x₇(t)) / cos(x₈(t)) + u₃(t) * sin(x₇(t)) / cos(x₈(t))
26+
(x₈)(t) == -u₂(t) * sin(x₇(t)) + u₃(t) * cos(x₇(t))
27+
(x₉)(t) ==
28+
u₂(t) * cos(x₇(t)) * tan(x₈(t)) + u₃(t) * sin(x₇(t)) * tan(x₈(t)) + u₄(t)
29+
30+
dt1 = sin(2π * t / T)
31+
df1 = 0
32+
dt3 = 2sin(4π * t / T)
33+
df3 = 0
34+
dt5 = 2t / T
35+
df5 = 2
36+
37+
0.5(
38+
(x₁(t) - dt1)^2 +
39+
(x₃(t) - dt3)^2 +
40+
(x₅(t) - dt5)^2 +
41+
x₇(t)^2 +
42+
x₈(t)^2 +
43+
x₉(t)^2 +
44+
r * (u₁(t)^2 + u₂(t)^2 + u₃(t)^2 + u₄(t)^2),
45+
) min
46+
end
47+
48+
init = @init ocp begin
49+
x(t) := 0.1 * ones(9)
50+
u(t) := 0.1 * ones(4)
51+
end
52+
53+
return (ocp=ocp, obj=4.2679623758, name="quadrotor", init=init)
54+
end

test/problems/transfer.jl

Lines changed: 92 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,92 @@
1+
# Transfer optimal control problem definition used by tests and examples.
2+
#
3+
# Returns a NamedTuple with fields:
4+
# - ocp :: the CTParser-defined optimal control problem
5+
# - obj :: reference optimal objective value (Ipopt / MadNLP, Collocation)
6+
# - name :: a short problem name
7+
# - init :: NamedTuple of components for CTSolvers.initial_guess
8+
9+
asqrt(x; ε=1e-9) = sqrt(sqrt(x^2+ε^2)) # Avoid issues with AD
10+
11+
const μ = 5165.8620912 # Earth gravitation constant
12+
13+
function F0(x)
14+
P, ex, ey, hx, hy, L = x
15+
pdm = asqrt(P / μ)
16+
cl = cos(L)
17+
sl = sin(L)
18+
w = 1 + ex * cl + ey * sl
19+
F = [0, 0, 0, 0, 0, w^2 / (P * pdm)]
20+
return F
21+
end
22+
23+
function F1(x)
24+
P, ex, ey, hx, hy, L = x
25+
pdm = asqrt(P / μ)
26+
cl = cos(L)
27+
sl = sin(L)
28+
F = pdm * [0, sl, -cl, 0, 0, 0]
29+
return F
30+
end
31+
32+
function F2(x)
33+
P, ex, ey, hx, hy, L = x
34+
pdm = asqrt(P / μ)
35+
cl = cos(L)
36+
sl = sin(L)
37+
w = 1 + ex * cl + ey * sl
38+
F = pdm * [2 * P / w, cl + (ex + cl) / w, sl + (ey + sl) / w, 0, 0, 0]
39+
return F
40+
end
41+
42+
function F3(x)
43+
P, ex, ey, hx, hy, L = x
44+
pdm = asqrt(P / μ)
45+
cl = cos(L)
46+
sl = sin(L)
47+
w = 1 + ex * cl + ey * sl
48+
pdmw = pdm / w
49+
zz = hx * sl - hy * cl
50+
uh = (1 + hx^2 + hy^2) / 2
51+
F = pdmw * [0, -zz * ey, zz * ex, uh * cl, uh * sl, zz]
52+
return F
53+
end
54+
55+
function Transfer(; Tmax=60)
56+
57+
cTmax = 3600^2 / 1e6; T = Tmax * cTmax # Conversion from Newtons to kg x Mm / h²
58+
mass0 = 1500 # Initial mass of the spacecraft
59+
β = 1.42e-02 # Engine specific impulsion
60+
P0 = 11.625 # Initial semilatus rectum
61+
ex0, ey0 = 0.75, 0 # Initial eccentricity
62+
hx0, hy0 = 6.12e-2, 0 # Initial ascending node and inclination
63+
L0 = π # Initial longitude
64+
Pf = 42.165 # Final semilatus rectum
65+
exf, eyf = 0, 0 # Final eccentricity
66+
hxf, hyf = 0, 0 # Final ascending node and inclination
67+
Lf = 3π # Estimation of final longitude
68+
x0 = [P0, ex0, ey0, hx0, hy0, L0] # Initial state
69+
xf = [Pf, exf, eyf, hxf, hyf, Lf] # Final state
70+
71+
ocp = @def begin
72+
tf R, variable
73+
t [0, tf], time
74+
x = (P, ex, ey, hx, hy, L) R⁶, state
75+
u R³, control
76+
x(0) == x0
77+
x[1:5](tf) == xf[1:5]
78+
mass = mass0 - β * T * t
79+
(t) == F0(x(t)) + T / mass * (u₁(t) * F1(x(t)) + u₂(t) * F2(x(t)) + u₃(t) * F3(x(t)))
80+
u₁(t)^2 + u₂(t)^2 + u₃(t)^2 1
81+
tf min
82+
end
83+
84+
init = @init ocp begin
85+
tf_i = 15
86+
x(t) := x0 + (xf - x0) * t / tf_i # Linear interpolation
87+
u(t) := [0.1, 0.5, 0.] # Initial guess for the control
88+
tf := tf_i # Initial guess for final time
89+
end
90+
91+
return (ocp=ocp, obj=14.79643132, name="transfer", init=init)
92+
end

0 commit comments

Comments
 (0)