OptimalControl.jl

Appearance

Shooting ​

The payoff of everything else in this section: turn a flow into a root-finding problem for the unknown initial costate (and switching times, and free final time), and solve it.

julia
using OptimalControl
using OrdinaryDiffEqTsit5
using NonlinearSolve

The shooting equation ​

The PMP gives necessary conditions but not directly. Shooting turns "integrate the flow and check the boundary/transversality conditions" into a root-finding problem: guess the unknowns (initial costate, switching times, free final time), integrate, measure how far the result is from satisfying the conditions, and let a nonlinear solver close the gap.

Worked example throughout: minimise the final time for  ,   , from    to — bang-bang, one switch, free . The reference solution is  , one switch at  ,  .

julia
t0 = 0.0
x0 = [-1.0, 0.0]
xf = [0.0, 0.0]
u_max = 1.0
u_min = -1.0

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

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

The shooting function ​

The problem has 4 unknowns — (2), the switching time , and — and 4 residuals: the target state (2), the switching condition (  where the bang-bang control switches, since   there), and the free-time transversality    (the Mayer cost is  , normal case). Written in place, the way NonlinearSolve wants it:

julia
H(x, p, u) = p[1] * x[2] + p[2] * u - 1

function shoot!(s, ξ)
    p0 = ξ[1:2]
    t1, tfv = ξ[3], ξ[4]
    x1, p1 = f_max(t0, x0, p0, t1; variable=tfv)
    xf_, pf = f_min(t1, x1, p1, tfv; variable=tfv)
    s[1:2] = xf_ - xf
    s[3] = p1[2]
    s[4] = H(xf_, pf, u_min)
    return nothing
end

s = zeros(4)
shoot!(s, [1.0, 1.0, 1.0, 2.0])   # residual at the reference solution
sqrt(sum(abs2, s))
4.456631126570476e-16

Solving it ​

NonlinearSolve closes the gap from a perturbed guess:

julia
ξ_guess = [1.0, 1.0, 1.0, 2.0] .* 1.1
prob = NonlinearProblem((s, ξ, _) -> shoot!(s, ξ), ξ_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)

unsafe=true inside the loop ​

A nonlinear solver explores guesses that don't correspond to a real solution — some of them can make the flow's integration blow up. The default behaviour is to throw, which would abort the whole solve on the first bad guess:

julia
julia> f_blowup = Flow(VectorField(x -> x^2));  # ẋ = x², diverges before t=1

julia> f_blowup(0.0, 10.0, 1.0)
┌ Warning: Verbosity toggle: dt_epsilon 
│  At t=0.10000000044405617, dt was forced below floating point epsilon 1.3877787807814457e-17. Aborting. There is either an error in your model specification or the true solution is unstable (or it cannot be represented in Float64 precision).
│ 
│ Diagnostics:
│ 
│ State Analysis:
│   u[1] = 7.141e+15 has grown >1e6× its initial value
│ 
│ Error Analysis:
│   step error estimate EEst = 1.125 (a step is accepted when EEst <= 1)
│   largest contributors to EEst = internalnorm(atmp), where atmp is the tolerance-weighted local error per state component:
│     atmp[1] = 1.125, u[1] = 7.141e+15, uprev[1] = 6.497e+15
└ @ DiffEqBase /projects/ctb/gha-runner/julia_depot/packages/DiffEqBase/2YsuI/src/check_error.jl:37
SolverFailure → _check_retcode, CTSolversSciMLIntegrator.jl:500
│
│  ODE integration failed
│
│  Retcode  Unstable
│
│  Context  SciML solve
│  Hint     Try tightening tolerances (reltol, abstol) or changing the solver algorithm.
└─

unsafe=true returns whatever the integrator produced instead of throwing — garbage, but a value, letting the shooting residual carry the failure forward as "very wrong" rather than crashing the solve:

julia
f_blowup(0.0, 10.0, 1.0; unsafe=true)
7.141089928711383e15

Inside a shoot!, wrap the flow calls with unsafe=true so an intermediate failure shows up as a large residual for the solver to step away from, not an exception that stops the search.

Free final time ​

The transversality residual H(xf_, pf, u_min) = -1 above is the free-final-time condition — no separate machinery needed beyond adding it as a residual. For a smooth (non-switching) free-final-time problem, variable_costate=true gives the extra adjoint directly if the transversality condition is stated in terms of instead:

julia
xf_v, pf_v, pvf = f_max(
    t0, x0, [1.0, 1.0], 1.0; variable=2.0, variable_costate=true
)
pvf
0.0

Switching times as unknowns ​

The example above already carries one: is solved for alongside and . Each additional switch adds one more unknown time and one more switching-condition residual — see Multi-phase flows for concatenating the corresponding flows once the times are known (or being solved for), and Turnpike (bang–singular–bang) for a worked case with two unknown switching times bracketing a singular arc.

Getting a starting point from a direct solve ​

Manufacturing a shooting guess by hand doesn't scale — the standard workflow is to solve the same problem directly first, then read costate(sol)(t0) off as the initial guess:

julia
using NLPModelsIpopt
direct_sol = solve(ocp; display=false)
p0_guess = costate(direct_sol)(0.0)
p0_guess
2-element Vector{Float64}:
 0.9999999963782983
 0.9942579851682065

Checking against the direct solution ​

Rebuild the full bang–bang trajectory from the shooting solution — f_max, then f_min at the solved switch — and read its endpoint (a concatenated flow returns the stacked state–costate vector [q, v, p_q, p_v] at the final time):

julia
p0_sol, t1_sol, tf_sol = sol.u[1:2], sol.u[3], sol.u[4]

f_bb = f_max * (t1_sol, f_min)
zf = f_bb(t0, x0, p0_sol, tf_sol; variable=tf_sol)

zf[1:2]     # final state ≈ [0, 0] — the indirect solution hits the target
2-element Vector{Float64}:
 -8.891964458725137e-17
 -3.218871821797339e-17

The two methods agree on the optimum. The cost here is the final time, so comparing objectives is comparing :

julia
tf_sol, objective(direct_sol)
(2.0000000000000004, 1.9999999927565821)

See also ​