Constrained arcs

Along a boundary arc, a state constraint is active, and the flow needs more than a control law: the constraint and its multiplier. This page builds such flows, and assembles them with the unconstrained arcs, on two problems with known solutions.

using OptimalControl
using OrdinaryDiffEqTsit5
using NLPModelsIpopt
using Plots

The setting

With the conventions of Notation and conventions, a state constraint is written  , and along a boundary arc, where  , the pseudo-Hamiltonian gains the term  , with  :

On the arc, differentiating   gives the control, and the maximisation condition then gives the multiplier, both in feedback form. Flow(ocp, law; constraint=c, multiplier=μ) integrates the corresponding Hamiltonian flow. The functions have the signatures

c([t, ]x, u[, v])    # the constraint, c ≥ 0
u([t, ]x, p[, v])    # the control law on the arc
μ([t, ]x, p[, v])    # the multiplier on the arc

with t for a non-autonomous problem and v for a problem with a variable.

A constraint of order one

This example comes from the v2.0.4 manual, after Bonnard, Faubourg, Launay and Trélat^[1]. Minimise for  ,  , from   to  , under the constraint   :

t0, tf = 0, 2
x0, xf = 1, 1 / 2
lb = 0.1

ocp = @def begin
    t ∈ [t0, tf], time
    x ∈ R, state
    u ∈ R, control
    -1 ≤ u(t) ≤ 1
    x(t0) == x0
    x(tf) == xf
    x(t) - lb ≥ 0, (x_con)
    ẋ(t) == u(t)
    ∫(x(t)^2) → min
end

Here    and       . Along a boundary arc,   is constant, so  . The maximisation condition with   gives   on the arc, and differentiating,    : the multiplier is  . The solution has three arcs:   , the boundary arc, then   , with  ,   and   .

c(x, u) = x - lb         # the constraint, c ≥ 0
u_b(x, p) = 0            # the control on the boundary arc
μ(x, p) = 2x             # the multiplier on the boundary arc

f_minus = Flow(ocp, (x, p) -> -1)
f_boundary = Flow(ocp, u_b; constraint=c, multiplier=μ)
f_plus = Flow(ocp, (x, p) -> 1)

t1, t2 = 0.9, 1.6
f = f_minus * (t1, f_boundary) * (t2, f_plus)

p0 = -0.982237546583301
zf = f(t0, x0, p0, tf)

2-element Vector{Float64}:
 0.5000000000000001
 0.2477624534167015

The concatenation reaches   (its point call returns [x; p], see Multi-phase flows). The direct method finds the same boundary arc, , and the same cost:

sol_d = solve(ocp; grid_size=400, display=false)
sol_f = f((t0, tf), x0, p0)
objective(sol_d), objective(sol_f)

(0.3813306239985142, 0.381333893408902)

constraint and multiplier go together: one without the other is rejected.

julia> Flow(ocp, u_b; constraint=c)
IncorrectArgument → top-level scope, REPL[1]:2
│
│  `constraint` and `multiplier` must be given together
│
│  Got       only `constraint`
│  Expected  both `constraint` and `multiplier`, or neither
│
│  Context   Flow(ocp, law; constraint=…, multiplier=…) — pairing check
└─

Other ways to give the constraint

The constraint can also be a typed StateConstraint (or ControlConstraint, MixedConstraint, PathConstraint for constraints that involve the control), with the same result:

f_typed = Flow(ocp, u_b; constraint=StateConstraint(x -> x - lb), multiplier=μ)
f_typed(t1, lb, 0.0, t2), f_boundary(t1, lb, 0.0, t2)

((0.1, 0.0), (0.1, 0.0))

Several constraints active on the same arc are given as tuples, one multiplier per constraint:

f_two = Flow(ocp, u_b;
    constraint=(c, (x, u) -> 10 - x),     # a second constraint, x ≤ 10, inactive
    multiplier=(μ, (x, p) -> 0))
f_two(t1, lb, 0.0, t2)

(0.1, 0.0)

Finally, constraint=:x_con names a path constraint of the problem by its label. It uses the function stored in the model as it is, without its bound: the result is right only when the constraint is written as here, c(x) ≥ 0. Written x(t) ≥ lb or lb - x(t) ≤ 0, the flow on the arc is wrong, without an error (CTFlows#440). The examples of this documentation use the function form.

A costate jump: a constraint of order two

This example also comes from the v2.0.4 manual. It is the Bryson–Denham problem: minimise for  , from   to   , with   :

l = 1 / 9

ocp_bd = @def begin
    t ∈ [0, 1], time
    x ∈ R², state
    u ∈ R, control
    x(0) == [0, 1]
    x(1) == [0, -1]
    x₁(t) ≤ l, (x_con)
    ẋ(t) == [x₂(t), u(t)]
    0.5∫(u(t)^2) → min
end

Here     and       . Along a boundary arc,  , so   , and differentiating again,  : the constraint is of order two. The maximisation condition gives   on the arc, then     and   . Off the arc,  . The solution has three arcs, with   ,     and    , and the costate jumps at and : increases by . On a Hamiltonian flow, a vector jump acts on the costate (see Multi-phase flows):

f_free = Flow(ocp_bd, (x, p) -> p[2])
f_arc = Flow(ocp_bd, (x, p) -> 0; constraint=(x, u) -> l - x[1], multiplier=(x, p) -> 0)

ν = 18
f_bd = f_free * (3l, [ν, 0], f_arc) * (1 - 3l, [ν, 0], f_free)
zf_bd = f_bd(0, [0, 1], [-18, -6], 1)

4-element Vector{Float64}:
 -4.378891374111009e-16
 -0.9999999999999994
 18.0
 -5.999999999999999

The flow reaches the target  . Its cost is  , the value the direct method finds too:

sol_bd = solve(ocp_bd; display=false)
sol_bd_flow = f_bd((0, 1), [0, 1], [-18, -6])
objective(sol_bd), objective(sol_bd_flow)

(4.000142991483484, 4.000000525894766)

To plot the flow against the direct solution, the same flows are built with saveat, so that the trajectory has enough points (see Plotting a flow). Use such flows for trajectories only: their point call returns the value at the last saveat time, not at the final time (CTFlows#441).

fine = (saveat=range(0, 1, 101), dense=false)
f_bd_fine = Flow(ocp_bd, (x, p) -> p[2]; fine...) *
    (3l, [ν, 0], Flow(ocp_bd, (x, p) -> 0; constraint=(x, u) -> l - x[1], multiplier=(x, p) -> 0, fine...)) *
    (1 - 3l, [ν, 0], Flow(ocp_bd, (x, p) -> p[2]; fine...))
sol_bd_fine = f_bd_fine((0, 1), [0, 1], [-18, -6])

plt = plot(sol_bd; label="direct", size=(800, 600))
plot!(plt, sol_bd_fine; label="flow", linestyle=:dash)

Partial Hamiltonian

hamiltonian_type=:partial (see Total or partial Hamiltonian) applies to a constrained arc too. Along an optimal boundary arc, the law and the multiplier make both forms give the same flow:

f_partial = Flow(ocp, u_b; constraint=c, multiplier=μ, hamiltonian_type=:partial)
f_partial(t1, lb, 0.0, t2)

(0.1, 0.0)

Control-free problems

A problem without a control has no boundary arcs in this sense: its flow rejects constraint and multiplier.

julia> ocp_cf = @def begin
           t ∈ [0, 1], time
           x ∈ R, state
           x(0) == 1
           ẋ(t) == -x(t)
           ∫(x(t)^2) → min
       end;

julia> Flow(ocp_cf; constraint=(x, u) -> 1, multiplier=(x, p) -> 1)
PreconditionError → top-level scope, REPL[2]:2
│
│  constrained flows are not supported for control-free problems
│
│  Reason   a `constraint`/`multiplier` pair augments the pseudo-Hamiltonian H̃ + μ·g with the control; a control-free Flow(ocp) has no control law and no pseudo-Hamiltonian to carry the constraint term
│
│  Context  Flow(ocp; constraint=…, multiplier=…) — control-free constraint guard
│  Hint     use Flow(ocp, law; constraint=…, multiplier=…) with a control law
└─

See also

• Multi-phase flows: concatenations and jumps.

• Shooting: solve for the junction times and the initial costate.

• State constraint: a speed constraint on the double integrator, from the direct solution to the shooting.

---

1. B. Bonnard, L. Faubourg, G. Launay and E. Trélat, Optimal control with state constraints and the space shuttle re-entry problem, J. Dyn. Control Syst. 9 (2003), no. 2, 155–199. ↩︎