Prepare 0.13 (#331)

benchmark/benchmark.jl

@@ -34,19 +34,20 @@ function bench_list(problem_list; verbose=2, nlp_solver, linear_solver, kwargs..
         if !isnothing(problem[:obj]) && !isapprox(sol.objective, problem[:obj], rtol = 5e-2)
             error("Objective mismatch for ",problem[:name],": ",sol.objective," instead of ",problem[:obj])
         else
-            verbose > 1 && @printf("%-30s completed after %4d iterations\n", problem[:name], sol.iterations)
+            verbose > 1 && @printf("%-30s: %4d iter ", problem[:name], sol.iterations)
         end
 
         # time
         t = @belapsed direct_solve($problem[:ocp], $nlp_solver; init=$problem[:init], display=false, $kwargs...)
         append!(t_list, t)
+        verbose > 1 && @printf("%7.2f s\n", t)
     end
 
     return sum(t_list)
 end
 
 
-function bench(;grid_size_list = [100, 250, 500, 1000], verbose = 1, nlp_solver=:ipopt, linear_solver=nothing, series = :default, kwargs...)
+function bench(;grid_size_list = [250, 500, 1000, 2500], verbose = 1, nlp_solver=:ipopt, linear_solver=nothing, names_list = :default, kwargs...)
 
     #######################################################
     # set (non) linear solvers and backends
@@ -63,8 +64,12 @@ function bench(;grid_size_list = [100, 250, 500, 1000], verbose = 1, nlp_solver=
     verbose > 1 && @printf("Blas config: %s\n", LinearAlgebra.BLAS.lbt_get_config())
 
     # load problems for benchmark
-    if series == :default
-        names_list = ["algal_bacterial", "beam", "fuller", "goddard", "jackson", "vanderpol"]
+    if names_list == :default
+        names_list = ["beam", "double_integrator_mintf", "double_integrator_minenergy", "double_integrator_freet0tf", "fuller", "goddard", "goddard_all", "jackson", "robbins", "simple_integrator", "vanderpol"]
+    elseif names_list == :all 
+        names_list = ["algal_bacterial", "beam", "bioreactor_1day", "bioreactor_Ndays", "bolza_freetf", "double_integrator_mintf", "double_integrator_minenergy", "double_integrator_freet0tf", "fuller", "goddard", "goddard_all", "insurance", "jackson", "robbins", "simple_integrator", "swimmer", "vanderpol"]
+    elseif names_list == :hard
+        names_list = ["algal_bacterial", "bioreactor_1day", "bioreactor_Ndays", "bolza_freetf", "goddard_all", "insurance", "swimmer"]
     end
     println("Problem list: ", names_list)
     problem_list = []

docs/Project.toml

@@ -5,7 +5,6 @@ Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 DocumenterMermaid = "a078cd44-4d9c-4618-b545-3ab9d77f9177"
 HSL = "34c5aeac-e683-54a6-a0e9-6e0fdc586c50"
 NLPModelsIpopt = "f4238b75-b362-5c4c-b852-0801c9a21d71"
-Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 
 [compat]
 CTBase = "0.14"
@@ -14,4 +13,3 @@ Documenter = "1"
 DocumenterMermaid = "0.1"
 HSL = "0.4"
 NLPModelsIpopt = "0.10"
-Plots = "1"

docs/make.jl

@@ -5,7 +5,6 @@ using CTBase
 
 using NLPModelsIpopt
 using HSL
-using Plots
 
 # to add docstrings from external packages
 DocMeta.setdocmeta!(CTBase, :DocTestSetup, :(using CTBase); recursive = true)

src/disc/irk.jl

@@ -6,18 +6,8 @@ Internal layout for NLP variables:
  X_N-1, U_N-1, K_N-1^1..K_N-1^s,
  X_N, U_N, V]
 with s the stage number and U given by linear interpolation in [t_i, t_i+1]
-NB. +++ could use constant interpolation for 1-stage methods (but U_N might end up unused)
+NB. 1-stage methods use constant interpolation instead (but U_N might end up unused)
 Path constraints are all evaluated at time steps
-
-+++LATER: full stage version with both state and control discretized at stage times !
-path constraints are evaluated at stage times; add x0 and xf for boundary conditions
-[X0, X01..X0s, U01..U0s, K01..K0s, ..., XN-11..XN-1s, UN-11..UN-1s, KN-11..KN-1s, XN, V]
-add stage_grid to time_grid since all interpolations will now use the stages instead of steps
-notes: 
-- Kij are just f(Xij) and may be omitted (smaller problem but maybe more nonlinear) ?
-- however we have to enforce the equations on the Xij in which the Kij appear.
-- avoid recomputing f(Xij) multiple times, use work array if needed
-Start with Midpoint_fullstage in midpoint.jl ?
 =#
 
 

src/disc/midpoint.jl

@@ -1,8 +1,9 @@
 #= Functions for implicit midpoint discretization scheme
 Internal layout for NLP variables: 
-[X_1,U_1, .., X_N,U_N, X_N+1, V]
+[X_1,U_1,K_1 .., X_N,U_N,K_N, X_N+1, V]
 with the convention u([t_i,t_i+1[) = U_i and u(tf) = U_N
-NB. stage variables are removed via the simplification x_s = (x_i + x_i+1) / 2
+NB. stage variables K_i can be removed via the simplification x_s = (x_i + x_i+1) / 2
+however this seems to give worse performance...
 =#
 
 struct Midpoint <: Discretization
@@ -11,19 +12,13 @@ struct Midpoint <: Discretization
     _step_variables_block::Int
     _state_stage_eqs_block::Int
     _step_pathcons_block::Int
-    _kvars::Bool
 
     # constructor
-    function Midpoint(dim_NLP_steps, dim_NLP_x, dim_NLP_u, dim_NLP_v, dim_u_cons, dim_x_cons, dim_xu_cons, dim_boundary_cons, dim_v_cons; kvars=false)
+    function Midpoint(dim_NLP_steps, dim_NLP_x, dim_NLP_u, dim_NLP_v, dim_u_cons, dim_x_cons, dim_xu_cons, dim_boundary_cons, dim_v_cons)
 
         # aux variables
-        if kvars
-            step_variables_block = dim_NLP_x * 2 + dim_NLP_u
-            state_stage_eqs_block = dim_NLP_x * 2
-        else
-            step_variables_block = dim_NLP_x + dim_NLP_u
-            state_stage_eqs_block = dim_NLP_x
-        end
+        step_variables_block = dim_NLP_x * 2 + dim_NLP_u
+        state_stage_eqs_block = dim_NLP_x * 2
 
         # NLP variables size ([state, control]_1..N, final state, variable)
         dim_NLP_variables = dim_NLP_steps * step_variables_block + dim_NLP_x + dim_NLP_v
@@ -34,7 +29,7 @@ struct Midpoint <: Discretization
         # NLP constraints size ([dynamics, path]_1..N, final path, boundary, variable)
         dim_NLP_constraints = dim_NLP_steps * (state_stage_eqs_block + step_pathcons_block) + step_pathcons_block + dim_boundary_cons + dim_v_cons
 
-        disc = new("Implicit Midpoint aka Gauss-Legendre collocation for s=1, 2nd order, symplectic", step_variables_block, state_stage_eqs_block, step_pathcons_block, kvars)
+        disc = new("Implicit Midpoint aka Gauss-Legendre collocation for s=1, 2nd order, symplectic", step_variables_block, state_stage_eqs_block, step_pathcons_block)
 
         return disc, dim_NLP_variables, dim_NLP_constraints
     end
@@ -169,39 +164,21 @@ function setStepConstraints!(docp::DOCP{Midpoint}, c, xu, v, time_grid, i, work)
             xs = work
         end
 
-        if docp.discretization._kvars
-            # state equation: midpoint rule
-            ki = get_stagevars_at_time_step(xu, docp, i, 1)
-            @views @. c[offset+1:offset+docp.dim_OCP_x] = xip1 - (xi + hi * ki[1:docp.dim_OCP_x])
-            if docp.is_lagrange
-            c[offset+docp.dim_NLP_x] = get_lagrange_state_at_time_step(xu, docp, i+1) - (get_lagrange_state_at_time_step(xu, docp, i) + hi * ki[docp.dim_NLP_x])
-            end
-            offset += docp.dim_NLP_x
-
-            # stage equation at mid step k_i = f(x_s)
-            docp.ocp.dynamics((@view c[offset+1:offset+docp.dim_OCP_x]), ts, xs, ui, v)
-            if docp.is_lagrange
-                docp.ocp.lagrange((@view c[offset+docp.dim_NLP_x:offset+docp.dim_NLP_x]), ts, xs, ui, v)
-            end            
-            @views @. c[offset+1:offset+docp.dim_NLP_x] = ki - c[offset+1:offset+docp.dim_NLP_x]
-            offset += docp.dim_NLP_x
-
-        else
-            # No stage variables, compute dynamics directly for state equation
-            # OCP dynamics
-            docp.ocp.dynamics((@view c[offset+1:offset+docp.dim_OCP_x]), ts, xs, ui, v)
-            # lagrange cost
-            if docp.is_lagrange
-                docp.ocp.lagrange((@view c[offset+docp.dim_NLP_x:offset+docp.dim_NLP_x]), ts, xs, ui, v)
-            end
-
-            # midpoint rule
-            @views @. c[offset+1:offset+docp.dim_OCP_x] = xip1 - (xi + hi * c[offset+1:offset+docp.dim_OCP_x])
-            if docp.is_lagrange
-            c[offset+docp.dim_NLP_x] = get_lagrange_state_at_time_step(xu, docp, i+1) - (get_lagrange_state_at_time_step(xu, docp, i) + hi * c[offset+docp.dim_NLP_x])
-            end
-            offset += docp.dim_NLP_x
+        # state equation: midpoint rule
+        ki = get_stagevars_at_time_step(xu, docp, i, 1)
+        @views @. c[offset+1:offset+docp.dim_OCP_x] = xip1 - (xi + hi * ki[1:docp.dim_OCP_x])
+        if docp.is_lagrange
+        c[offset+docp.dim_NLP_x] = get_lagrange_state_at_time_step(xu, docp, i+1) - (get_lagrange_state_at_time_step(xu, docp, i) + hi * ki[docp.dim_NLP_x])
         end
+        offset += docp.dim_NLP_x
+
+        # stage equation at mid step k_i = f(x_s)
+        docp.ocp.dynamics((@view c[offset+1:offset+docp.dim_OCP_x]), ts, xs, ui, v)
+        if docp.is_lagrange
+            docp.ocp.lagrange((@view c[offset+docp.dim_NLP_x:offset+docp.dim_NLP_x]), ts, xs, ui, v)
+        end            
+        @views @. c[offset+1:offset+docp.dim_NLP_x] = ki - c[offset+1:offset+docp.dim_NLP_x]
+        offset += docp.dim_NLP_x
 
     end
 

src/docp.jl

@@ -177,9 +177,7 @@ struct DOCP{T <: Discretization, X <: ScalVect, U <: ScalVect, V <: ScalVect}
         if disc_method == :trapeze
             discretization, dim_NLP_variables, dim_NLP_constraints = CTDirect.Trapeze(dim_NLP_steps, dim_NLP_x, dim_NLP_u, dim_NLP_v, dim_u_cons, dim_x_cons, dim_xu_cons, dim_boundary_cons, dim_v_cons)
         elseif disc_method == :midpoint
-            discretization, dim_NLP_variables, dim_NLP_constraints = CTDirect.Midpoint(dim_NLP_steps, dim_NLP_x, dim_NLP_u, dim_NLP_v, dim_u_cons, dim_x_cons, dim_xu_cons, dim_boundary_cons, dim_v_cons)          
-        elseif disc_method == :midpoint_kvars
-            discretization, dim_NLP_variables, dim_NLP_constraints = CTDirect.Midpoint(dim_NLP_steps, dim_NLP_x, dim_NLP_u, dim_NLP_v, dim_u_cons, dim_x_cons, dim_xu_cons, dim_boundary_cons, dim_v_cons; kvars=true)
+            discretization, dim_NLP_variables, dim_NLP_constraints = CTDirect.Midpoint(dim_NLP_steps, dim_NLP_x, dim_NLP_u, dim_NLP_v, dim_u_cons, dim_x_cons, dim_xu_cons, dim_boundary_cons, dim_v_cons)
         elseif disc_method == :gauss_legendre_1
                 discretization, dim_NLP_variables, dim_NLP_constraints = CTDirect.Gauss_Legendre_1(dim_NLP_steps, dim_NLP_x, dim_NLP_u, dim_NLP_v, dim_u_cons, dim_x_cons, dim_xu_cons, dim_boundary_cons, dim_v_cons)
         elseif disc_method == :gauss_legendre_2

test/problems/bolza.jl

@@ -6,6 +6,8 @@ function bolza_freetf()
         t ∈ [0, tf], time
         x ∈ R, state
         u ∈ R, control
+        tf >= 0.1
+        x(t) >= 0
         ẋ(t) == tf * u(t)
         x(0) == 0
         x(tf) == 1

test/problems/double_integrator.jl

@@ -20,7 +20,7 @@ end
 
 
 # min energy with fixed tf
-function double_integrator_minenergy(T)
+function double_integrator_minenergy(T=2)
     @def ocp begin
         t ∈ [0, T], time
         x ∈ R², state

test/problems/swimmer.jl

@@ -2,9 +2,8 @@
 
 # +++ make 2 versions: 1 stroke periodic and free N strokes
 
-function swimmer()
+function swimmer(tf = 25)
     @def swimmer begin
-        tf = 25
         t ∈ [0, tf], time
         x ∈ R^5, state
         u ∈ R^2, control
@@ -31,8 +30,8 @@ function swimmer()
         th = x[3](t)
         b1 = x[4](t)
         b3 = x[5](t)
-        u1 = u[1](t)
-        u2 = u[2](t)
+        a1 = u[1](t)
+        a2 = u[2](t)
 
         aux =
             543 +
@@ -126,9 +125,9 @@ function swimmer()
                 2 * cos(2 * b3) - 6 * cos(b1 + b3) - 6 * cos(2 * b1 + b3)
             ) / (2 * aux)
 
-        ẋ(t) == [g11 * u1 + g21 * u2, g12 * u1 + g22 * u2, g13 * u1 + g23 * u2, u1, u2]
+        ẋ(t) == [g11 * a1 + g21 * a2, g12 * a1 + g22 * a2, g13 * a1 + g23 * a2, a1, a2]
         x[1](tf) → max
-        #∫(u1^2 + u2^2) → min
+        #∫(a1^2 + a2^2) → min
     end
 
     return ((ocp = swimmer, obj = 0.984273, name = "swimmer", init = nothing))

test/suite/test_all_ocp.jl

@@ -1,4 +1,4 @@
-println("Test: OCP definitions")
+println("testing: OCP definitions")
 
 # beam
 if !isdefined(Main, :beam)

test/suite/test_objective.jl

@@ -33,7 +33,7 @@ if !isdefined(Main, :bolza_freetf)
 end
 prob = bolza_freetf()
 @testset verbose = true showtiming = true ":bolza :tf_in_dyn_and_cost" begin
-    sol = direct_solve(prob.ocp, display = false, grid_size=100)
+    sol = direct_solve(prob.ocp, display = false)
     @test sol.objective ≈ prob.obj rtol = 1e-2
 end