Merge pull request #252 from control-toolbox/0.10.0

0.10.0

Author: PierreMartinon

Project.toml

@@ -1,7 +1,7 @@
 name = "OptimalControl"
 uuid = "5f98b655-cc9a-415a-b60e-744165666948"
 authors = ["Olivier Cots <olivier.cots@toulouse-inp.fr>"]
-version = "0.9.4"
+version = "0.10.0"
 
 [deps]
 CTBase = "54762871-cc72-4466-b8e8-f6c8b58076cd"
@@ -12,7 +12,7 @@ DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 
 [compat]
 CTBase = "0.11"
-CTDirect = "0.9"
+CTDirect = "0.10"
 CTFlows = "0.4"
 CommonSolve = "0.2"
 DocStringExtensions = "0.9"

docs/src/tutorial-nlp.md

@@ -8,7 +8,7 @@ We describe here some more advanced operations related to the discretized optima
 When calling `solve(ocp)` three steps are performed internally:
 
 - first, the OCP is discretized into a DOCP (a nonlinear optimization problem) with [`direct_transcription`](@ref),
-- then, this DOCP is solved,
+- then, this DOCP is solved (with the internal function `solve_docp`),
 - finally, a functional solution of the OCP is rebuilt from the solution of the discretized problem, with [`build_solution`](@ref).
 
 These steps can also be done separately, for instance if you want to use your own NLP solver. Let us load the modules
@@ -37,23 +37,18 @@ nothing # hide
 First let us discretize the problem.
 
 ```@example main
-docp = direct_transcription(ocp)
+docp, nlp = direct_transcription(ocp)
 nothing # hide
 ```
 
-The DOCP contains a copy of the original OCP, and the resulting discretized problem, in our case an `ADNLPModel`.
-You can extract this raw NLP problem  with the [`get_nlp`](@ref) method.
-
-```@example main
-nlp = get_nlp(docp)
-```
+The DOCP contains information related to the transcription, including a copy of the original OCP, and the NLP is the resulting discretized problem, in our case an `ADNLPModel`.
 
 You could then use the solver of your choice to solve it.
 For an example we use the `ipopt` solver from [`NLPModelsIpopt.jl`](https://github.com/JuliaSmoothOptimizers/NLPModelsIpopt.jl) package to solve the NLP problem.
 
 ```@example main
 using NLPModelsIpopt
-nlp_sol = ipopt(get_nlp(docp); print_level=4, mu_strategy="adaptive", tol=1e-8, sb="yes")
+nlp_sol = ipopt(nlp; print_level=4, mu_strategy="adaptive", tol=1e-8, sb="yes")
 nothing # hide
 ```
 
@@ -67,13 +62,13 @@ plot(sol)
 An initial guess, including warm start, can be passed to [`direct_transcription`](@ref) the same way as for `solve`.
 
 ```@example main
-docp = direct_transcription(ocp, init=sol)
+docp, nlp = direct_transcription(ocp, init=sol)
 nothing # hide
 ```
 
 It can also be changed after the transcription is done, with  [`set_initial_guess`](@ref).
 
 ```@example main
-set_initial_guess(docp, sol)
+set_initial_guess(docp, nlp, sol)
 nothing # hide
 ```

src/OptimalControl.jl

@@ -42,15 +42,13 @@ export *
 
 # CTDirect
 export direct_transcription
-export get_nlp
 export set_initial_guess
 export build_solution
 export save
 export load
 export export_ocp_solution
 export import_ocp_solution
 
-
 # CTBase
 export Index
 export Autonomous, NonAutonomous

src/solve.jl

@@ -66,14 +66,6 @@ julia> sol = solve(ocp, init=(state=t->[-1+t, t*(t-1)], control=t->6-12*t))
 ```
 
 """
-#=
-function solve(ocp::OptimalControlModel, description::Symbol...; 
-    display::Bool=__display(),
-    init=__ocp_init(),
-    kwargs...)
-    =#
-
-
 function solve(ocp::OptimalControlModel, description::Symbol...;
     init=__ocp_init(),
     grid_size::Integer=CTDirect.__grid_size_direct(),
@@ -88,17 +80,16 @@ function solve(ocp::OptimalControlModel, description::Symbol...;
 
     # print chosen method
     method = getFullDescription(description, available_methods())
-    display ? println("Method = ", method) : nothing
+    #display ? println("Method = ", method) : nothing
 
     # if no error before, then the method is correct: no need of else
     if :direct ∈ method
-        #return CTDirect.solve(ocp, clean(method)...; display=display, init=init, kwargs...)
 
         # build discretized OCP
-        docp = direct_transcription(ocp, description, init=init, grid_size=grid_size, time_grid=time_grid)
+        docp, nlp = direct_transcription(ocp, description, init=init, grid_size=grid_size, time_grid=time_grid)
 
         # solve DOCP (NB. init is already embedded in docp)
-        docp_solution = CTDirect.solve_docp(docp, display=display, print_level=print_level, mu_strategy=mu_strategy, tol=tol, max_iter=max_iter, linear_solver=linear_solver; kwargs...)
+        docp_solution = CTDirect.solve_docp(docp, nlp,  display=display, print_level=print_level, mu_strategy=mu_strategy, tol=tol, max_iter=max_iter, linear_solver=linear_solver; kwargs...)
 
         # build and return OCP solution
         return build_solution(docp, docp_solution)

test/problems/beam.jl

@@ -0,0 +1,18 @@
+# Beam example from bocop
+
+function 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

test/problems/beam.png

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

test/problems/bioreactor.jl

@@ -0,0 +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
+
+    return((ocp=fuller, obj=2.683944e-1, name="fuller", init=nothing))
+end