Skip to content

Commit 0d3fa6e

Browse files
ocotsocots
committed
quick beam
1 parent a43cd71 commit 0d3fa6e

File tree

5 files changed

+166
-70
lines changed

5 files changed

+166
-70
lines changed

‎test/figures/init/beam.pdf‎

48.6 KB
Binary file not shown.

‎test/runtests.jl‎

Lines changed: 17 additions & 14 deletions
Original file line numberDiff line numberDiff line change
@@ -32,14 +32,17 @@ for file in files
3232
end
3333

3434
# Remove from the tests the following problems
35-
# problems_to_exclude = [
36-
# :quadrotor
37-
# ]
38-
# list_of_problems = setdiff(list_of_problems, problems_to_exclude)
35+
problems_to_exclude = [
36+
:bioreactor,
37+
:cart_pendulum,
38+
:dielectrophoretic_particle,
39+
:moon_lander,
40+
]
41+
list_of_problems = setdiff(list_of_problems, problems_to_exclude)
3942

40-
# list_of_problems = [
41-
# :truck_trailer
42-
# ]
43+
list_of_problems = [
44+
:beam
45+
]
4346

4447
# The list of all the problems to test
4548
const LIST_OF_PROBLEMS = deepcopy(list_of_problems)
@@ -53,13 +56,13 @@ const VERBOSE = true # print or not details during tests
5356
@testset "OptimalControlProblems tests" verbose=VERBOSE showtiming=true begin
5457
for name in (
5558
#:aqua,
56-
#:kwargs,
57-
#:JuMP, # convergence tests for JuMP models
58-
#:OptimalControl, # convergence tests for OptimalControl models
59-
:OptimalControl_s, # convergence tests for OptimalControl models
60-
#:init, # comparison between OptimalControl and JuMP: init
61-
#:solution, # comparison between OptimalControl and JuMP: solution
62-
#:quick, # quick comparison: objective rel error only
59+
# :kwargs,
60+
# :JuMP, # convergence tests for JuMP models
61+
# :OptimalControl, # convergence tests for OptimalControl models
62+
# :OptimalControl_s, # convergence tests for OptimalControl models
63+
:init, # comparison between OptimalControl and JuMP: init
64+
# :solution, # comparison between OptimalControl and JuMP: solution
65+
# :quick, # quick comparison: objective rel error only
6366
)
6467
@testset "$(name)" verbose=VERBOSE begin
6568
test_name = Symbol(:test_, name)

‎test/test_OptimalControl_s.jl‎

Lines changed: 8 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -1,12 +1,13 @@
11
# test_OptimalControl_optimality
22
function test_OptimalControl_s()
3-
kwargs = Dict(
4-
:print_level => 0,
3+
kwargs_madnlp = Dict(
4+
:print_level => ERROR,
55
:tol => TOL,
6-
:mu_strategy => MU_STRATEGY,
7-
:sb => SB,
6+
#:mu_strategy => MU_STRATEGY,
7+
#:sb => SB,
88
:max_iter => MAX_ITER,
99
:max_wall_time => MAX_WALL_TIME,
10+
linear_solver => MumpsSolver,
1011
)
1112

1213
for f in LIST_OF_PROBLEMS
@@ -17,7 +18,7 @@ function test_OptimalControl_s()
1718
keep_problem = true
1819

1920
#
20-
DEBUG && println("\n", "┌─ ", string(f), " (OptimalControl)")
21+
DEBUG && println("\n", "┌─ ", string(f), " (OptimalControl_s)")
2122
DEBUG && println("│")
2223

2324
# Set up the model
@@ -28,9 +29,9 @@ function test_OptimalControl_s()
2829
DEBUG && println("├─ Solve")
2930
DEBUG && println("│")
3031
print(" First solve: ");
31-
@time sol = madnlp(nlp; linear_solver=MumpsSolver) # debug , kwargs...)
32+
@time sol = madnlp(nlp; kwargs...)
3233
print(" Second solve: ");
33-
@time sol = madnlp(nlp; linear_solver=MumpsSolver) # debug , kwargs...)
34+
@time sol = madnlp(nlp; kwargs...)
3435
DEBUG && println("│")
3536

3637
# Infos

‎test/test_quick.jl‎

Lines changed: 63 additions & 19 deletions
Original file line numberDiff line numberDiff line change
@@ -7,7 +7,7 @@ function test_quick()
77
ε_abs_objective = 1e-6
88

99
# options for solvers
10-
kwargs = Dict(
10+
kwargs_ipopt = Dict(
1111
:print_level => 0,
1212
:tol => TOL,
1313
:mu_strategy => MU_STRATEGY,
@@ -16,6 +16,16 @@ function test_quick()
1616
:max_wall_time => MAX_WALL_TIME,
1717
)
1818

19+
kwargs_madnlp = Dict(
20+
:print_level => MadNLP.ERROR,
21+
:tol => TOL,
22+
#:mu_strategy => MU_STRATEGY,
23+
#:sb => SB,
24+
:max_iter => MAX_ITER,
25+
:max_wall_time => MAX_WALL_TIME,
26+
:linear_solver => MumpsSolver,
27+
)
28+
1929
max_r_err = -Inf # relative error max
2030

2131
for f in LIST_OF_PROBLEMS
@@ -25,13 +35,6 @@ function test_quick()
2535
DEBUG && println("\n", "┌─ ", string(f))
2636
DEBUG && println("│")
2737

28-
########## OptimalControl ##########
29-
docp = OptimalControlProblems.eval(f)(OptimalControlBackend(); N=N)
30-
nlp = nlp_model(docp)
31-
nlp_sol = NLPModelsIpopt.ipopt(nlp; kwargs...)
32-
sol = build_ocp_solution(docp, nlp_sol)
33-
o_oc = objective(sol)
34-
3538
############### JuMP ###############
3639
nlp = OptimalControlProblems.eval(f)(JuMPBackend(); N=N)
3740
set_optimizer(nlp, Ipopt.Optimizer)
@@ -45,26 +48,67 @@ function test_quick()
4548
optimize!(nlp)
4649
o_jp = objective_value(nlp)
4750

48-
############### TEST ###############
49-
# objective
50-
o_di = abs(o_oc-o_jp)
51-
o_bd = max(0.5*(abs(o_oc) + abs(o_jp))*ε_rel_objective, ε_abs_objective)
51+
########## OptimalControl ##########
52+
docp = OptimalControlProblems.eval(f)(OptimalControlBackend(); N=N)
53+
nlp = nlp_model(docp)
54+
nlp_sol = NLPModelsIpopt.ipopt(nlp; kwargs_ipopt...)
55+
sol = build_ocp_solution(docp, nlp_sol)
56+
o_oc = objective(sol)
5257

58+
########## OptimalControl_s ##########
59+
docp = OptimalControlProblems.eval(Symbol(f, :_s))(OptimalControlBackend(), :madnlp, :exa; N=N)
60+
nlp = nlp_model(docp)
61+
nlp_sol = madnlp(nlp; kwargs_madnlp...)
62+
sol = build_ocp_solution(docp, nlp_sol)
63+
o_os = objective(sol)
64+
65+
############### TEST ###############
5366
DEBUG && println("├─ objective")
5467
DEBUG && println("│")
55-
DEBUG && println("│ o_oc = ", o_oc)
5668
DEBUG && println("│ o_jp = ", o_jp)
57-
DEBUG && println("│ r_err = ", o_di/(0.5*(abs(o_oc) + abs(o_jp))))
58-
DEBUG && println("│ a_err = ", o_di)
59-
DEBUG && println("│ bound = ", o_bd)
69+
DEBUG && println("│ o_oc = ", o_oc)
70+
DEBUG && println("│ o_os = ", o_os)
71+
DEBUG && println("│")
72+
73+
# comparison JuMP and OptimalControl
74+
A = o_oc
75+
B = o_jp
76+
77+
o_di = abs(A-B)
78+
o_bd = max(0.5*(abs(A) + abs(B))*ε_rel_objective, ε_abs_objective)
79+
max_r_err = max(max_r_err, o_di/(0.5*(abs(A) + abs(B))))
80+
81+
DEBUG && println("│ JuMP vs OptimalControl")
82+
DEBUG && println("│")
83+
DEBUG && println("│ r_err = ", o_di/(0.5*(abs(A) + abs(B))))
84+
DEBUG && println("│ a_err = ", o_di)
85+
DEBUG && println("│ bound = ", o_bd)
6086

6187
res = @my_test_broken o_di < o_bd
6288

63-
DEBUG && res && println("│ \033[1;32mTest Passed\033[0m")
64-
DEBUG && !res && println("│ \033[1;31mTest Failed\033[0m")
89+
DEBUG && res && println("│ \033[1;32mTest Passed\033[0m")
90+
DEBUG && !res && println("│ \033[1;31mTest Failed\033[0m")
91+
DEBUG && println("│")
92+
93+
# comparison JuMP and OptimalControl_s
94+
A = o_os
95+
B = o_jp
96+
97+
o_di = abs(A-B)
98+
o_bd = max(0.5*(abs(A) + abs(B))*ε_rel_objective, ε_abs_objective)
99+
max_r_err = max(max_r_err, o_di/(0.5*(abs(A) + abs(B))))
100+
101+
DEBUG && println("│ JuMP vs OptimalControl_s")
65102
DEBUG && println("│")
103+
DEBUG && println("│ r_err = ", o_di/(0.5*(abs(A) + abs(B))))
104+
DEBUG && println("│ a_err = ", o_di)
105+
DEBUG && println("│ bound = ", o_bd)
106+
107+
res = @my_test_broken o_di < o_bd
66108

67-
max_r_err = max(max_r_err, o_di/(0.5*(abs(o_oc) + abs(o_jp))))
109+
DEBUG && res && println("│ \033[1;32mTest Passed\033[0m")
110+
DEBUG && !res && println("│ \033[1;31mTest Failed\033[0m")
111+
DEBUG && println("│")
68112

69113
#
70114
DEBUG && println("└─")

‎test/utils.jl‎

Lines changed: 78 additions & 30 deletions
Original file line numberDiff line numberDiff line change
@@ -97,8 +97,8 @@ function comparison(; max_iter, test_name)
9797
ε_rel_control = 1e-1
9898
ε_abs_control = 1e-6
9999

100-
# options for solvers
101-
Options = Dict(
100+
# options_ipopt for solvers
101+
options_ipopt = Dict(
102102
:print_level => 0,
103103
:tol => TOL,
104104
:mu_strategy => MU_STRATEGY,
@@ -107,6 +107,16 @@ function comparison(; max_iter, test_name)
107107
:max_wall_time => MAX_WALL_TIME,
108108
)
109109

110+
options_madnlp = Dict(
111+
:print_level => MadNLP.ERROR,
112+
:tol => TOL,
113+
#:mu_strategy => MU_STRATEGY,
114+
#:sb => SB,
115+
:max_iter => max_iter,
116+
:max_wall_time => MAX_WALL_TIME,
117+
:linear_solver => MumpsSolver,
118+
)
119+
110120
# we loop over the problems
111121
for f in LIST_OF_PROBLEMS
112122
N = metadata[f][:N] # get default N
@@ -119,30 +129,16 @@ function comparison(; max_iter, test_name)
119129
DEBUG && println("\n┌─ ", string(f), " (", string(test_name), ")")
120130
DEBUG && println("│")
121131

122-
########## OptimalControl ##########
123-
docp = OptimalControlProblems.eval(f)(OptimalControlBackend(); N=N)
124-
nlp_oc = nlp_model(docp)
125-
nlp_sol = NLPModelsIpopt.ipopt(nlp_oc; Options...)
126-
sol = build_ocp_solution(docp, nlp_sol)
127-
sol_oc = deepcopy(sol) # for plotting
128-
129-
t_oc = time_grid(sol)
130-
x_oc = state(sol).(t_oc)
131-
u_oc = control(sol).(t_oc)
132-
o_oc = objective(sol)
133-
i_oc = iterations(sol)
134-
v_oc = variable(sol)
135-
136132
############### JuMP ###############
137133
nlp_jp = OptimalControlProblems.eval(f)(JuMPBackend(); N=N)
138134
set_optimizer(nlp_jp, Ipopt.Optimizer)
139135
set_silent(nlp_jp)
140-
set_optimizer_attribute(nlp_jp, "tol", Options[:tol])
141-
set_optimizer_attribute(nlp_jp, "max_iter", Options[:max_iter])
142-
set_optimizer_attribute(nlp_jp, "mu_strategy", Options[:mu_strategy])
136+
set_optimizer_attribute(nlp_jp, "tol", options_ipopt[:tol])
137+
set_optimizer_attribute(nlp_jp, "max_iter", options_ipopt[:max_iter])
138+
set_optimizer_attribute(nlp_jp, "mu_strategy", options_ipopt[:mu_strategy])
143139
set_optimizer_attribute(nlp_jp, "linear_solver", "mumps")
144-
set_optimizer_attribute(nlp_jp, "max_wall_time", Options[:max_wall_time])
145-
set_optimizer_attribute(nlp_jp, "sb", Options[:sb])
140+
set_optimizer_attribute(nlp_jp, "max_wall_time", options_ipopt[:max_wall_time])
141+
set_optimizer_attribute(nlp_jp, "sb", options_ipopt[:sb])
146142
optimize!(nlp_jp)
147143

148144
t_jp = time_grid(f, nlp_jp)
@@ -153,8 +149,34 @@ function comparison(; max_iter, test_name)
153149
v_jp = variable(f, nlp_jp)
154150
p_jp = costate(f, nlp_jp).(t_jp)
155151

152+
########## OptimalControl ##########
153+
docp = OptimalControlProblems.eval(f)(OptimalControlBackend(); N=N)
154+
nlp_oc = nlp_model(docp)
155+
nlp_sol = NLPModelsIpopt.ipopt(nlp_oc; options_ipopt...)
156+
sol_oc = build_ocp_solution(docp, nlp_sol)
157+
158+
t_oc = time_grid(sol_oc)
159+
x_oc = state(sol_oc).(t_oc)
160+
u_oc = control(sol_oc).(t_oc)
161+
o_oc = objective(sol_oc)
162+
i_oc = iterations(sol_oc)
163+
v_oc = variable(sol_oc)
164+
165+
########## OptimalControl_s ##########
166+
docp = OptimalControlProblems.eval(Symbol(f, :_s))(OptimalControlBackend(), :madnlp, :exa; N=N)
167+
nlp_os = nlp_model(docp)
168+
nlp_sol = madnlp(nlp_os; options_madnlp...)
169+
sol_os = build_ocp_solution(docp, nlp_sol)
170+
171+
t_os = time_grid(sol_os)
172+
x_os = state(sol_os).(t_os)
173+
u_os = control(sol_os).(t_os)
174+
o_os = objective(sol_os)
175+
i_os = iterations(sol_os)
176+
v_os = variable(sol_os)
177+
156178
########## Iterations ##########
157-
DEBUG && println("├─ Iterations → OC: ", i_oc, ", JP: ", i_jp)
179+
DEBUG && println("├─ Iterations → JP: ", i_jp, ", OC: ", i_oc, ", OS: ", i_os)
158180

159181
keep_problem = true
160182

@@ -303,18 +325,43 @@ function comparison(; max_iter, test_name)
303325
@assert(length(p_vars)==n)
304326

305327
# OptimalControl
328+
color = 1
306329
labelOC = if (test_name == :solution)
307330
"OptimalControl: " * string(i_oc) * " it"
308331
else
309332
"OptimalControl"
310333
end
311334
plt = plot(
312335
sol_oc;
313-
state_style=(color=1,),
314-
costate_style=(color=1, legend=:none),
315-
control_style=(color=1, legend=:none),
316-
path_style=(color=1, legend=:none),
317-
dual_style=(color=1, legend=:none),
336+
state_style=(color=color,),
337+
costate_style=(color=color, legend=:none),
338+
control_style=(color=color, legend=:none),
339+
path_style=(color=color, legend=:none),
340+
dual_style=(color=color, legend=:none),
341+
size=(900, 220*(n+m)),
342+
label=labelOC,
343+
leftmargin=20mm,
344+
)
345+
for i in 2:n
346+
plot!(plt[i]; legend=:none)
347+
end
348+
349+
# OptimalControl_s
350+
color = 2
351+
labelOC = if (test_name == :solution)
352+
"OptimalControl_s: " * string(i_oc) * " it"
353+
else
354+
"OptimalControl_s"
355+
end
356+
plot!(
357+
plt,
358+
sol_os;
359+
linestyle=:dot,
360+
state_style=(color=color,),
361+
costate_style=(color=color, legend=:none),
362+
control_style=(color=color, legend=:none),
363+
path_style=(color=color, legend=:none),
364+
dual_style=(color=color, legend=:none),
318365
size=(900, 220*(n+m)),
319366
label=labelOC,
320367
leftmargin=20mm,
@@ -324,21 +371,22 @@ function comparison(; max_iter, test_name)
324371
end
325372

326373
# JuMP
374+
color = 3
327375
labelJP = (test_name == :solution) ? "JuMP: " * string(i_jp) * " it" : "JuMP"
328376
for i in eachindex(x_vars) # state
329377
xi_jp = [x_jp[k][i] for k in eachindex(t_jp)]
330378
label = i == 1 ? labelJP : :none
331-
plot!(plt[i], t_jp, xi_jp; color=2, linestyle=:dash, label=label)
379+
plot!(plt[i], t_jp, xi_jp; color=color, linestyle=:dash, label=label)
332380
end
333381

334382
for i in eachindex(p_vars) # costate
335383
pi_jp = [p_jp[k][i] for k in eachindex(t_jp)]
336-
plot!(plt[n + i], t_jp, -pi_jp; color=2, linestyle=:dash, label=:none)
384+
plot!(plt[n + i], t_jp, -pi_jp; color=color, linestyle=:dash, label=:none)
337385
end
338386

339387
for i in eachindex(u_vars) # control
340388
ui_jp = [u_jp[k][i] for k in eachindex(t_jp)]
341-
plot!(plt[2n + i], t_jp, ui_jp; color=2, linestyle=:dash, label=:none)
389+
plot!(plt[2n + i], t_jp, ui_jp; color=color, linestyle=:dash, label=:none)
342390
end
343391

344392
# save figure

0 commit comments

Comments
 (0)