OptimalControl.jl

Appearance

Shooting ​

The maximum principle gives the extremals, but not the initial costate, the switching times or a free final time. Shooting finds them: write the conditions the extremal must satisfy as equations in these unknowns, integrate the flows to evaluate them, and solve the equations with a nonlinear solver. The direct method gives the starting point.

julia
using OptimalControl
using OrdinaryDiffEqTsit5
using NonlinearSolve
using NLPModelsIpopt

The problem ​

We take the minimum-time transfer of the double integrator: minimise for  ,   , from    to :

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

ocp = @def begin
    tf ∈ R, variable
    t ∈ [0, tf], time
    x = (q, v) ∈ R², state
    u ∈ R, control
    -1 ≤ u(t) ≤ 1
    x(0) == x0, (start)
    x(tf) == xf
    ẋ(t) == [v(t), u(t)]
    tf → min
end

With    and the Mayer cost , the pseudo-Hamiltonian is   , maximised by  : the switching function is . As on Multi-phase flows, the solution is    then   , with  ,   and  . Shooting will find these values.

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

The shooting function ​

There are four unknowns,  , and , and four conditions:

  • the target,   (two equations);

  • the switch: the switching function vanishes at ,  ;

  • the transversality condition of the free final time: for the Mayer cost , the maximised Hamiltonian equals   at , that is    (see Notation and conventions).

The shooting function integrates the two arcs and returns the four residuals. It is written in place, as NonlinearSolve expects:

julia
H(x, p, u) = p[1] * x[2] + p[2] * u   # p ⋅ f(x, u)

function shoot!(s, ξ, _)
    p0, t1, tf = ξ[1:2], ξ[3], ξ[4]
    x1, p1 = f_max(t0, x0, p0, t1; variable=tf, unsafe=true)
    x2, p2 = f_min(t1, x1, p1, tf; variable=tf, unsafe=true)
    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

At the solution, the residuals vanish:

julia
s = zeros(4)
shoot!(s, [1, 1, 1, 2], nothing)
s
4-element Vector{Float64}:
 -2.313802281140895e-16
 -3.425505599574657e-17
  1.8158145930425736e-16
 -3.3306690738754696e-16

unsafe=true ​

While it searches, a nonlinear solver tries unknowns that are far from the solution, and the integration can fail on them. By default a flow then raises an error, which would stop the search. With unsafe=true, the flow returns what the integrator computed instead, and the residual is just large. On  , which blows up before   from  :

julia
julia> f_blowup = Flow(VectorField(x -> x^2));

julia> f_blowup(0, 10.0, 1)
SolverFailure → #67, REPL[2]:3
│
│  ODE integration failed
│
│  Retcode  Unstable
│
│  Context  SciML solve
│  Hint     Try tightening tolerances (reltol, abstol) or changing the solver algorithm.
└─
julia
f_blowup(0, 10.0, 1; unsafe=true)
7.141089928711383e15

This is why shoot! above calls the flows with unsafe=true.

Starting from the direct solution ​

A shooting method converges from a good initial guess. The direct method provides one:

  • the initial costate is the multiplier of the initial condition, dual(sol, ocp, :start) (more accurate than costate(sol)(0), see Solution object);

  • the switching time is read off the control, where it changes sign;

  • the final time is the variable.

julia
sol_d = solve(ocp; display=false)

p0_guess = dual(sol_d, ocp, :start)
tg = time_grid(sol_d)
t1_guess = tg[findfirst(t -> control(sol_d)(t) < 0, tg)]
tf_guess = variable(sol_d)

ξ_guess = [p0_guess; t1_guess; tf_guess]
4-element Vector{Float64}:
 0.9999999963782983
 0.9982579851391655
 0.9999999963782911
 1.9999999927565821

Solving the shooting equations ​

julia
prob = NonlinearProblem(shoot!, ξ_guess)
sol = NonlinearSolve.solve(prob, SimpleNewtonRaphson(); abstol=1e-10, reltol=1e-10)
sol.u, sol.retcode
([1.0, 1.0, 1.0000000000000002, 2.0000000000000004], SciMLBase.ReturnCode.Success)

The solver finds  ,   and  .

Checking the solution ​

The two arcs, concatenated at the switching time found, reach the target. The point call of a concatenation returns [x; p] (see Multi-phase flows):

julia
p0_s, t1_s, tf_s = sol.u[1:2], sol.u[3], sol.u[4]

f_bb = f_max * (t1_s, f_min)
zf = f_bb(t0, x0, p0_s, tf_s; variable=tf_s)
zf[1:2]
2-element Vector{Float64}:
 -8.891964458725137e-17
 -3.218871821797339e-17

The final time is the cost, so the indirect and the direct solutions agree when their final times do:

julia
tf_s, objective(sol_d)
(2.0000000000000004, 1.9999999927565821)

More switches, other unknowns ​

Each additional switch adds an unknown time and a switching condition; a boundary arc adds its entry and exit times, with the junction conditions (see Constrained arcs). When the variable enters the dynamics or the cost, its transversality condition uses the costate of the variable, computed with variable_costate=true (see From an OCP). The examples solve several such problems, for instance Turnpike (bang–singular–bang), with two switching times around a singular arc.

See also ​