OptimalControl.jl

Appearance

Writing a shooting function ​

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)

A simple shooting function ​

The simplest case has no switching — a single arc, fixed time. Out-of-place, the residual is just the target-state gap:

julia
function shoot_simple(p0)
    xf_, _ = f_max(t0, x0, p0, 1.0; variable=2.0)   # arbitrary tf here — illustration only
    return xf_ - xf
end

In-place, for the solver ​

The real 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):

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.195339464508916e-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 = 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.10000000044405623, 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.412e+15 has grown >1e6× its initial value
│ 
│ Error Analysis:
│   step error estimate EEst = 1.346 (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.346, u[1] = 7.412e+15, uprev[1] = 6.721e+15
└ @ DiffEqBase ~/.julia/packages/DiffEqBase/mLfsc/src/check_error.jl:20
SolverFailure → top-level scope, REPL[2]:2
│
│  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.412281287077651e15

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

Multiple shooting and switching times ​

The example above already has one: is an unknown 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).

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 ​

julia
objective(direct_sol)
1.9999999927565821
julia
indirect_sol = f_max((t0, sol.u[3]), x0, sol.u[1:2]; variable=sol.u[4])

Compare state(indirect_sol)/objective-derived quantities against the direct solve the same way Solution object already teaches — both describe the same optimal trajectory, found by two different methods.

See also ​