integrate CTBase update

Author: ocots

Project.toml

@@ -1,7 +1,7 @@
 name = "OptimalControl"
 uuid = "5f98b655-cc9a-415a-b60e-744165666948"
 authors = ["Olivier Cots <[email protected]>"]
-version = "0.2.2"
+version = "0.3.0"
 
 [deps]
 CTBase = "54762871-cc72-4466-b8e8-f6c8b58076cd"
@@ -11,7 +11,7 @@ HamiltonianFlows = "5fb78580-10a0-4606-82bc-07a60f425ab3"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 
 [compat]
-CTBase = "0.2"
+CTBase = "0.4"
 CTDirect = "0.1"
 CTDirectShooting = "0.1"
 HamiltonianFlows = "1.0"

examples/basic_f.ipynb

@@ -18,8 +18,8 @@ A = [ 0 1
       0 0 ]
 B = [ 0
       1 ]
-constraint!(ocp, :dynamics, (x, u) -> A*x + B*u[1])
-objective!(ocp, :lagrange, (x, u) -> 0.5u[1]^2)
+constraint!(ocp, :dynamics, (x, u) -> A*x + B*u)
+objective!(ocp, :lagrange, (x, u) -> 0.5u^2)
 
 # replace default callback
 function myprint(i, sᵢ, dᵢ, Uᵢ, gᵢ, fᵢ)
@@ -52,4 +52,4 @@ function myprint3(i, sᵢ, dᵢ, xᵢ, gᵢ, fᵢ)
     old_f = fᵢ
 end
 cbs = (PrintCallback(myprint3, priority=0),)
-sol = solve(ocp, :direct, :shooting, callbacks=cbs);
+sol = solve(ocp, :direct, :shooting, callbacks=cbs)

examples/callbacks-print.jl

@@ -18,8 +18,8 @@ A = [ 0 1
       0 0 ]
 B = [ 0
       1 ]
-constraint!(ocp, :dynamics, (x, u) -> A*x + B*u[1])
-objective!(ocp, :lagrange, (x, u) -> 0.5u[1]^2)
+constraint!(ocp, :dynamics, (x, u) -> A*x + B*u)
+objective!(ocp, :lagrange, (x, u) -> 0.5u^2)
 
 # replace default stop callbacks
 function mystop(i, sᵢ, dᵢ, xᵢ, gᵢ, fᵢ, ng₀, oTol, aTol, sTol, iterations)

examples/callbacks-stop.jl

@@ -16,8 +16,8 @@ A = [ 0 1
       0 0 ]
 B = [ 0
       1 ]
-constraint!(ocp, :dynamics, (x, u) -> A*x + B*u[1])
-objective!(ocp, :lagrange, (x, u) -> 0.5u[1]^2)
+constraint!(ocp, :dynamics, (x, u) -> A*x + B*u)
+objective!(ocp, :lagrange, (x, u) -> 0.5u^2)
 
 # initial iterate
 u_sol(t) = 6.0-12.0*t # solution
@@ -34,4 +34,4 @@ ps = plot(sol, size=(800, 400))
 # plot target
 point_style = (color=:black, seriestype=:scatter, markersize=3, markerstrokewidth=0, label="")
 plot!(ps[1], [tf], [xf[1]]; point_style...)
-plot!(ps[1], [tf], [xf[2]]; point_style...)
+plot!(ps[2], [tf], [xf[2]]; point_style...)

examples/direct-shooting.jl

@@ -1,10 +1,9 @@
 ### Goddard
 
 ## Direct solve
-
 using OptimalControl
 using Plots
-unicodeplots()
+#unicodeplots()
 
 # Parameters
 const Cd = 310
@@ -46,37 +45,36 @@ function F1(x)
     return F
 end
 
-f(x, u) = F0(x) + u[1]*F1(x)
+f(x, u) = F0(x) + u*F1(x)
 
 constraint!(ocp, :dynamics, f)
 
 # Solve
-sol = solve(ocp)
-plot(sol)
+direct_sol = solve(ocp)
+plot(direct_sol)
 
 ## Indirect solve
-
 using MINPACK
 
 # Bang controls
-u0(x, p) = [0.]
-u1(x, p) = [1.]
+u0(x, p) = 0.
+u1(x, p) = 1.
 
 # Computation of singular control of order 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)
-us(x, p) = [-H001(x, p) / H101(x, p)]
+us(x, p) = -H001(x, p) / H101(x, p)
 
 # Computation of boundary control
 remove_constraint!(ocp, :eq1)
 remove_constraint!(ocp, :eq3)
 constraint!(ocp, :boundary, (t0, x0, tf, xf) -> xf[3], mf, :eq4) # one value => equality (not boxed inequality); changed to equality constraint for shooting
 #
 g(x) = constraint(ocp, :eq2, :upper)(x, 0) # g(x, u) ≥ 0 (cf. nonnegative multiplier)
-ub(x, _) = [-Ad(F0, g)(x) / Ad(F1, g)(x)]
+ub(x, _) = -Ad(F0, g)(x) / Ad(F1, g)(x)
 μb(x, p) = H01(x, p) / Ad(F1, g)(x)
 
 f0 = Flow(ocp, u0)
@@ -101,20 +99,21 @@ function shoot!(s, p0, t1, t2, t3, tf) # B+ S C B0 structure
 end
 
 # Initialisation from direct solution
-t = sol.T; t = (t[1:end-1] + t[2:end]) / 2
-x = [ sol.X[i, 1:3] for i ∈ 1:N+1 ]; x = (x[1:end-1] + x[2:end]) / 2
-u = [ sol.U[i, 1  ] for i ∈ 1:N+1 ]; u = (u[1:end-1] + u[2:end]) / 2
-p = [ sol.P[i, 1:3] for i ∈ 1:N   ]
+t = direct_sol.times
+x = direct_sol.state
+u = direct_sol.control
+p = direct_sol.adjoint
+φ(t) = H1(x(t), p(t))
 
-u_plot = plot(t, u, xlabel = "t", ylabel = "u", legend = false)
-H1_plot = plot(t, H1.(x, p), xlabel = "t", ylabel = "H1", legend = false)
-g_plot = plot(t, g.(x), xlabel = "t", ylabel = "g", legend = false)
-display(plot(u_plot, H1_plot, g_plot, layout=(3,1), size=(800,600)))
+u_plot = plot(t, t -> u(t)[1], xlabel = "t", ylabel = "u", legend = false)
+H1_plot = plot(t, φ, xlabel = "t", ylabel = "H1", legend = false)
+g_plot = plot(t, g∘x, xlabel = "t", ylabel = "g", legend = false)
+display(plot(u_plot, H1_plot, g_plot, layout=(3,1), size=(400,300)))
 
 η = 1e-3
-t13 = t[ abs.(H1.(x, p)) .≤ η ]
-t23 = t[ 0 .≤ g.(x) .≤ η ]
-p0 = p[1]
+t13 = t[ abs.(φ.(t)) .≤ η ]
+t23 = t[ 0 .≤ (g∘x).(t) .≤ η ]
+p0 = p(t0)
 t1 = min(t13...)
 t2 = min(t23...)
 t3 = max(t23...)
@@ -125,22 +124,22 @@ println("Initial guess:\n", ξ)
 
 # Solve
 nle = (s, ξ) -> shoot!(s, ξ[1:3], ξ[4], ξ[5], ξ[6], ξ[7])
-sol = fsolve(nle, ξ, show_trace=true)
-println(sol)
+indirect_sol = fsolve(nle, ξ, show_trace=true)
+println(indirect_sol)
 
 # Plots
-p0 = sol.x[1:3]
-t1 = sol.x[4]
-t2 = sol.x[5]
-t3 = sol.x[6]
-tf = sol.x[7]
+p0 = indirect_sol.x[1:3]
+t1 = indirect_sol.x[4]
+t2 = indirect_sol.x[5]
+t3 = indirect_sol.x[6]
+tf = indirect_sol.x[7]
 
 f1sb0 = f1 * (t1, fs) * (t2, fb) * (t3, f0) # concatenation of the Hamiltonian flows
-sol = f1sb0((t0, tf), x0, p0)
-r_plot = plot(sol, idxs=(0, 1), xlabel="t", label="r")
-v_plot = plot(sol, idxs=(0, 2), xlabel="t", label="v")
-m_plot = plot(sol, idxs=(0, 3), xlabel="t", label="m")
-pr_plot = plot(sol, idxs=(0, 4), xlabel="t", label="p_r")
-pv_plot = plot(sol, idxs=(0, 5), xlabel="t", label="p_v")
-pm_plot = plot(sol, idxs=(0, 6), xlabel="t", label="p_m")
-plot(r_plot, pr_plot, v_plot, pv_plot, m_plot, pv_plot, layout=(3, 2))
+flow_sol = f1sb0((t0, tf), x0, p0)
+r_plot = plot(flow_sol, idxs=(0, 1), xlabel="t", label="r")
+v_plot = plot(flow_sol, idxs=(0, 2), xlabel="t", label="v")
+m_plot = plot(flow_sol, idxs=(0, 3), xlabel="t", label="m")
+pr_plot = plot(flow_sol, idxs=(0, 4), xlabel="t", label="p_r")
+pv_plot = plot(flow_sol, idxs=(0, 5), xlabel="t", label="p_v")
+pm_plot = plot(flow_sol, idxs=(0, 6), xlabel="t", label="p_m")
+plot(r_plot, pr_plot, v_plot, pv_plot, m_plot, pv_plot, layout=(3, 2), size=(600, 300))

examples/goddard.jl

@@ -2,6 +2,7 @@ using OptimalControl
 #
 using Printf
 using LinearAlgebra
+#
 
 # description of the optimal control problem
 t0 = 0
@@ -20,8 +21,8 @@ A = [ 0 1
       0 0 ]
 B = [ 0
       1 ]
-constraint!(ocp, :dynamics, (x, u) -> A*x + B*u[1])
-objective!(ocp, :lagrange, (x, u) -> 0.5u[1]^2)
+constraint!(ocp, :dynamics, (x, u) -> A*x + B*u)
+objective!(ocp, :lagrange, (x, u) -> 0.5u^2)
 
 # solution
 u_sol(t) = 6.0-12.0*t # solution
@@ -57,19 +58,19 @@ end
 # --------------------------------------------------------------------------------------------------
 # resolution
 
-T_default = OptimalControl.__grid(t0, tf)
+T_default = range(t0, tf, CTBase.__grid_size_direct_shooting())
 
 # init=nothing, grid=nothing
 sol = solve(ocp, :direct, :shooting, callbacks=(PrintCallback(my_cb_print(T_default)),))
 sol_1 = sol
-println(sol.T)
+println(time_steps(sol))
 
 # init=nothing, grid=T
 # the grid T may be non uniform
 N = 101
 T = range(t0, tf, N)
 sol = solve(ocp, :direct, :shooting, grid=T, callbacks=(PrintCallback(my_cb_print(T)),))
-println(sol.T)
+println(time_steps(sol))
 
 # init=U, grid=nothing
 # we assume U is given on a uniform grid
@@ -78,15 +79,15 @@ N = 101
 T = range(t0, tf, N)
 U = [[u_sol(T[i])-1.0] for i = 1:N-1]
 sol = solve(ocp, :direct, :shooting, init=U, callbacks=(PrintCallback(my_cb_print(T_default)),))
-println(sol.T)
+println(time_steps(sol))
 
 # init=U, grid=T
 # the grid T may be non uniform
 N = 101
 T = range(t0, tf, N)
 U = [[u_sol(T[i])-1.0] for i = 1:N-1]
 sol = solve(ocp, :direct, :shooting, init=U, grid=T, callbacks=(PrintCallback(my_cb_print(T)),))
-println(sol.T)
+println(time_steps(sol))
 
 # init=(T,U), grid=nothing
 # if U is not given on a uniform grid, you must give the grid
@@ -96,7 +97,7 @@ T = range(t0, tf, N)
 T = T.^2 # t0 = 0 and tf = 1 so it is ok
 U = [[u_sol(T[i])-1.0] for i = 1:N-1]
 sol = solve(ocp, :direct, :shooting, init=(T, U), callbacks=(PrintCallback(my_cb_print(T_default)),))
-println(sol.T)
+println(time_steps(sol))
 
 # init=(T,U), grid=T_
 # if U is not given on a uniform grid, you must give the grid
@@ -106,33 +107,33 @@ T = T.^2 # t0 = 0 and tf = 1 so it is ok
 U = [[u_sol(T[i])-1.0] for i = 1:N-1]
 T_ = range(t0, tf, 301)
 sol = solve(ocp, :direct, :shooting, init=(T, U), grid=T_, callbacks=(PrintCallback(my_cb_print(T_)),))
-println(sol.T)
+println(time_steps(sol))
 
 # init=t->u(t), grid=nothing
 u(t) = [u_sol(t)-1.0]
 sol = solve(ocp, :direct, :shooting, init=u, callbacks=(PrintCallback(my_cb_print(T_default)),))
-println(sol.T)
+println(time_steps(sol))
 
 # init=t->u(t), grid=T
 N = 101
 T = range(t0, tf, N)
 u(t) = [u_sol(t)-1.0]
 sol = solve(ocp, :direct, :shooting, init=u, grid=T, callbacks=(PrintCallback(my_cb_print(T)),))
-println(sol.T)
+println(time_steps(sol))
 
 # init=sol, grid=nothing
 sol = solve(ocp, :direct, :shooting, init=sol_1, callbacks=(PrintCallback(my_cb_print(T_default)),))
-println(sol.T)
+println(time_steps(sol))
 
 # init=sol, grid=T
 N = 101
 T = range(t0, tf, N)
 sol = solve(ocp, :direct, :shooting, init=sol_1, grid=T, callbacks=(PrintCallback(my_cb_print(T)),))
-println(sol.T)
+println(time_steps(sol))
 
 # make direct solve to init a direct shooting
 N = 201
 T = range(t0, tf, N)
 sol = solve(ocp, :direct, grid_size=N)
 sol = solve(ocp, :direct, :shooting, init=sol, grid=T, callbacks=(PrintCallback(my_cb_print(T)),))
-println(sol.T)
+println(time_steps(sol))

examples/goddard_f.ipynb

@@ -19,9 +19,13 @@ const __display = CTBase.__display
 include("flows.jl")
 include("solve.jl")
 
+# ----------------------------------------
+# to remove when put in the right package
+include("CTBase.jl")
+# ----------------------------------------
+
 # export functions only for user
 export solve
-export plot, plot!
 export Flow
 
 end

examples/initialization.jl

@@ -17,7 +17,7 @@ time!(ocp, :initial, t0) # if not provided, final time is free
 state!(ocp, 3) # state dim
 control!(ocp, 1) # control dim
 constraint!(ocp, :initial, x0)
-constraint!(ocp, :control, u -> u[1], 0., 1.)
+constraint!(ocp, :control, u -> u, 0., 1.)
 constraint!(ocp, :mixed, (x, u) -> x[1], r0, Inf, :state_con1)
 constraint!(ocp, :mixed, (x, u) -> x[2], 0., vmax, :state_con2)
 constraint!(ocp, :mixed, (x, u) -> x[3], m0, mf, :state_con3)
@@ -27,31 +27,31 @@ objective!(ocp, :mayer,  (t0, x0, tf, xf) -> xf[1], :max)
 D(x) = Cd * x[2]^2 * exp(-β*(x[1]-1))
 F0(x) = [ x[2], -D(x)/x[3]-1/x[1]^2, 0 ]
 F1(x) = [ 0, Tmax/x[3], -b*Tmax ]
-f(x, u) = F0(x) + u[1]*F1(x)
+f(x, u) = F0(x) + u*F1(x)
 constraint!(ocp, :dynamics, f)
 
 # --------------------------------------------------------
 # Indirect
 
 # Bang controls
-u0(x, p) = [0.]
-u1(x, p) = [1.]
+u0(x, p) = 0.
+u1(x, p) = 1.
 
 # Computation of singular control of order 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)
-us(x, p) = [-H001(x, p) / H101(x, p)]
+us(x, p) = -H001(x, p) / H101(x, p)
 
 # Computation of boundary control
 remove_constraint!(ocp, :state_con1)
 remove_constraint!(ocp, :state_con3)
 constraint!(ocp, :boundary, (t0, x0, tf, xf) -> xf[3], mf, :final_con) # one value => equality (not boxed inequality)
 
 g(x) = constraint(ocp, :state_con2, :upper)(x, 0.0) # g(x, u) ≥ 0 (cf. nonnegative multiplier)
-ub(x, _) = [-Ad(F0, g)(x) / Ad(F1, g)(x)]
+ub(x, _) = -Ad(F0, g)(x) / Ad(F1, g)(x)
 μb(x, p) = H01(x, p) / Ad(F1, g)(x)
 
 f0 = Flow(ocp, u0)