fixed insurance picture; more tests added

PierreMartinon · PierreMartinon · commit d203fbd3e01c · 2024-07-31T09:58:14.000+02:00
diff --git a/test/problems/fuller.jl.4148.mem b/test/problems/fuller.jl.4148.mem
diff --git a/test/problems/insurance.png b/test/problems/insurance.png
diff --git a/test/runtests.jl b/test/runtests.jl
@@ -13,7 +13,9 @@ include("problems/goddard.jl")
 #
 @testset verbose = true showtiming = true "Optimal control tests" begin
     for name ∈ (
+        :abstract_ocp,
         :basic,
+        :continuation,
         :goddard_direct,
         :goddard_indirect,
         :grid,
diff --git a/test/test_abstract_ocp.jl b/test/test_abstract_ocp.jl
@@ -0,0 +1,29 @@
+# for parser testing
+function test_abstract_ocp()
+
+# 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 
+    ẋ(t) == [ x₂(t), u(t) ] 
+    tf → min
+end
+
+@testset verbose = true showtiming = true ":double_integrator :min_tf :abstract" begin
+    sol1 = solve(ocp1, print_level=0, tol=1e-12)
+    @test sol1.objective ≈ 2.0 rtol=1e-2
+end
+
+@testset verbose = true showtiming = true ":goddard :max_rf :abstract :constr" begin
+    ocp3 = goddard_a().ocp
+    sol3 = solve(ocp3, print_level=0, tol=1e-12)    
+    @test sol3.objective ≈ 1.0125 rtol=1e-2
+end
+
+end
diff --git a/test/test_continuation.jl b/test/test_continuation.jl
@@ -0,0 +1,107 @@
+# discrete continuation
+function test_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
+            ẋ(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
+            ẋ(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
+
+
+# parametric ocp definition
+if test3
+    sol0 = solve(goddard().ocp, print_level=0)
+
+    @testset verbose = true showtiming = true ":global_variable :warm_start" begin
+        sol = sol0
+        Tmax_list = []
+        obj_list = []
+        for Tmax=3.5:-0.5:1 
+            sol = solve(goddard(Tmax=Tmax).ocp, 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
+            # 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
+
+end
diff --git a/test/test_grid.jl b/test/test_grid.jl
@@ -16,4 +16,84 @@ function test_grid()
     sol = solve(ocp; time_grid=time_grid, display=false)
     @test sol.objective ≈ 0.309 rtol=1e-2
 
+
+    # 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,101)
+        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,101), 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
+        ẋ(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,101),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
+
+
 end
