Singular control

For a control-affine system

the pseudo-Hamiltonian is   , where   are the Hamiltonian lifts of the vector fields and . The maximising control is given by the sign of the switching function . When vanishes on a time interval, the arc is singular: the maximisation does not give the control, which is computed by differentiating along the extremal.

This page computes a singular control by hand, then with the tools of Geometry, and checks it with the direct and the indirect methods. For a simpler case, where the singular control follows from the optimality conditions without brackets, see Turnpike (bang–singular–bang).

using OptimalControl
using NLPModelsIpopt
using OrdinaryDiffEqTsit5
using NonlinearSolve
using Plots

The problem

A vehicle moves in the plane, in a drift field. Its state is  , where is the position and the heading, and it turns at a bounded rate:

We look for the fastest transfer from the origin to the point , with free initial and final headings:

ocp = @def begin
    tf ∈ R, variable
    t ∈ [0, tf], time
    q = (x, y, θ) ∈ R³, state
    u ∈ R, control

    -1 ≤ u(t) ≤ 1
    -π / 2 ≤ θ(t) ≤ π / 2   # helps the direct method

    x(0) == 0
    y(0) == 0
    x(tf) == 1
    y(tf) == 0

    ∂(q)(t) == [cos(θ(t)), sin(θ(t)) + x(t), u(t)]

    tf → min
end

The bound on keeps the direct method away from trajectories that turn the long way round. It is not active at the solution, so the indirect method, which ignores it, finds the same solution. The system is control-affine, with

Direct method

direct_sol = solve(ocp; display=false)
variable(direct_sol)

1.1497309627876084

opt = (state_bounds_style=:none, control_bounds_style=:none)
plt = plot(direct_sol; label="direct", size=(800, 800), opt...)

The control stays strictly inside  : the arc is singular.

The singular control by hand

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

On a singular arc,  , and so are its time derivatives. Along an extremal, the time derivative of a function is , a Poisson bracket (see Time derivative along an extremal), with     .

First derivative. Since  ,

depends only on , so only the pair contributes:

On the singular arc,  , that is,  .

Second derivative.

and it vanishes on the arc, so wherever  ,

In  , only the pair contributes:     . In  , only the pair contributes:    . Therefore

The denominator does not vanish

If    were zero on the arc, then with  ,

This matrix has determinant , so   , and with  ,  . Then  , which is impossible: for a minimum-time problem,   along the extremal.

Simplification. Multiply the numerator and the denominator by , and use   in the denominator:

The singular control depends on the heading only. Along the direct solution, it matches the control:

tg = time_grid(direct_sol)
θ(t) = state(direct_sol)(t)[3]
plt_u = plot(direct_sol, :control; label="direct", size=(800, 300), opt...)
plot!(plt_u, tg, t -> sin(θ(t))^2; label="sin²θ", linestyle=:dash)

The singular control with brackets

The tools of Geometry compute the same brackets from the vector fields: their lifts, then the iterated Poisson brackets with @Lie:

F0(q) = [cos(q[3]), sin(q[3]) + q[1], 0]
F1(q) = [0, 0, 1]

H0 = Lift(F0)
H1 = Lift(F1)

H01 = @Lie {H0, H1}
H001 = @Lie {H0, H01}
H101 = @Lie {H1, H01}

us_bracket(q, p) = -H001(q, p) / H101(q, p)

At a point where   , that is,   and  , the formula gives :

q_s = [0.3, -0.1, 0.5]
p_s = [1.7, 1.7 * tan(q_s[3]), 0]
us_bracket(q_s, p_s), sin(q_s[3])^2

(0.22984884706593012, 0.22984884706593015)

Indirect method

Here the whole extremal is singular: the headings are free at both ends, so the transversality conditions give   , and   vanishes on the whole interval. The flow uses the singular control over the whole horizon:

u_s(q) = sin(q[3])^2
f = Flow(ocp, (q, p, tf) -> u_s(q))

There are five unknowns,  , the initial heading and the final time , and five conditions: the target (two equations), the transversality conditions of the headings, and the transversality condition of the final time,  , where   since  :

t0 = 0

function shoot!(s, ξ, _)
    p0, θ0, tf = ξ[1:3], ξ[4], ξ[5]
    q_tf, p_tf = f(t0, [0, 0, θ0], p0, tf; variable=tf)
    s[1] = q_tf[1] - 1               # x(tf) = 1
    s[2] = q_tf[2]                   # y(tf) = 0
    s[3] = p0[3]                     # free initial heading
    s[4] = p_tf[3]                   # free final heading
    s[5] = H0(q_tf, p_tf) - 1        # free final time
    return nothing
end

The direct solution gives the initial guess:

ξ_guess = [costate(direct_sol)(t0)..., state(direct_sol)(t0)[3], variable(direct_sol)]
prob = NonlinearProblem(shoot!, ξ_guess)
shooting_sol = NonlinearSolve.solve(prob; show_trace=Val(false))
p0_sol, θ0_sol, tf_sol = shooting_sol.u[1:3], shooting_sol.u[4], shooting_sol.u[5]

([0.7826328346326575, -0.6224836112086409, 0.0], -0.6719121189244032, 1.149730885756097)

The residual is zero:

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

5-element Vector{Float64}:
  0.0
 -2.396519136091636e-17
  0.0
 -2.197316037274865e-15
 -9.992007221626409e-16

Comparison

The two methods find the same final time, up to the discretisation error of the direct method:

indirect_sol = f((t0, tf_sol), [0, 0, θ0_sol], p0_sol; variable=tf_sol)
variable(direct_sol), tf_sol

(1.1497309627876084, 1.149730885756097)

For the plot, the flow is built with saveat, so that the trajectory has 101 points (see Plotting a flow):

f_plot = Flow(ocp, (q, p, tf) -> u_s(q); saveat=range(t0, tf_sol, 101), dense=false)
plot!(plt, f_plot((t0, tf_sol), [0, 0, θ0_sol], p0_sol; variable=tf_sol);
      label="indirect", linestyle=:dash, opt...)

Along the extremal, the switching function   stays zero, stays inside  , and the control computed with the brackets is :

qs, ps = state(indirect_sol), costate(indirect_sol)
ts = range(t0, tf_sol, 5)
[(ps(t)[3], us_bracket(qs(t), ps(t)) - u_s(qs(t))) for t in ts]

5-element Vector{Tuple{Float64, Float64}}:
 (0.0, 2.933525644621682e-11)
 (-1.691108885040827e-12, 2.024452827598111e-11)
 (2.3879003709128852e-11, 1.7184004219572557e-11)
 (1.8146338587845905e-11, -1.5047740831164447e-11)
 (-2.197317854924036e-15, -9.895695374240177e-12)

See also

• Poisson bracket: the brackets and the singular-control formula, in general.

• The @Lie macro: the bracket notation used above.

• Lift: the Hamiltonian lift of a vector field.

• Shooting: the shooting method.

• The Goddard tutorial: a problem with bang, singular and boundary arcs.