From eb58261dd99dc2245c384df0f7fbfcb49aa4a7f3 Mon Sep 17 00:00:00 2001 From: Pierre Martinon Date: Mon, 22 Jul 2024 17:10:51 +0200 Subject: [PATCH] unified test files (#182) * unified: abstract_ocp test * unified: constraints test * unified: objective test * unified: misc test * unified: initial guess test * unified: continuation test * unified test for grid --- test/{test_grid.jl => refine_grid.jl} | 83 +------ test/suite/continuation.jl | 82 ------- test/suite/grid.jl | 55 ----- .../{abstract_ocp.jl => test_abstract_ocp.jl} | 0 .../{constraints.jl => test_constraints.jl} | 0 test/suite/test_continuation.jl | 107 ++++++++ test/suite/test_grid.jl | 79 ++++++ ...initial_guess.jl => test_initial_guess.jl} | 0 test/suite/{misc.jl => test_misc.jl} | 2 +- .../suite/{objective.jl => test_objective.jl} | 0 test/test_abstract_ocp.jl | 45 ---- test/test_constraints.jl | 29 --- test/test_continuation.jl | 90 ------- test/test_initial_guess.jl | 231 ------------------ test/test_misc.jl | 64 ----- test/test_objective.jl | 53 ---- 16 files changed, 189 insertions(+), 731 deletions(-) rename test/{test_grid.jl => refine_grid.jl} (63%) delete mode 100644 test/suite/continuation.jl delete mode 100644 test/suite/grid.jl rename test/suite/{abstract_ocp.jl => test_abstract_ocp.jl} (100%) rename test/suite/{constraints.jl => test_constraints.jl} (100%) create mode 100644 test/suite/test_continuation.jl create mode 100644 test/suite/test_grid.jl rename test/suite/{initial_guess.jl => test_initial_guess.jl} (100%) rename test/suite/{misc.jl => test_misc.jl} (97%) rename test/suite/{objective.jl => test_objective.jl} (100%) delete mode 100644 test/test_abstract_ocp.jl delete mode 100644 test/test_constraints.jl delete mode 100644 test/test_continuation.jl delete mode 100644 test/test_initial_guess.jl delete mode 100644 test/test_misc.jl delete mode 100644 test/test_objective.jl diff --git a/test/test_grid.jl b/test/refine_grid.jl similarity index 63% rename from test/test_grid.jl rename to test/refine_grid.jl index a9b134d..2113d75 100644 --- a/test/test_grid.jl +++ b/test/refine_grid.jl @@ -2,89 +2,10 @@ include("deps.jl") using Plots using Printf -println("Test: grid options") -test1 = false -test2 = false -test3 = false +println("Test: refined grid options") test4 = false # use time_grid=:refined to do this internally, ie in solve call a continuation_steps(...) ? test5 = false # time_grid=:optimized ? -test6 = true # time_grid=:optimized ? - -# 1. simple integrator min energy (dual control for test) -if test1 - ocp = Model() - state!(ocp, 1) - control!(ocp, 2) - time!(ocp, t0=0, tf=1) - constraint!(ocp, :initial, lb=-1, ub=-1) - constraint!(ocp, :final, lb=0, ub=0) - constraint!(ocp, :control, lb=[0,0], ub=[Inf, Inf]) - dynamics!(ocp, (x, u) -> -x - u[1] + u[2]) - objective!(ocp, :lagrange, (x, u) -> (u[1]+u[2])^2) - sol0 = solve(ocp, print_level=0) - - # solve with explicit and non uniform time grid - time_grid = LinRange(0,1,CTDirect.__grid_size_direct()+1) - sol5 = solve(ocp, time_grid=time_grid, print_level=0) - println((sol5.objective == sol0.objective) && (sol5.iterations == sol0.iterations)) - - time_grid = [0,0.1,0.3,0.6,0.98,0.99,1] - sol6 = solve(ocp, time_grid=time_grid, print_level=0) - println(sol6.objective) -end - -# 2. integrator free times -if test2 - ocp = Model(variable=true) - state!(ocp, 2) - control!(ocp, 1) - variable!(ocp, 2) - time!(ocp, ind0=1, indf=2) - constraint!(ocp, :initial, lb=[0,0], ub=[0,0]) - constraint!(ocp, :final, lb=[1,0], ub=[1,0]) - constraint!(ocp, :control, lb=-1, ub=1) - constraint!(ocp, :variable, lb=[0.1,0.1], ub=[10,10]) - constraint!(ocp, :variable, f=v->v[2]-v[1], lb=0.1, ub=Inf) - dynamics!(ocp, (x, u, v) -> [x[2], u]) - objective!(ocp, :mayer, (x0, xf, v) -> v[1], :max) - - sol = solve(ocp, time_grid=LinRange(0,1,CTDirect.__grid_size_direct()+1), print_level=0, tol=1e-12) - println("Target 8.0, found ", sol.objective) - - sol = solve(ocp, time_grid=[0,0.1,0.3,0.5,0.6,0.8,0.95,1], print_level=0) - plot(sol, show=true) - println("Target 8.0, coarse grid ", sol.objective) -end - -# 3. parametric ocp -if test3 - function ocp_T(T) - @def ocp begin - t ∈ [ 0, T ], time - x ∈ R², state - u ∈ R, control - q = x₁ - v = x₂ - q(0) == 0 - v(0) == 0 - q(T) == 1 - v(T) == 0 - ẋ(t) == [ v(t), u(t) ] - ∫(u(t)^2) → min - end - return ocp - end - - ocpT2 = ocp_T(2) - solT2 = solve(ocpT2, print_level=0) - - solT2_exp = solve(ocpT2, time_grid=LinRange(0,1,CTDirect.__grid_size_direct()+1),print_level=0) - println("T=2 Check explicit grid ", (solT2.objective==solT2_exp.objective) && (solT2.iterations==solT2_exp.iterations)) - - solT2_nonunif = solve(ocpT2, time_grid=[0,0.3,1,1.9,2],print_level=0) - println("T=2 with non-uniform grid ", solT2_nonunif.objective) - plot(solT2_nonunif, show=true) -end +test6 = false # time_grid=:optimized ? # 4. pseudo grid refinement with manual grid input if test4 diff --git a/test/suite/continuation.jl b/test/suite/continuation.jl deleted file mode 100644 index ba6378a..0000000 --- a/test/suite/continuation.jl +++ /dev/null @@ -1,82 +0,0 @@ -println("Test: discrete continuation") - -# continuation on final time with parametric ocp definition -function ocp_T(T) - @def ocp begin - t ∈ [ 0, T ], time - x ∈ R², state - u ∈ R, control - q = x₁ - v = x₂ - q(0) == 0 - v(0) == 0 - q(T) == 1 - v(T) == 0 - ẋ(t) == [ v(t), u(t) ] - ∫(u(t)^2) → min - end - return ocp -end -@testset verbose = true showtiming = true ":parametric_ocp :warm_start" begin - init = () - obj_list = [] - for T=1:5 - ocp = ocp_T(T) - sol = solve(ocp, print_level=0, init=init) - init = sol - push!(obj_list, sol.objective) - end - @test obj_list ≈ [12, 1.5, 0.44, 0.19, 0.096] rtol=1e-2 -end - -#= -# continuation with parametric definition of the ocp -relu(x) = max(0, x) -μ = 10 -p_relu(x) = log(abs(1 + exp(μ*x)))/μ -f(x) = 1-x -m(x) = (p_relu∘f)(x) -T = 2 -function myocp(ρ) - @def ocp begin - τ ∈ R, variable - s ∈ [ 0, 1 ], time - x ∈ R², state - u ∈ R², control - x₁(0) == 0 - x₂(0) == 1 - x₁(1) == 1 - ẋ(s) == [τ*(u₁(s)+2), (T-τ)*u₂(s)] - -1 ≤ u₁(s) ≤ 1 - -1 ≤ u₂(s) ≤ 1 - 0 ≤ τ ≤ T - -(x₂(1)-2)^3 - ∫( ρ * ( τ*m(x₁(s))^2 + (T-τ)*m(x₂(s))^2 ) ) → min - end - return ocp -end -@testset verbose = true showtiming = true ":parametric_ocp :warm_start" begin - init = () - obj_list = [] - for ρ in [0.1, 5, 10, 30, 100] - ocp = myocp(ρ) - sol = solve(ocp, print_level=0, init=init) - init = sol - push!(obj_list, sol.objective) - end - @test obj_list ≈ [-0.034, -1.7, -6.2, -35, -148] rtol=1e-2 -end -=# - -# goddard -sol0 = solve(goddard, print_level=0) - -@testset verbose = true showtiming = true ":global_variable :warm_start" begin - sol = sol0 - obj_list = [] - for Tmax_local=3.5:-0.5:1 - global Tmax = Tmax_local - sol = solve(goddard, print_level=0, init=sol) - push!(obj_list, sol.objective) - end - @test obj_list ≈ [1.0125, 1.0124, 1.0120, 1.0112, 1.0092, 1.0036] rtol=1e-2 -end diff --git a/test/suite/grid.jl b/test/suite/grid.jl deleted file mode 100644 index c3a4b2c..0000000 --- a/test/suite/grid.jl +++ /dev/null @@ -1,55 +0,0 @@ -println("Test: grid options") - -# simple integrator min energy -# control split as positive/negative parts for m=2 tets case -ocp = Model() -state!(ocp, 1) -control!(ocp, 2) -time!(ocp, t0=0, tf=1) -constraint!(ocp, :initial, lb=-1, ub=-1) -constraint!(ocp, :final, lb=0, ub=0) -constraint!(ocp, :control, lb=[0,0], ub=[Inf, Inf]) -dynamics!(ocp, (x, u) -> -x - u[1] + u[2]) -objective!(ocp, :lagrange, (x, u) -> (u[1]+u[2])^2) -sol0 = solve(ocp, print_level=0) - -# solve with explicit and non uniform time grid -@testset verbose = true showtiming = true ":explicit_grid" begin - time_grid = LinRange(0,1,CTDirect.__grid_size_direct()+1) - sol5 = solve(ocp, time_grid=time_grid, print_level=0) - @test (sol5.objective == sol0.objective) && (sol5.iterations == sol0.iterations) -end - -@testset verbose = true showtiming = true ":non_uniform_grid" begin - time_grid = [0,0.1,0.3,0.6,0.98,0.99,1] - sol6 = solve(ocp, time_grid=time_grid, print_level=0) - @test sol6.objective ≈ 0.309 rtol=1e-2 -end - - -# recheck solution (T=2) with explicit / non-uniform grid -@def ocpT2 begin - t ∈ [ 0, 2 ], time - x ∈ R², state - u ∈ R, control - q = x₁ - v = x₂ - q(0) == 0 - v(0) == 0 - q(2) == 1 - v(2) == 0 - ẋ(t) == [ v(t), u(t) ] - ∫(u(t)^2) → min -end - -solT2 = solve(ocpT2, print_level=0) - -@testset verbose = true showtiming = true ":explicit_grid" begin - solT2_exp = solve(ocpT2, time_grid=LinRange(0,1,CTDirect.__grid_size_direct()+1),print_level=0) - @test (solT2_exp.objective == solT2.objective) && (solT2_exp.iterations == solT2.iterations) -end - -@testset verbose = true showtiming = true ":non_uniform_grid" begin - solT2_nonunif = solve(ocpT2, time_grid=[0,0.3,1,1.9,2],print_level=0) - @test solT2_nonunif.objective ≈ 2.43 rtol=1e-2 -end \ No newline at end of file diff --git a/test/suite/abstract_ocp.jl b/test/suite/test_abstract_ocp.jl similarity index 100% rename from test/suite/abstract_ocp.jl rename to test/suite/test_abstract_ocp.jl diff --git a/test/suite/constraints.jl b/test/suite/test_constraints.jl similarity index 100% rename from test/suite/constraints.jl rename to test/suite/test_constraints.jl diff --git a/test/suite/test_continuation.jl b/test/suite/test_continuation.jl new file mode 100644 index 0000000..3a66b95 --- /dev/null +++ b/test/suite/test_continuation.jl @@ -0,0 +1,107 @@ +println("Test: discrete continuation") + +test1 = true +test2 = true +test3 = true +draw_plot = false + +# parametric ocp definition +if test1 + function ocp_T(T) + @def ocp begin + t ∈ [ 0, T ], time + x ∈ R², state + u ∈ R, control + q = x₁ + v = x₂ + q(0) == 0 + v(0) == 0 + q(T) == 1 + v(T) == 0 + ẋ(t) == [ v(t), u(t) ] + ∫(u(t)^2) → min + end + return ocp + end + @testset verbose = true showtiming = true ":parametric_ocp :warm_start" begin + init = () + obj_list = [] + for T=1:5 + ocp = ocp_T(T) + sol = solve(ocp, print_level=0, init=init) + init = sol + push!(obj_list, sol.objective) + end + @test obj_list ≈ [12, 1.5, 0.44, 0.19, 0.096] rtol=1e-2 + end +end + + +# parametric ocp definition +if test2 + relu(x) = max(0, x) + μ = 10 + p_relu(x) = log(abs(1 + exp(μ*x)))/μ + f(x) = 1-x + m(x) = (p_relu∘f)(x) + T = 2 + function myocp(ρ) + @def ocp begin + τ ∈ R, variable + s ∈ [ 0, 1 ], time + x ∈ R², state + u ∈ R², control + x₁(0) == 0 + x₂(0) == 1 + x₁(1) == 1 + ẋ(s) == [τ*(u₁(s)+2), (T-τ)*u₂(s)] + -1 ≤ u₁(s) ≤ 1 + -1 ≤ u₂(s) ≤ 1 + 0 ≤ τ ≤ T + -(x₂(1)-2)^3 - ∫( ρ * ( τ*m(x₁(s))^2 + (T-τ)*m(x₂(s))^2 ) ) → min + end + return ocp + end + + @testset verbose = true showtiming = true ":parametric_ocp :warm_start" begin + init = () + obj_list = [] + for ρ in [0.1, 5, 10, 30, 100] + ocp = myocp(ρ) + sol = solve(ocp, print_level=0, init=init) + init = sol + push!(obj_list, sol.objective) + end + @test obj_list ≈ [-0.034, -1.7, -6.2, -35, -148] rtol=1e-2 + end +end + + +# global variable used in ocp +if test3 + Tmax = 3.5 + sol0 = solve(goddard, print_level=0) + + @testset verbose = true showtiming = true ":global_variable :warm_start" begin + sol = sol0 + Tmax_list = [] + obj_list = [] + for Tmax_local=3.5:-0.5:1 + global Tmax = Tmax_local + sol = solve(goddard, print_level=0, init=sol) + push!(Tmax_list, Tmax) + push!(obj_list, sol.objective) + end + @test obj_list ≈ [1.0125, 1.0124, 1.0120, 1.0112, 1.0092, 1.0036] rtol=1e-2 + + if draw_plot + using Plots + # plot obj(vmax) + pobj = plot(Tmax_list, obj_list, label="r(tf)", xlabel="Maximal thrust (Tmax)", ylabel="Maximal altitude r(tf)",seriestype=:scatter) + # plot multiple solutions + plot(sol0) + p = plot!(sol) + display(plot(pobj, p, layout=2, reuse=false, size=(1000,500))) + end + end +end diff --git a/test/suite/test_grid.jl b/test/suite/test_grid.jl new file mode 100644 index 0000000..af52220 --- /dev/null +++ b/test/suite/test_grid.jl @@ -0,0 +1,79 @@ +println("Test: grid options") + +# 1. simple integrator min energy (dual control for test) +ocp = Model() +state!(ocp, 1) +control!(ocp, 2) +time!(ocp, t0=0, tf=1) +constraint!(ocp, :initial, lb=-1, ub=-1) +constraint!(ocp, :final, lb=0, ub=0) +constraint!(ocp, :control, lb=[0,0], ub=[Inf, Inf]) +dynamics!(ocp, (x, u) -> -x - u[1] + u[2]) +objective!(ocp, :lagrange, (x, u) -> (u[1]+u[2])^2) +sol0 = solve(ocp, print_level=0) + +# solve with explicit and non uniform time grid +@testset verbose = true showtiming = true ":explicit_grid" begin + time_grid = LinRange(0,1,CTDirect.__grid_size_direct()+1) + sol = solve(ocp, time_grid=time_grid, print_level=0) + @test (sol.objective == sol0.objective) && (sol.iterations == sol0.iterations) +end + +@testset verbose = true showtiming = true ":non_uniform_grid" begin + time_grid = [0,0.1,0.3,0.6,0.98,0.99,1] + sol = solve(ocp, time_grid=time_grid, print_level=0) + @test sol.objective ≈ 0.309 rtol=1e-2 +end + + +# 2. integrator free times +ocp = Model(variable=true) +state!(ocp, 2) +control!(ocp, 1) +variable!(ocp, 2) +time!(ocp, ind0=1, indf=2) +constraint!(ocp, :initial, lb=[0,0], ub=[0,0]) +constraint!(ocp, :final, lb=[1,0], ub=[1,0]) +constraint!(ocp, :control, lb=-1, ub=1) +constraint!(ocp, :variable, lb=[0.1,0.1], ub=[10,10]) +constraint!(ocp, :variable, f=v->v[2]-v[1], lb=0.1, ub=Inf) +dynamics!(ocp, (x, u, v) -> [x[2], u]) +objective!(ocp, :mayer, (x0, xf, v) -> v[1], :max) +sol0 = solve(ocp, print_level=0) + +@testset verbose = true showtiming = true ":explicit_grid" begin + sol = solve(ocp, time_grid=LinRange(0,1,CTDirect.__grid_size_direct()+1), print_level=0) + @test (sol.objective == sol0.objective) && (sol.iterations == sol0.iterations) +end + +@testset verbose = true showtiming = true ":non_uniform_grid" begin + sol = solve(ocp, time_grid=[0,0.1,0.3,0.5,0.6,0.8,0.95,1], print_level=0) + #plot(sol, show=true) + @test sol.objective ≈ 7.96 rtol=1e-2 +end + +# 3. parametric ocp (T=2) with explicit / non-uniform grid +@def ocpT2 begin + t ∈ [ 0, 2 ], time + x ∈ R², state + u ∈ R, control + q = x₁ + v = x₂ + q(0) == 0 + v(0) == 0 + q(2) == 1 + v(2) == 0 + ẋ(t) == [ v(t), u(t) ] + ∫(u(t)^2) → min +end +sol0 = solve(ocpT2, print_level=0) + +@testset verbose = true showtiming = true ":explicit_grid" begin + sol = solve(ocpT2, time_grid=LinRange(0,1,CTDirect.__grid_size_direct()+1),print_level=0) + @test (sol.objective == sol0.objective) && (sol.iterations == sol0.iterations) +end + +@testset verbose = true showtiming = true ":non_uniform_grid" begin + sol = solve(ocpT2, time_grid=[0,0.3,1,1.9,2],print_level=0) + @test sol.objective ≈ 2.43 rtol=1e-2 +end diff --git a/test/suite/initial_guess.jl b/test/suite/test_initial_guess.jl similarity index 100% rename from test/suite/initial_guess.jl rename to test/suite/test_initial_guess.jl diff --git a/test/suite/misc.jl b/test/suite/test_misc.jl similarity index 97% rename from test/suite/misc.jl rename to test/suite/test_misc.jl index a5a55d9..e99c2f5 100644 --- a/test/suite/misc.jl +++ b/test/suite/test_misc.jl @@ -32,7 +32,7 @@ end @test sol.iterations == 5 end -# build solution from NLP (+++ add actual test ?) +# build solution from NLP (+++ add more tests ?) @testset verbose = true showtiming = true ":build_solution_from_nlp" begin docp, nlp = direct_transcription(ocp, init=sol0) dsol = CTDirect.solve_docp(docp, nlp, print_level=0, tol=1e-12) diff --git a/test/suite/objective.jl b/test/suite/test_objective.jl similarity index 100% rename from test/suite/objective.jl rename to test/suite/test_objective.jl diff --git a/test/test_abstract_ocp.jl b/test/test_abstract_ocp.jl deleted file mode 100644 index 04d4989..0000000 --- a/test/test_abstract_ocp.jl +++ /dev/null @@ -1,45 +0,0 @@ -include("deps.jl") -include("../problems/goddard.jl") - -println("Test: abstract OCP definition") - -# double integrator min tf, abstract definition -@def ocp1 begin - tf ∈ R, variable - t ∈ [ 0, tf ], time - x ∈ R², state - u ∈ R, control - -1 ≤ u(t) ≤ 1 - x(0) == [ 0, 0 ] - x(tf) == [ 1, 0 ] - 0.1 ≤ tf ≤ Inf - ẋ(t) == [ x₂(t), u(t) ] - tf → min -end - -sol1 = solve(ocp1, print_level=0, tol=1e-12) -println("Target 2.0, found ", sol1.objective) - - -# goddard -# NB. the ≤ is not the same as <= (parse error for <=) -@def ocp3 begin - tf ∈ R, variable - t ∈ [ 0, tf ], time - x ∈ R^3, state - u ∈ R, control - 0.1 ≤ tf ≤ Inf - r = x[1] - v = x[2] - m = x[3] - x(0) == [1, 0, 1] - m(tf) == 0.6 - 1 ≤ r(t) ≤ 1.1 - 0 ≤ v(t) ≤ vmax - 0 ≤ u(t) ≤ 1 - ẋ(t) == F0(x(t)) + u(t)*F1(x(t)) - r(tf) → max -end - -sol3 = solve(ocp3, print_level=0, tol=1e-12) -println("Target 1.0125, found ", sol3.objective) \ No newline at end of file diff --git a/test/test_constraints.jl b/test/test_constraints.jl deleted file mode 100644 index 1fad95e..0000000 --- a/test/test_constraints.jl +++ /dev/null @@ -1,29 +0,0 @@ -include("deps.jl") -include("../problems/goddard.jl") - -println("Test: constraint types") - - -# goddard, box constraints -ocp = goddard -sol = solve(ocp, print_level=0, tol=1e-8) -println("Target 1.0125, found ", sol.objective, " at ", sol.iterations, " iterations") - -# functional constraints -ocp1 = Model(variable=true) -state!(ocp1, 3) -control!(ocp1, 1) -variable!(ocp1, 1) -time!(ocp1, t0=0, indf=1) -constraint!(ocp1, :initial, lb=x0, ub=x0) -constraint!(ocp1, :final, rg=3, lb=mf, ub=Inf) -constraint!(ocp1, :state, f=(x,v)->x, lb=[r0,v0,mf], ub=[r0+0.2,vmax,m0]) -constraint!(ocp1, :control, f=(u,v)->u, lb=0, ub=1) -constraint!(ocp1, :variable, f=v->v, lb=0.01, ub=Inf) -objective!(ocp1, :mayer, (x0, xf, v) -> xf[1], :max) -dynamics!(ocp1, (x, u, v) -> F0(x) + u*F1(x) ) - -# note: the equations do not handle r<1 well -# without the box constraint on x, the default init (0.1) is not suitable -sol = solve(ocp1, print_level=0, tol=1e-8, init=(state=[1.01,0.05,0.75],)) -println("Target 1.0125, found ", sol.objective, " at ", sol.iterations, " iterations") diff --git a/test/test_continuation.jl b/test/test_continuation.jl deleted file mode 100644 index 0e9415a..0000000 --- a/test/test_continuation.jl +++ /dev/null @@ -1,90 +0,0 @@ -include("deps.jl") -using Printf -using Plots - -test1 = true -test2 = false -test3 = true - -# continuation on fixed final time -# NB. time! can be called only once, so we redefine the ocp -if test1 - println("\nDiscrete continuation on final time") - # parametric ocp - function ocp_T(T) - @def ocp begin - t ∈ [ 0, T ], time - x ∈ R², state - u ∈ R, control - q = x₁ - v = x₂ - q(0) == 0 - v(0) == 0 - q(T) == 1 - v(T) == 0 - ẋ(t) == [ v(t), u(t) ] - ∫(u(t)^2) → min - end - return ocp - end - - # continuation on final time - init1 = () - for T=1:5 - local ocp1 = ocp_T(T) - local sol1 = solve(ocp1, print_level=0, init=init1) - global init1 = sol1 - @printf("T %.2f objective %9.6f iterations %d\n", T, sol1.objective, sol1.iterations) - end - - - -end - -# continuation with parametric definition of the ocp -if test2 - println("\nDiscrete continuation on parametric ocp") - include("../problems/parametric.jl") #does not call function! - error("needs to be redone") - - # continuation on rho - init2 = () - ρs = [0.1, 5, 10, 30, 100] - for ρ in ρs - local ocp2 = myocp(ρ) - local sol2 = solve(ocp2, print_level=0, init=init2) - global init2 = sol2 - @printf("Rho %5.1f objective %9.4f iterations %d\n", ρ, sol2.objective, sol2.iterations) - end -end - - -# goddard max final altitude -if (test3) - include("../problems/goddard.jl") - ocp = goddard - # solve unconstrained problem - sol0 = solve(ocp, print_level=0) - #@printf("\nObjective for goddard reference solution %.6f", sol0.objective) - - # using a global variable in ocp definition - println("\nDiscrete continuation on maximal thrust") - # continuation on Tmax (using default init is slower) - Tmax_list = [] - obj_list = [] - sol3 = sol0 - for Tmax_local=3.5:-0.5:1 - global Tmax = Tmax_local - global sol3 = solve(ocp, print_level=0, init=sol3) - @printf("Tmax %.2f objective %.6f iterations %d\n", Tmax, sol3.objective, sol3.iterations) - push!(Tmax_list, Tmax) - push!(obj_list, sol3.objective) - end - - # plot obj(vmax) - pobj = plot(Tmax_list, obj_list, label="r(tf)", xlabel="Maximal thrust (Tmax)", ylabel="Maximal altitude r(tf)",seriestype=:scatter) - # plot multiple solutions - plot(sol0) - p = plot!(sol3) - display(plot(pobj, p, layout=2, reuse=false, size=(1000,500))) -end diff --git a/test/test_initial_guess.jl b/test/test_initial_guess.jl deleted file mode 100644 index 2cd2e90..0000000 --- a/test/test_initial_guess.jl +++ /dev/null @@ -1,231 +0,0 @@ -include("deps.jl") -using Printf -using Plots - -println("Test: initial guess options\n") - -# use 0 iterations to check initial guess, >0 to check cv -maxiter = 0 - -# test functions -function check_xf(sol, xf) - return xf == sol.state(sol.times[end]) -end -function check_uf(sol, uf) - return uf == sol.control(sol.times[end]) -end -function check_v(sol, v) - return v == sol.variable -end - -# reference solution -@def ocp begin - tf ∈ R, variable - t ∈ [ 0, tf ], time - x ∈ R², state - u ∈ R, control - -1 ≤ u(t) ≤ 1 - x(0) == [ 0, 0 ] - x(tf) == [ 1, 0 ] - 0.05 ≤ tf ≤ Inf - ẋ(t) == [ x₂(t), u(t) ] - tf → min -end -sol0 = solve(ocp, print_level=0) - -# constant initial guess -x_const = [0.5, 0.2] -u_const = 0.5 -v_const = 0.15 - -# functional initial guess -x_func = t->[t^2, sqrt(t)] -u_func = t->(cos(10*t)+1)*0.5 - -# interpolated initial gues -x_vec = [[0, 0], [1, 2], [5, -1]] -x_matrix = [0 0; 1 2; 5 -1] -u_vec = [0, 0.3, .1] - - -################################################# -# 1 Pass initial guess to all-in-one solve call -println("1. Passing the initial guess at the main solve level") - -# 1.a default initial guess -sol = solve(ocp, print_level=0, max_iter=maxiter) -if maxiter > 0 - @printf("%-56s %.3f at %d iterations\n", "Default initial guess", sol.objective, sol.iterations) -else - println(check_xf(sol, [0.1, 0.1]) && check_uf(sol, 0.1) && check_v(sol, 0.1)) -end -sol = solve(ocp, print_level=0, init=(), max_iter=maxiter) -if maxiter > 0 - @printf("%-56s %.3f at %d iterations\n", "Default initial guess", sol.objective, sol.iterations) -else - println(check_xf(sol, [0.1, 0.1]) && check_uf(sol, 0.1) && check_v(sol, 0.1)) -end -sol = solve(ocp, print_level=0, init=nothing, max_iter=maxiter) -if maxiter > 0 - @printf("%-56s %.3f at %d iterations\n", "Default initial guess", sol.objective, sol.iterations) -else - println(check_xf(sol, [0.1, 0.1]) && check_uf(sol, 0.1) && check_v(sol, 0.1)) -end - -# 1.b constant initial guess -sol = solve(ocp, print_level=0, init=(state=x_const,), max_iter=maxiter) -if maxiter > 0 - @printf("%-56s %.3f at %d iterations\n", "Constant x; default u,v", sol.objective, sol.iterations) -else - println(check_xf(sol, x_const)) -end - -sol = solve(ocp, print_level=0, init=(control=u_const,), max_iter=maxiter) -if maxiter > 0 - @printf("%-56s %.3f at %d iterations\n", "Constant u; default x,v", sol.objective, sol.iterations) -else - println(check_uf(sol, u_const)) -end - -sol = solve(ocp, print_level=0, init=(variable=v_const,), max_iter=maxiter) -if maxiter > 0 - @printf("%-56s %.3f at %d iterations\n", "Constant v; default x,u", sol.objective, sol.iterations) -else - println(check_v(sol, v_const)) -end - -sol = solve(ocp, print_level=0, init=(state=x_const, control=u_const), max_iter=maxiter) -if maxiter > 0 - @printf("%-56s %.3f at %d iterations\n", "Constant x,u; default v", sol.objective, sol.iterations) -else - println(check_xf(sol, x_const) && check_uf(sol, u_const)) -end - -sol = solve(ocp, print_level=0, init=(state=x_const, variable=v_const), max_iter=maxiter) -if maxiter > 0 - @printf("%-56s %.3f at %d iterations\n", "Constant x,v; default u", sol.objective, sol.iterations) -else - println(check_xf(sol, x_const) && check_v(sol, v_const)) -end - -sol = solve(ocp, print_level=0, init=(control=u_const, variable=v_const), max_iter=maxiter) -if maxiter > 0 - @printf("%-56s %.3f at %d iterations\n", "Constant u,v; default x", sol.objective, sol.iterations) -else - println(check_uf(sol, u_const) && check_v(sol, v_const)) -end - -sol = solve(ocp, print_level=0, init=(state=x_const, control=u_const, variable=v_const), max_iter=maxiter) -if maxiter > 0 - @printf("%-56s %.3f at %d iterations\n", "Constant x,u,v", sol.objective, sol.iterations) -else - println(check_xf(sol, x_const) && check_uf(sol, u_const) && check_v(sol, v_const)) -end - -# 1.c functional initial guess -sol = solve(ocp, print_level=0, init=(state=x_func,), max_iter=maxiter) -if maxiter > 0 - @printf("%-56s %.3f at %d iterations\n", "Functional x; default u,v", sol.objective, sol.iterations) -else - println(check_xf(sol, x_func(sol.times[end]))) -end - -sol = solve(ocp, print_level=0, init=(control=u_func,), max_iter=maxiter) -if maxiter > 0 - @printf("%-56s %.3f at %d iterations\n", "Functional u; default x,v", sol.objective, sol.iterations) -else - println(check_uf(sol, u_func(sol.times[end]))) -end - -sol = solve(ocp, print_level=0, init=(state=x_func, control=u_func), max_iter=maxiter) -if maxiter > 0 - @printf("%-56s %.3f at %d iterations\n", "Functional x,u; default v", sol.objective, sol.iterations) -else - println(check_xf(sol, x_func(sol.times[end])) && check_uf(sol, u_func(sol.times[end]))) -end - -# 1.d interpolated initial guess -t_vec = [0, .1, v_const] -sol = solve(ocp, print_level=0, init=(time=t_vec, state=x_vec, control=u_vec, variable=v_const), max_iter=maxiter) -if maxiter > 0 - @printf("%-56s %.3f at %d iterations\n", "Vector t,x,u; constant v", sol.objective, sol.iterations) -else - println(check_xf(sol, x_vec[end]) && check_uf(sol, u_vec[end]) && check_v(sol, v_const)) -end - -t_matrix = [0 .1 v_const] -sol = solve(ocp, print_level=0, init=(time=t_matrix, state=x_vec, control=u_vec, variable=v_const), max_iter=maxiter) -if maxiter > 0 - @printf("%-56s %.3f at %d iterations\n", "Matrix t; vector x,u; constant v", sol.objective, sol.iterations) -else - println(check_xf(sol, x_vec[end]) && check_uf(sol, u_vec[end]) && check_v(sol, v_const)) -end - -sol = solve(ocp, print_level=0, init=(time=t_vec, state=x_matrix, control=u_vec, variable=v_const), max_iter=maxiter) -if maxiter > 0 - @printf("%-56s %.3f at %d iterations\n", "Matrix x; vector t,u; constant v", sol.objective, sol.iterations) -else - println(check_xf(sol, x_vec[end]) && check_uf(sol, u_vec[end]) && check_v(sol, v_const)) -end - -# 1.e mixed initial guess -sol = solve(ocp, print_level=0, init=(time=t_vec, state=x_vec, control=u_func, variable=v_const), max_iter=maxiter) -if maxiter > 0 - @printf("%-56s %.3f at %d iterations\n", "Mixed: vector t,x; functional u; constant v", sol.objective, sol.iterations) -else - println(check_xf(sol, x_vec[end]) && check_uf(sol, u_func(sol.times[end])) && check_v(sol, v_const)) -end - -# 1.f warm start -sol = solve(ocp, print_level=0, init=sol0, max_iter=maxiter) -if maxiter > 0 - @printf("%-56s %.3f at %d iterations\n", "Warm start from reference solution", sol.objective, sol.iterations) -else - println(check_xf(sol, sol.state(sol.times[end])) && check_uf(sol, sol.control(sol.times[end])) && check_v(sol, sol.variable)) -end - -################################################# -# 2 Setting the initial guess at the DOCP level -println("\n2. Setting the initial guess at the DOCP level") -docp, nlp = direct_transcription(ocp) -# mixed init -set_initial_guess(docp, nlp, (time=t_vec, state=x_vec, control=u_func, variable=v_const)) -dsol = CTDirect.solve_docp(docp, nlp, print_level=0, max_iter=maxiter) -sol = build_solution(docp, dsol) -if maxiter > 0 - @printf("%-56s %.3f at %d iterations\n", "Mixed initial guess set in DOCP", sol.objective, sol.iterations) -else - println(check_xf(sol, x_vec[end]) && check_uf(sol, u_func(sol.times[end])) && check_v(sol, v_const)) -end - -# warm start -set_initial_guess(docp, nlp, sol0) -dsol = CTDirect.solve_docp(docp, nlp, print_level=0, max_iter=maxiter) -sol = build_solution(docp, dsol) -if maxiter > 0 - @printf("%-56s %.3f at %d iterations\n", "Warm start set in DOCP", sol.objective, sol.iterations) -else - println(check_xf(sol, sol.state(sol.times[end])) && check_uf(sol, sol.control(sol.times[end])) && check_v(sol, sol.variable)) -end - -################################################# -# 3 Passing the initial guess to solve call -println("\n3. Passing the initial guess to solve call") -set_initial_guess(docp, nlp, ()) # reset init in docp -# mixed init -dsol = CTDirect.solve_docp(docp, nlp, init=(time=t_vec, state=x_vec, control=u_func, variable=v_const), print_level=0, max_iter=maxiter) -sol = build_solution(docp, dsol) -if maxiter > 0 - @printf("%-56s %.3f at %d iterations\n", "Mixed initial guess passed to solve", sol.objective, sol.iterations) -else - println(check_xf(sol, x_vec[end]) && check_uf(sol, u_func(sol.times[end])) && check_v(sol, v_const)) -end - -# warm start -dsol = CTDirect.solve_docp(docp, nlp, init=sol0, print_level=0, max_iter=maxiter) -sol = build_solution(docp, dsol) -if maxiter > 0 - @printf("%-56s %.3f at %d iterations\n", "Warm start passed to solve", sol.objective, sol.iterations) -else - println(check_xf(sol, sol.state(sol.times[end])) && check_uf(sol, sol.control(sol.times[end])) && check_v(sol, sol.variable)) -end \ No newline at end of file diff --git a/test/test_misc.jl b/test/test_misc.jl deleted file mode 100644 index 683b025..0000000 --- a/test/test_misc.jl +++ /dev/null @@ -1,64 +0,0 @@ -include("deps.jl") -using Plots - -println("Test: misc") - -# simple integrator min energy -# control split as positive/negative parts for m=2 tets case -ocp = Model() -state!(ocp, 1) -control!(ocp, 2) -time!(ocp, t0=0, tf=1) -constraint!(ocp, :initial, lb=-1, ub=-1) -constraint!(ocp, :final, lb=0, ub=0) -constraint!(ocp, :control, lb=[0,0], ub=[Inf, Inf]) -dynamics!(ocp, (x, u) -> -x - u[1] + u[2]) -objective!(ocp, :lagrange, (x, u) -> (u[1]+u[2])^2) - -# all-in-one solve call -println(available_methods()) -println(is_solvable(ocp)) -println("Test simple integrator: all in one solve call") -sol = solve(ocp, print_level=0, tol=1e-12) -println("Target 0.313, found ", sol.objective, " at ", sol.iterations, " iterations") - -# split calls -println("Test simple integrator: split calls") -println("Direct transcription with default init") -docp, nlp = direct_transcription(ocp) -dsol = CTDirect.solve_docp(docp, nlp, print_level=0, tol=1e-12) -sol1 = build_solution(docp, dsol) -println("Target 0.313, found ", sol1.objective, " at ", sol1.iterations, " iterations") - -# test NLP getter -#nlp = get_nlp(docp) - -# warm start in directTranscription -println("\nDirect transcription with warm start") -docp2, nlp2 = direct_transcription(ocp, init=sol) -dsol2 = CTDirect.solve_docp(docp2, nlp2, print_level=0, tol=1e-12) - -# build_solution from NLP solution (primal only) -println("Rebuild OCP solution from NLP solution (primal only)") -sol3 = build_solution(docp2, primal=dsol2.solution) -plot(sol3, show=true) - -# build_solution from NLP solution (primal only) -println("Rebuild OCP solution from NLP solution (primal,dual)") -sol4 = build_solution(docp2, primal=dsol2.solution, dual=dsol2.multipliers) -plot(sol4, show=true) - -# save / load solution in JLD2 format -println("\nSave / load solution in JLD2 format") -save(sol, filename_prefix="./test/solution_test") -sol_reloaded = load("./test/solution_test") -println("Check JLD2 solution ", sol.objective == sol_reloaded.objective) - -# export / read discrete solution in JSON format -# NB. we recover here a JSON Object... -println("Export / import solution in JSON format") -export_ocp_solution(sol, filename_prefix="./test/solution_test") -sol_disc_reloaded = import_ocp_solution("./test/solution_test") -println("Check JSON solution ", sol.objective == sol_disc_reloaded.objective) - -println("") diff --git a/test/test_objective.jl b/test/test_objective.jl deleted file mode 100644 index ae22281..0000000 --- a/test/test_objective.jl +++ /dev/null @@ -1,53 +0,0 @@ -include("deps.jl") - -println("Test: objective") - - -# min tf -ocp = Model(variable=true) -state!(ocp, 2) -control!(ocp, 1) -variable!(ocp, 1) -time!(ocp, t0=0, indf=1) -constraint!(ocp, :initial, lb=[0,0], ub=[0,0]) -constraint!(ocp, :final, lb=[1,0], ub=[1,0]) -constraint!(ocp, :control, lb=-1, ub=1) -constraint!(ocp, :variable, lb=0.1, ub=10) -dynamics!(ocp, (x, u, v) -> [x[2], u]) -objective!(ocp, :mayer, (x0, xf, v) -> v) -sol = solve(ocp, print_level=0, tol=1e-12) -println("Target 2.0, found ", sol.objective) - - -# min tf (lagrange) -ocp = Model(variable=true) -state!(ocp, 2) -control!(ocp, 1) -variable!(ocp, 1) -time!(ocp, t0=0, indf=1) -constraint!(ocp, :initial, lb=[0,0], ub=[0,0]) -constraint!(ocp, :final, lb=[1,0], ub=[1,0]) -constraint!(ocp, :control, lb=-1, ub=1) -constraint!(ocp, :variable, lb=0.05, ub=10) -dynamics!(ocp, (x, u, v) -> [x[2], u]) -objective!(ocp, :lagrange, (x, u, v) -> 1) -sol = solve(ocp, print_level=0, tol=1e-12) -println("Target 2.0, found ", sol.objective) -plot(sol, show=true) - - -# max t0 (free t0 and tf) -ocp = Model(variable=true) -state!(ocp, 2) -control!(ocp, 1) -variable!(ocp, 2) -time!(ocp, ind0=1, indf=2) -constraint!(ocp, :initial, lb=[0,0], ub=[0,0]) -constraint!(ocp, :final, lb=[1,0], ub=[1,0]) -constraint!(ocp, :control, lb=-1, ub=1) -constraint!(ocp, :variable, lb=[0.1,0.1], ub=[10,10]) -constraint!(ocp, :variable, f=v->v[2]-v[1], lb=0.1, ub=Inf) -dynamics!(ocp, (x, u, v) -> [x[2], u]) -objective!(ocp, :mayer, (x0, xf, v) -> v[1], :max) -sol = solve(ocp, print_level=0, tol=1e-12) -println("Target 8.0, found ", sol.objective)