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 → 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

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