Skip to content

Commit 862251d

Browse files
ocotsocots
committed
remove CTProblems dependence
1 parent 64e0904 commit 862251d

File tree

8 files changed

+210
-133
lines changed

8 files changed

+210
-133
lines changed

‎test/Goddard.jl‎

Lines changed: 70 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,70 @@
1+
function Goddard()
2+
3+
# parameters
4+
Cd = 310
5+
Tmax = 3.5
6+
β = 500
7+
b = 2
8+
t0 = 0
9+
r0 = 1
10+
v0 = 0
11+
vmax = 0.1
12+
m0 = 1
13+
mf = 0.6
14+
x0 = [ r0, v0, m0 ]
15+
16+
# the model
17+
@def ocp begin
18+
# parameters
19+
Cd = 310
20+
Tmax = 3.5
21+
β = 500
22+
b = 2
23+
t0 = 0
24+
r0 = 1
25+
v0 = 0
26+
vmax = 0.1
27+
m0 = 1
28+
mf = 0.6
29+
x0 = [ r0, v0, m0 ]
30+
31+
# variables
32+
tf ∈ R, variable
33+
t ∈ [ t0, tf ], time
34+
x ∈ R³, state
35+
u ∈ R, control
36+
r = x₁
37+
v = x₂
38+
m = x₃
39+
40+
# constraints
41+
0 ≤ u(t) ≤ 1, (u_con)
42+
r(t) ≥ r0, (x_con_rmin)
43+
0 ≤ v(t) ≤ vmax, (x_con_vmax)
44+
x(t0) == x0, (initial_con)
45+
m(tf) == mf, (final_con)
46+
47+
# dynamics
48+
ẋ(t) == F0(x(t)) + u(t)*F1(x(t))
49+
50+
# objective
51+
r(tf) → max
52+
end
53+
54+
# dynamics
55+
function F0(x)
56+
r, v, m = x
57+
D = Cd * v^2 * exp(-β*(r - 1))
58+
return [ v, -D/m - 1/r^2, 0 ]
59+
end
60+
function F1(x)
61+
r, v, m = x
62+
return [ 0, Tmax/m, -b*Tmax ]
63+
end
64+
65+
#
66+
objective = 1.0125762700748699
67+
68+
return ocp, objective
69+
70+
end

‎test/Project.toml‎

Lines changed: 1 addition & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1,9 +1,8 @@
11
[deps]
2-
CTProblems = "45d9ea3f-a92f-411f-833f-222dd4fb9cd8"
3-
DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
42
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
53
NLPModelsIpopt = "f4238b75-b362-5c4c-b852-0801c9a21d71"
64
NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
5+
OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
76
SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
87
Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
98

‎test/runtests.jl‎

Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1,12 +1,13 @@
11
using OptimalControl
22
using Test
33
using LinearAlgebra
4-
using CTProblems
54
using SciMLBase
65
using NonlinearSolve
7-
using DifferentialEquations
6+
using OrdinaryDiffEq
87
using NLPModelsIpopt
98

9+
include("Goddard.jl")
10+
1011
#
1112
@testset verbose = true showtiming = true "Optimal control tests" begin
1213
for name ∈ (

‎test/test_basic.jl‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -10,7 +10,7 @@ function test_basic()
1010
∫( 0.5u(t)^2 ) → min
1111
end
1212

13-
sol = solve(ocp)
13+
sol = solve(ocp; display=false)
1414
@test sol.objective ≈ 6 atol=1e-2
1515

1616
end

‎test/test_goddard_direct.jl‎

Lines changed: 3 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -1,12 +1,11 @@
11
function test_goddard_direct()
22

33
# goddard with state constraint - maximize altitude
4-
prob = Problem(:goddard, :classical, :altitude, :x_dim_3, :u_dim_1, :mayer, :x_cons, :u_cons, :singular_arc)
5-
ocp = prob.model
4+
ocp, objective = Goddard()
65

76
# initial guess (constant state and control functions)
87
init = (state=[1.01, 0.05, 0.8], control=0.1, variable=0.2)
9-
sol = solve(ocp, grid_size=10, print_level=0, init=init)
10-
@test sol.objective ≈ prob.solution.objective atol=5e-3
8+
sol = solve(ocp; grid_size=10, display=false, init=init)
9+
@test sol.objective ≈ objective atol=5e-3
1110

1211
end

‎test/test_goddard_indirect.jl‎

Lines changed: 111 additions & 112 deletions
Original file line numberDiff line numberDiff line change
@@ -1,116 +1,115 @@
11
function test_goddard_indirect()
22

3-
prob = Problem(:goddard, :classical, :altitude, :x_dim_3, :u_dim_1, :mayer, :x_cons, :u_cons, :singular_arc)
4-
ocp = prob.model
5-
sol = prob.solution
6-
title = prob.title
7-
8-
# parameters
9-
Cd = 310
10-
Tmax = 3.5
11-
β = 500
12-
b = 2
13-
t0 = 0
14-
r0 = 1
15-
v0 = 0
16-
vmax = 0.1
17-
m0 = 1
18-
mf = 0.6
19-
x0 = [ r0, v0, m0 ]
20-
21-
#
22-
remove_constraint!(ocp, :x_con_rmin)
23-
g(x) = vmax-constraint(ocp, :x_con_vmax)(x, Real[]) # g(x, u) ≥ 0 (cf. nonnegative multiplier)
24-
final_mass_cons(xf) = constraint(ocp, :final_con)(x0, xf, Real[])-mf
25-
26-
function F0(x)
27-
r, v, m = x
28-
D = Cd * v^2 * exp(-β*(r - 1))
29-
return [ v, -D/m - 1/r^2, 0 ]
30-
end
31-
32-
function F1(x)
33-
r, v, m = x
34-
return [ 0, Tmax/m, -b*Tmax ]
35-
end
36-
37-
# bang controls
38-
u0 = 0
39-
u1 = 1
40-
41-
# singular control
42-
H0 = Lift(F0)
43-
H1 = Lift(F1)
44-
H01 = @Lie {H0, H1}
45-
H001 = @Lie {H0, H01}
46-
H101 = @Lie {H1, H01}
47-
us(x, p) = -H001(x, p) / H101(x, p)
48-
49-
# boundary control
50-
ub(x) = -Lie(F0, g)(x) / Lie(F1, g)(x)
51-
μ(x, p) = H01(x, p) / Lie(F1, g)(x)
52-
53-
# flows
54-
f0 = Flow(ocp, (x, p, v) -> u0)
55-
f1 = Flow(ocp, (x, p, v) -> u1)
56-
fs = Flow(ocp, (x, p, v) -> us(x, p))
57-
fb = Flow(ocp, (x, p, v) -> ub(x), (x, u, v) -> g(x), (x, p, v) -> μ(x, p))
58-
59-
# shooting function
60-
function shoot!(s, p0, t1, t2, t3, tf) # B+ S C B0 structure
61-
62-
x1, p1 = f1(t0, x0, p0, t1)
63-
x2, p2 = fs(t1, x1, p1, t2)
64-
x3, p3 = fb(t2, x2, p2, t3)
65-
xf, pf = f0(t3, x3, p3, tf)
66-
s[1] = final_mass_cons(xf)
67-
s[2:3] = pf[1:2] - [ 1, 0 ]
68-
s[4] = H1(x1, p1)
69-
s[5] = H01(x1, p1)
70-
s[6] = g(x2)
71-
s[7] = H0(xf, pf) # free tf
72-
73-
end
74-
75-
# tests
76-
t1 = 0.025246759388000528
77-
t2 = 0.061602092906721286
78-
t3 = 0.10401664867856217
79-
tf = 0.20298394547952422
80-
p0 = [3.9428857983400074, 0.14628855388160236, 0.05412448008321635]
81-
82-
# test shooting function
83-
s = zeros(eltype(p0), 7)
84-
shoot!(s, p0, t1, t2, t3, tf)
85-
s_guess_sol = [-0.02456074767656735,
86-
-0.05699760226157302,
87-
0.0018629693253921868,
88-
-0.027013078908634858,
89-
-0.21558816838342798,
90-
-0.0121146739026253,
91-
0.015713236406057297]
92-
@test s ≈ s_guess_sol atol=1e-6
93-
94-
#
95-
ξ0 = [ p0 ; t1 ; t2 ; t3 ; tf ]
96-
97-
#
98-
foo!(s, ξ, λ) = shoot!(s, ξ[1:3], ξ[4], ξ[5], ξ[6], ξ[7])
99-
#sol = fsolve(foo!, ξ0, show_trace=true); println(sol)
100-
prob = NonlinearProblem(foo!, ξ0)
101-
#sol = NonlinearSolve.solve(prob)
102-
sol = solve(prob)
103-
104-
p0 = sol.u[1:3]
105-
t1 = sol.u[4]
106-
t2 = sol.u[5]
107-
t3 = sol.u[6]
108-
tf = sol.u[7];
109-
110-
shoot!(s, p0, t1, t2, t3, tf)
111-
112-
#@test sol.converged
113-
@test SciMLBase.successful_retcode(sol)
114-
@test norm(s) < 1e-6
3+
#
4+
include("Goddard.jl")
5+
ocp, obj = Goddard()
6+
7+
# parameters
8+
Cd = 310
9+
Tmax = 3.5
10+
β = 500
11+
b = 2
12+
t0 = 0
13+
r0 = 1
14+
v0 = 0
15+
vmax = 0.1
16+
m0 = 1
17+
mf = 0.6
18+
x0 = [ r0, v0, m0 ]
19+
20+
#
21+
remove_constraint!(ocp, :x_con_rmin)
22+
g(x) = vmax-constraint(ocp, :x_con_vmax)(x, Real[]) # g(x, u) ≥ 0 (cf. nonnegative multiplier)
23+
final_mass_cons(xf) = constraint(ocp, :final_con)(x0, xf, Real[])-mf
24+
25+
function F0(x)
26+
r, v, m = x
27+
D = Cd * v^2 * exp(-β*(r - 1))
28+
return [ v, -D/m - 1/r^2, 0 ]
29+
end
30+
31+
function F1(x)
32+
r, v, m = x
33+
return [ 0, Tmax/m, -b*Tmax ]
34+
end
35+
36+
# bang controls
37+
u0 = 0
38+
u1 = 1
39+
40+
# singular control
41+
H0 = Lift(F0)
42+
H1 = Lift(F1)
43+
H01 = @Lie {H0, H1}
44+
H001 = @Lie {H0, H01}
45+
H101 = @Lie {H1, H01}
46+
us(x, p) = -H001(x, p) / H101(x, p)
47+
48+
# boundary control
49+
ub(x) = -Lie(F0, g)(x) / Lie(F1, g)(x)
50+
μ(x, p) = H01(x, p) / Lie(F1, g)(x)
51+
52+
# flows
53+
f0 = Flow(ocp, (x, p, v) -> u0)
54+
f1 = Flow(ocp, (x, p, v) -> u1)
55+
fs = Flow(ocp, (x, p, v) -> us(x, p))
56+
fb = Flow(ocp, (x, p, v) -> ub(x), (x, u, v) -> g(x), (x, p, v) -> μ(x, p))
57+
58+
# shooting function
59+
function shoot!(s, p0, t1, t2, t3, tf) # B+ S C B0 structure
60+
61+
x1, p1 = f1(t0, x0, p0, t1)
62+
x2, p2 = fs(t1, x1, p1, t2)
63+
x3, p3 = fb(t2, x2, p2, t3)
64+
xf, pf = f0(t3, x3, p3, tf)
65+
s[1] = final_mass_cons(xf)
66+
s[2:3] = pf[1:2] - [ 1, 0 ]
67+
s[4] = H1(x1, p1)
68+
s[5] = H01(x1, p1)
69+
s[6] = g(x2)
70+
s[7] = H0(xf, pf) # free tf
71+
72+
end
73+
74+
# tests
75+
t1 = 0.025246759388000528
76+
t2 = 0.061602092906721286
77+
t3 = 0.10401664867856217
78+
tf = 0.20298394547952422
79+
p0 = [3.9428857983400074, 0.14628855388160236, 0.05412448008321635]
80+
81+
# test shooting function
82+
s = zeros(eltype(p0), 7)
83+
shoot!(s, p0, t1, t2, t3, tf)
84+
s_guess_sol = [-0.02456074767656735,
85+
-0.05699760226157302,
86+
0.0018629693253921868,
87+
-0.027013078908634858,
88+
-0.21558816838342798,
89+
-0.0121146739026253,
90+
0.015713236406057297]
91+
@test s ≈ s_guess_sol atol=1e-6
92+
93+
#
94+
ξ0 = [ p0 ; t1 ; t2 ; t3 ; tf ]
95+
96+
#
97+
foo!(s, ξ, λ) = shoot!(s, ξ[1:3], ξ[4], ξ[5], ξ[6], ξ[7])
98+
#sol = fsolve(foo!, ξ0, show_trace=true); println(sol)
99+
prob = NonlinearProblem(foo!, ξ0)
100+
#sol = NonlinearSolve.solve(prob)
101+
sol = solve(prob)
102+
103+
p0 = sol.u[1:3]
104+
t1 = sol.u[4]
105+
t2 = sol.u[5]
106+
t3 = sol.u[6]
107+
tf = sol.u[7];
108+
109+
shoot!(s, p0, t1, t2, t3, tf)
110+
111+
#@test sol.converged
112+
@test SciMLBase.successful_retcode(sol)
113+
@test norm(s) < 1e-6
115114

116115
end

‎test/test_grid.jl‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -13,7 +13,7 @@ function test_grid()
1313
objective!(ocp, :lagrange, (x, u) -> (u[1]+u[2])^2)
1414

1515
time_grid = [0,0.1,0.3,0.6,0.98,0.99,1]
16-
sol = solve(ocp, time_grid=time_grid, print_level=0)
16+
sol = solve(ocp; time_grid=time_grid, display=false)
1717
@test sol.objective ≈ 0.309 rtol=1e-2
1818

1919
end

0 commit comments

Comments
 (0)