problem library now defines functions (#186)

* problem library now defines functions that return named tuples for problems and relevant data, including init

* updated tests

benchmark/benchmark.jl

@@ -17,23 +17,25 @@ blas_config = LinearAlgebra.BLAS.lbt_get_config()
 #######################################################
 # set parameters
 tol = 1e-8
-grid_size = 5000
+grid_size = 500
 precompile = true
 @printf("Settings: tol=%g grid_size=%d precompile=%s\n\n", tol, grid_size, precompile)
 
 #######################################################
-# load examples
-# +++ could add :init field in the problems !
+# load examples library
+problem_path = pwd()*"/problems"
+for problem_file in filter(contains(r".jl$"), readdir(problem_path; join=true))
+    include(problem_file)
+end
+
+# load problems for benchmark
 print("Loading problems: ")
 names_list = ["beam", "bioreactor_1day_periodic", "fuller", "goddard", "insurance", "jackson", "swimmer", "vanderpol"]
 problem_list = []
-problem_path = pwd()*"/problems"
-for problem_file in filter(contains(r".jl$"), readdir(problem_path; join=true))
-    ocp_data = include(problem_file)
-    if ocp_data.name in names_list
-        @printf("%s ", ocp_data.name)
-        push!(problem_list,ocp_data)
-    end
+for problem_name in names_list
+    ocp_data = getfield(Main, Symbol(problem_name))()
+    @printf("%s ", ocp_data.name)
+    push!(problem_list,ocp_data)
 end
 println("")
 
@@ -55,8 +57,8 @@ end
 t_list = []
 println("Benchmark step")
 for problem in problem_list
-    t = @elapsed local sol = solve(problem[:ocp], display=false, linear_solver=linear_solver, grid_size=grid_size, tol=tol)
-    if !isapprox(sol.objective, problem[:obj], rtol=5e-2)
+    t = @elapsed local sol = solve(problem[:ocp], init=problem[:init], display=false, linear_solver=linear_solver, grid_size=grid_size, tol=tol)
+    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
         @printf("%-30s completed in %6.2f s after %4d iterations\n",problem[:name],t,sol.iterations)

problems/beam.jl

@@ -1,15 +1,18 @@
 # Beam example from bocop
 
-@def beam begin
-    t ∈ [ 0, 1 ], time
-    x ∈ R², state
-    u ∈ R, control
-    x(0) == [0,1]
-    x(1) == [0,-1]
-    ẋ(t) == [x₂(t), u(t)]
-    0 ≤ x₁(t) ≤ 0.1    
-    -10 ≤ u(t) ≤ 10
-    ∫(u(t)^2) → min
-end
+function beam()
 
-return ((ocp=beam, obj=8.898598, name="beam"))
+    @def beam begin
+        t ∈ [ 0, 1 ], time
+        x ∈ R², state
+        u ∈ R, control
+        x(0) == [0,1]
+        x(1) == [0,-1]
+        ẋ(t) == [x₂(t), u(t)]
+        0 ≤ x₁(t) ≤ 0.1    
+        -10 ≤ u(t) ≤ 10
+        ∫(u(t)^2) → min
+    end
+
+    return ((ocp=beam, obj=8.898598, name="beam", init=nothing))
+end

problems/bioreactor.jl

@@ -0,0 +1,96 @@
+# Bioreactor example from bocop
+
+# METHANE PROBLEM
+# mu2 according to growth model
+# mu according to light model
+# time scale is [0,10] for 24h (day then night)
+
+# growth model MONOD
+function growth(s, mu2m, Ks)
+    return mu2m * s / (s+Ks)
+end
+
+# light model: max^2 (0,sin) * mubar
+# DAY/NIGHT CYCLE: [0,2 halfperiod] rescaled to [0,2pi]
+function light(time, halfperiod)
+    pi = 3.141592653589793
+    days = time / (halfperiod*2)
+    tau = (days - floor(days)) * 2*pi
+    return max(0,sin(tau))^2
+end
+
+# 1 day periodic problem
+function bioreactor_1day_periodic()
+    @def bioreactor_1 begin
+        # constants
+        beta = 1
+        c = 2
+        gamma = 1
+        Ks = 0.05
+        mu2m = 0.1
+        mubar = 1
+        r = 0.005
+        halfperiod = 5
+        T = halfperiod * 2
+
+        # ocp
+        t ∈ [ 0, T ], time
+        x ∈ R³, state
+        u ∈ R, control
+        y = x[1]
+        s = x[2]
+        b = x[3]
+        mu = light(t, halfperiod) * mubar
+        mu2 = growth(s(t), mu2m, Ks)
+        [0,0,0.001] ≤ x(t) ≤ [Inf, Inf, Inf]
+        0 ≤ u(t) ≤ 1
+        1 ≤ y(0) ≤ Inf
+        1 ≤ b(0) ≤ Inf
+        x(0)- x(T) == [0,0,0]
+        ẋ(t) == [mu*y(t)/(1+y(t))-(r+u(t))*y(t),
+                -mu2*b(t) + u(t)*beta*(gamma*y(t)-s(t)),
+                (mu2-u(t)*beta)*b(t)]
+        ∫(mu2*b(t)/(beta+c)) → max
+    end
+
+    return ((ocp=bioreactor_1, obj=0.614134, name="bioreactor_1day_periodic", init=nothing))
+end
+
+
+# N days (non periodic)
+function bioreactor_Ndays()
+    @def bioreactor_N begin
+        # constants
+        beta = 1
+        c = 2
+        gamma = 1
+        halfperiod = 5
+        Ks = 0.05
+        mu2m = 0.1
+        mubar = 1
+        r = 0.005
+        N = 30
+        T = 10 * N
+
+        # ocp
+        t ∈ [ 0, T ], time
+        x ∈ R³, state
+        u ∈ R, control
+        y = x[1]
+        s = x[2]
+        b = x[3]
+        mu = light(t, halfperiod) * mubar
+        mu2 = growth(s(t), mu2m, Ks)
+        [0,0,0.001] ≤ x(t) ≤ [Inf, Inf, Inf]
+        0 ≤ u(t) ≤ 1
+        0.05 ≤ y(0) ≤ 0.25
+        0.5 ≤ s(0) ≤ 5
+        0.5 ≤ b(0) ≤ 3
+        ẋ(t) == [mu*y(t)/(1+y(t))-(r+u(t))*y(t),
+                -mu2*b(t) + u(t)*beta*(gamma*y(t)-s(t)),
+                (mu2-u(t)*beta)*b(t)]
+        ∫(mu2*b(t)/(beta+c)) → max
+    end
+
+    return ((ocp=bioreactor_N, obj=19.0745, init=(state=[50,50,50],), name="bioreactor_Ndays"))
+end

problems/fuller.jl

@@ -1,17 +1,16 @@
 # Fuller example
+function fuller()
 
-@def fuller begin
-    t ∈ [ 0, 3.5 ], time
-    x ∈ R², state
-    u ∈ R, control
-    -1 ≤ u(t) ≤ 1
-    x(0) == [ 0, 1 ]
-    x(3.5) == [ 0, 0 ]
-    ẋ(t) == [ x₂(t), u(t) ] 
-    ∫(x₁(t)^2) → min
-end
+    @def fuller begin
+        t ∈ [ 0, 3.5 ], time
+        x ∈ R², state
+        u ∈ R, control
+        -1 ≤ u(t) ≤ 1
+        x(0) == [ 0, 1 ]
+        x(3.5) == [ 0, 0 ]
+        ẋ(t) == [ x₂(t), u(t) ] 
+        ∫(x₁(t)^2) → min
+    end
 
-#sol = solve(fuller, grid_size=1000)
-#plot(sol)
-
-return((ocp=fuller, obj=2.683944e-1, name="fuller"))
+    return((ocp=fuller, obj=2.683944e-1, name="fuller", init=nothing))
+end