Time minimisation (bang–bang)

The wagon of Energy minimisation again, with a bounded force this time,  : we transfer it from rest at    to rest at   as fast as possible. The dynamics are those of the double integrator,

with the limit conditions    and  , where   and the final time is free.

using OptimalControl
using NLPModelsIpopt
using OrdinaryDiffEqTsit5
using NonlinearSolve
using Plots

The problem

The final time is a variable of the problem, and the cost:

t0 = 0
x0, xf = [-1, 0], [0, 0]

ocp = @def begin
    tf ∈ R, variable
    t ∈ [t0, tf], time
    x = (q, v) ∈ R², state
    u ∈ R, control

    -1 ≤ u(t) ≤ 1

    x(t0) == x0
    x(tf) == xf

    ẋ(t) == [v(t), u(t)]

    tf → min
end

Direct method

We solve it with 200 time steps:

direct_sol = solve(ocp; grid_size=200, display=false)
variable(direct_sol)

1.9999999920052862

plt = plot(direct_sol; label="direct", size=(800, 600))

The control is , then : the wagon accelerates as much as it can, then brakes as much as it can.

The exact solution

With the conventions of Notation and conventions, the cost is a Mayer cost, so the pseudo-Hamiltonian is

where   is the costate. It is linear in , and maximised by  : the control is bang–bang, and is the switching function. Since   and   , is affine, so it changes sign at most once, at a time , from    to   .

On ,   and    . On , the wagon brakes, and stops at  , at    . The target gives   and  .

The switching condition   gives  . For a free final time with a Mayer cost, the transversality condition is  , that is,    , with  . Since     , we get  :

Indirect method

The two arcs have their flows, with the control and . The final time is a variable, so it is passed to the flows by keyword, variable=tf:

f_max = Flow(ocp, (x, p, v) -> 1)
f_min = Flow(ocp, (x, p, v) -> -1)

The shooting function has four unknowns, , and , and four conditions: the target, the switch, and the transversality condition:

H(x, p, u) = p[1] * x[2] + p[2] * u   # pseudo-Hamiltonian, Mayer form

function shoot!(s, ξ, _)
    p0, t1, tf = ξ[1:2], ξ[3], ξ[4]
    x1, p1 = f_max(t0, x0, p0, t1; variable=tf)
    x2, p2 = f_min(t1, x1, p1, tf; variable=tf)
    s[1:2] = x2 - xf             # target
    s[3] = p1[2]                 # switching function at t1
    s[4] = H(x2, p2, -1) - 1     # transversality, free final time
    return nothing
end

The direct solution gives the initial guess: its costate at , its final time, and the switching time, where its switching function changes sign:

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

p0_guess = p_d(t0)
t1_guess = tg[findfirst(t -> p_d(t)[2] < 0, tg)]
tf_guess = variable(direct_sol)

ξ_guess = [p0_guess..., t1_guess, tf_guess]

4-element Vector{Float64}:
 0.999999996002643
 0.993999298650952
 0.9999999960026431
 1.9999999920052862

We solve the shooting equations:

prob = NonlinearProblem(shoot!, ξ_guess)
shooting_sol = NonlinearSolve.solve(prob; show_trace=Val(false))
p0_sol, t1_sol, tf_sol = shooting_sol.u[1:2], shooting_sol.u[3], shooting_sol.u[4]

([1.0, 1.0], 1.0000000000000002, 2.0000000000000004)

These are the exact values,  ,   and  , and the residual is zero:

s = zeros(4)
shoot!(s, shooting_sol.u, nothing)
s

4-element Vector{Float64}:
 -8.891964458725137e-17
 -3.218871821797339e-17
 -4.0463145620773896e-17
  0.0

Comparison

The concatenation of the two flows at gives the whole extremal (see Multi-phase flows):

φ = f_max * (t1_sol, f_min)
indirect_sol = φ((t0, tf_sol), x0, p0_sol; variable=tf_sol)
state(indirect_sol)(tf_sol)

2-element Vector{Float64}:
 -8.891964458725137e-17
 -3.218871821797339e-17

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

fine = (saveat=range(t0, tf_sol, 201), dense=false)
φ_plot = Flow(ocp, (x, p, v) -> 1; fine...) * (t1_sol, Flow(ocp, (x, p, v) -> -1; fine...))
plot!(plt, φ_plot((t0, tf_sol), x0, p0_sol; variable=tf_sol); label="indirect", linestyle=:dash)

The two final times agree up to the discretisation error of the direct method:

variable(direct_sol), tf_sol

(1.9999999920052862, 2.0000000000000004)

Notes

• Shooting solves this problem step by step, starting from the multipliers of the direct solution.

• MINPACK.jl can replace NonlinearSolve to solve the shooting equations.

• The Goddard tutorial uses the same steps on a problem with bang, singular and boundary arcs.