OptimalControl.jl

Appearance

Turnpike (bang–singular–bang) ​

A scalar system,   with   , is driven between two states over a fixed horizon while minimising . The dynamics are affine in and the cost does not depend on , so the pseudo-Hamiltonian is linear in : the optimal control is bang,   , except on an interval where the switching function vanishes, a singular arc. This is the smallest problem with one, and unlike Singular control, the singular control follows from the optimality conditions without Poisson brackets.

julia
using OptimalControl
using NLPModelsIpopt
using OrdinaryDiffEqTsit5
using NonlinearSolve
using Plots

The problem ​

The horizon is fixed,  , and the transfer goes from   to  . The state and the control are scalars, not vectors of length one:

julia
t0, tf = 0.0, 2.0
x0, xf = 1.0, 0.5

ocp = @def begin
    t ∈ [t0, tf], time
    x ∈ R, state
    u ∈ R, control
    -1 ≤ u(t) ≤ 1
    x(t0) == x0
    x(tf) == xf
    ẋ(t) == u(t)
    ∫(x(t)^2) → min
end

Direct method ​

julia
direct_sol = solve(ocp; grid_size=100, display=false)
plt = plot(direct_sol, :state, :control; label="direct", size=(800, 400))

The state goes down to the origin as fast as it can, stays there, then climbs to the target. The flat middle part, the turnpike, is the singular arc, where the control is , inside its bounds.

The singular control ​

With    (see Notation and conventions), the pseudo-Hamiltonian is

It is linear in , and maximised by  : the switching function is the costate . Where vanishes on an interval, and not only at isolated times, this rule does not give the control, which is singular. On such an interval,  ; the adjoint equation is    , so   there, and then   gives

The extremal is therefore bang–singular–bang:

ArcInterval
bangfrom to
singular
bangfrom to

The state decreases at unit rate from , so it reaches at  ; it increases at unit rate to , so it leaves the arc at  . On the first arc,     and  , so     and   . The cost is         .

Indirect method ​

Each arc has its flow, with a constant control:

julia
f_minus = Flow(ocp, (x, p) -> -1)
f_sing = Flow(ocp, (x, p) -> 0)
f_plus = Flow(ocp, (x, p) -> 1)

The unknowns are the initial costate and the two switching times; the horizon is fixed. Three conditions close the system: the trajectory enters the singular arc at  , with the switching function  , and it reaches the target at :

julia
function shoot!(s, ξ, _)
    p0, t1, t2 = ξ
    x1, p1 = f_minus(t0, x0, p0, t1)
    x2, p2 = f_sing(t1, x1, p1, t2)
    x3, _ = f_plus(t2, x2, p2, tf)
    s[1] = x1          # enter the singular arc at x = 0
    s[2] = p1          # where the switching function vanishes
    s[3] = x3 - xf     # reach the target
    return nothing
end

The direct solution gives the initial guess. The singular arc is where its switching function, the costate, stays close to zero, as in the Goddard tutorial:

julia
tg = time_grid(direct_sol)
p_d = costate(direct_sol)

t12 = tg[abs.(p_d.(tg)) .≤ 1e-3]   # the grid points on the singular arc
ξ_guess = [p_d(t0), minimum(t12), maximum(t12)]
3-element Vector{Float64}:
 -0.9801999456479715
  0.96
  1.52
julia
prob = NonlinearProblem(shoot!, ξ_guess)
shooting_sol = NonlinearSolve.solve(prob; show_trace=Val(false))
p0_sol, t1_sol, t2_sol = shooting_sol.u
3-element Vector{Float64}:
 -1.0000000000000036
  1.0000000000000002
  1.4999999999999998

The solver finds the values computed by hand,   ,   and  , with a zero residual:

julia
s = zeros(3)
shoot!(s, shooting_sol.u, nothing)
s
3-element Vector{Float64}:
 8.034495536742733e-17
 4.696458266519542e-16
 2.220446049250313e-16

Comparison ​

The concatenation of the three flows at the switching times gives the whole extremal, and its cost is , as is the cost of the direct solution, up to the discretisation error:

julia
φ = f_minus * (t1_sol, f_sing) * (t2_sol, f_plus)
indirect_sol = φ((t0, tf), x0, p0_sol)
objective(direct_sol), objective(indirect_sol)
(0.3749500107307793, 0.3749998537000679)

For the plot, the flows are built with saveat, so that the trajectory has 201 points (see Plotting a flow):

julia
fine = (saveat=range(t0, tf, 201), dense=false)
φ_plot = Flow(ocp, (x, p) -> -1; fine...) *
    (t1_sol, Flow(ocp, (x, p) -> 0; fine...)) *
    (t2_sol, Flow(ocp, (x, p) -> 1; fine...))
plot!(plt, φ_plot((t0, tf), x0, p0_sol), :state, :control; label="indirect", linestyle=:dash)

See also ​