Merge branch 'main' into 222-manipulation-of-the-nlp

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/make.jl

@@ -20,6 +20,7 @@ makedocs(
         prettyurls = false,
         size_threshold_ignore = [
             "api-ctbase/types.md",
+            "tutorial-plot.md",
         ],
         assets = [
             asset("https://control-toolbox.org/assets/css/documentation.css"),

docs/src/tutorial-iss.md

@@ -1,4 +1,4 @@
-# Indirect simple shooting
+# [Indirect simple shooting](@id iss)
 
 In this tutorial we present the indirect simple shooting method on a simple example.
 
@@ -160,7 +160,16 @@ end # hide
 nothing # hide
 ```
 
-We get:
+The solution can be plot calling first the flow.
+
+```@example main
+sol = φ( (t0, tf), x0, p0_sol )
+plot(sol, size=(800, 600))
+```
+
+In the indirect shooting method, the research of the optimal control is replaced by the computation
+of its associated extremal. This computation is equivalent to finding the initial covector solution
+to the shooting function.
 
 ```@example main
 exp(p0; saveat=[]) = φ((t0, tf), x0, p0, saveat=saveat).ode_sol # hide

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. 
@@ -48,16 +48,11 @@ nothing # hide
 We 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`.
 
 We can now use the solver of our choice to solve it.
 
@@ -67,8 +62,7 @@ For a first example we use the `ipopt` solver from [`NLPModelsIpopt.jl`](https:/
 
 ```@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
 ```
 
@@ -84,14 +78,14 @@ 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
 ```
 
@@ -101,7 +95,6 @@ with as initial guess the solution from the first resolution.
 ```@example main
 using Percival
 
-nlp = get_nlp(docp)
 output = percival(nlp)
 print(output)
 ```

docs/src/tutorial-plot.md

@@ -66,6 +66,33 @@ plot(sol;
      control_style = (color=:red, linewidth=2))
 ```
 
+## From Flow
+
+The previous resolution of the optimal control problem was done with the `solve` function.
+If you use an indirect shooting method and solve shooting equations, you may want to plot the 
+associated solution. To do so, you need to use the `Flow` function to 
+reconstruct the solution. See the [Indirect Simple Shooting](@ref iss) tutorial for an example.
+In our example, you must provide the maximising control $(x, p) \mapsto p_2$ together with the 
+optimal control problem.
+
+```@example main
+t0 = 0
+tf = 1
+x0 = [ -1, 0 ]
+p0 = sol.costate(t0)
+f  = Flow(ocp, (x, p) -> p[2])
+sol_flow = f( (t0, tf), x0, p0 )
+plot(sol_flow)
+```
+
+You can notice that the time grid has very few points. To have a better visualisation (the accuracy 
+won't change), you can give a finer grid.
+
+```@example main
+sol_flow = f( (t0, tf), x0, p0; saveat=range(t0, tf, 100) )
+plot(sol_flow)
+```
+
 ## Split versus group layout
 
 If you prefer to get a more compact figure, you can use the `layout` optional keyword argument with `:group` value. It will group the state, costate and control trajectories in one subplot for each.

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