Merge branch 'dev' into revise

Author: ocots

Project.toml

@@ -1,22 +1,16 @@
 name = "OptimalControl"
 uuid = "5f98b655-cc9a-415a-b60e-744165666948"
 authors = ["Olivier Cots <[email protected]>"]
-version = "0.5.0"
+version = "0.6.0"
 
 [deps]
 CTBase = "54762871-cc72-4466-b8e8-f6c8b58076cd"
 CTDirect = "790bbbee-bee9-49ee-8912-a9de031322d5"
-CTDirectShooting = "e0e6c04b-5022-4cd2-bea2-4a09fff39444"
 CTFlows = "1c39547c-7794-42f7-af83-d98194f657c2"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
-MINPACK = "4854310b-de5a-5eb6-a2a5-c1dee2bd17f9"
-Markdown = "d6f4376e-aef5-505a-96c1-9c027394607a"
-Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
-Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 
 [compat]
-CTBase = "0.5"
-CTDirect = "0.2"
-CTDirectShooting = "0.1"
-Reexport = "1.2"
+CTBase = "0.7"
+CTDirect = "0.3"
+CTFlows = "0.3"
 julia = "1.8"

examples/basic.jl

@@ -58,9 +58,9 @@ dynamics!(ocp_f, (x, u) -> A*x + B*u)
 objective!(ocp_f, :lagrange, (x, u) -> 0.5u^2)
 
 # Solve
-ocp = ocp_f
+ocp = ocp_a
 N = 50
-sol = solve(ocp_f, grid_size=N)
+sol = solve(ocp, grid_size=N)
 
 # Plot
 plot(sol)

examples/basic_fig.png

@@ -19,6 +19,8 @@ v0 = 0
 vmax = 0.1
 m0 = 1
 mf = 0.6
+
+# Initial state
 x0 = [ r0, v0, m0 ]
 
 # Abstract model
@@ -35,9 +37,9 @@ x0 = [ r0, v0, m0 ]
 
     x(t0) == [ r0, v0, m0 ]
     0  ≤ u(t) ≤ 1
-    r0  ≤ r(t) ≤ Inf,    (1)
-    0  ≤ v(t) ≤ vmax,    (2)
-    mf ≤ m(t) ≤ m0,      (3)
+    r0 ≤ r(t) ≤ Inf,    (1)
+    0  ≤ v(t) ≤ vmax,   (2)
+    mf ≤ m(t) ≤ m0,     (3)
 
     ẋ(t) == F0(x(t)) + u(t) * F1(x(t))
 
@@ -48,14 +50,12 @@ end
 F0(x) = begin
     r, v, m = x
     D = Cd * v^2 * exp(-β*(r - 1))
-    F = [ v, -D/m - 1/r^2, 0 ]
-    return F
+    return [ v, -D/m - 1/r^2, 0 ]
 end
 
 F1(x) = begin
     r, v, m = x
-    F = [ 0, Tmax/m, -b*Tmax ]
-    return F
+    return [ 0, Tmax/m, -b*Tmax ]
 end
 
 # Functional model
@@ -87,27 +87,27 @@ savefig("goddard_fig1.png")
 md"![fig](goddard_fig1.png)"
 
 # Shooting function
-u0(x, p) = 0.
-u1(x, p) = 1.
-
-H0(x, p) = p' * F0(x)
-H1(x, p) = p' * F1(x)
-H01 = Poisson(H0, H1)
-H001 = Poisson(H0, H01)
-H101 = Poisson(H1, H01)
+u0 = 0
+u1 = 1
+
+H0 = Lift(F0)
+H1 = Lift(F1)
+H01  = @Poisson {H0, H1}
+H001 = @Poisson {H0, H01}
+H101 = @Poisson {H1, H01}
 us(x, p) = -H001(x, p) / H101(x, p)
 
-remove_constraint!(ocp, :eq1)
-remove_constraint!(ocp, :eq3)
-constraint!(ocp, :final, Index(3), mf, :eq4)
-g(x) = vmax - constraint(ocp, :eq2)(x)
-ub(x, _) = -Ad(F0, g)(x) / Ad(F1, g)(x)
-μb(x, p) = H01(x, p) / Ad(F1, g)(x)
+#remove_constraint!(ocp, :eq1)
+#remove_constraint!(ocp, :eq3)
+#constraint!(ocp, :final, Index(3), mf, :eq4)
+g(x) = vmax - constraint(ocp, :eq2)(x, Real[])
+ub(x)   = -Lie(F0, g)(x) / Lie(F1, g)(x)
+μ(x, p) = H01(x, p) / Lie(F1, g)(x)
 
-f0 = Flow(ocp, u0)
-f1 = Flow(ocp, u1)
-fs = Flow(ocp, us)
-fb = Flow(ocp, ub, (x, _) -> g(x), μb)
+f0 = Flow(ocp, (x, p, tf) -> u0)
+f1 = Flow(ocp, (x, p, tf) -> u1)
+fs = Flow(ocp, (x, p, tf) -> us(x, p))
+fb = Flow(ocp, (x, p, tf) -> ub(x), (x, u, tf) -> g(x), (x, p, tf) -> μ(x, p))
 
 shoot!(s, p0, t1, t2, t3, tf) = begin
 
@@ -129,11 +129,10 @@ t = direct_sol.times
 x = direct_sol.state
 u = direct_sol.control
 p = direct_sol.adjoint
-H1(t) = H1(x(t), p(t))
 
-u_plot = plot(t, t -> u(t)[1], xlabel = "t", ylabel = "u", legend = false)
-H1_plot = plot(t, H1, xlabel = "t", ylabel = "H1", legend = false)
-g_plot = plot(t, g ∘ x, xlabel = "t", ylabel = "g", legend = false)
+u_plot  = plot(t, u,                    label = "u(t)")
+H1_plot = plot(t, t -> H1(x(t), p(t)),  label = "H₁(x(t), p(t))")
+g_plot  = plot(t, g ∘ x,                label = "g(x(t))")
 display(plot(u_plot, H1_plot, g_plot, layout=(3,1), size=(400,300)))
 savefig("goddard_fig2.png")
 md"![fig](goddard_fig2.png)"