CTFlows.jl

Appearance

Constrained flows ​

Flow(ocp, law; constraint, multiplier) augments the pseudo-Hamiltonian with a path constraint term, for OCPs whose Pontryagin's Maximum Principle carries a state or mixed constraint alongside the control. This chapter covers the API, the sign convention, the hamiltonian_type semantics for constrained flows (:total vs :partial, and why they agree on a constrained arc), and two worked examples.

The augmented pseudo-Hamiltonian ​

Given the OCP's own pseudo-Hamiltonian       (see Optimal control), a path constraint   with multiplier augments it as

Pass constraint/multiplier to Flow(ocp, law; …) alongside a DynClosedLoop law — they must be given together:

julia
# OCP: ẋ = -x + u, min ∫ 0.5*u^2 dt  (autonomous, fixed, 1-D — scalar convention)
pre = CTModels.Building.PreModel()
CTModels.Building.time_dependence!(pre; autonomous=true)
CTModels.Building.time!(pre; t0=0.0, tf=1.0)
CTModels.Building.state!(pre, 1)
CTModels.Building.control!(pre, 1)
CTModels.Building.dynamics!(pre, (r, t, x, u, v) -> (r[1] = -x + u; nothing))
CTModels.Building.objective!(pre, :min; lagrange=(t, x, u, v) -> 0.5 * u^2)
ocp = CTModels.Building.build(pre)

law = Data.DynClosedLoop((x, p) -> p)   # u = p (from ∂H̃/∂u = p - u = 0)
c = 0.3
f = Flows.Flow(
    ocp, law;
    constraint=Data.StateConstraint(x -> x),
    multiplier=Data.Multiplier((x, p) -> c),
    reltol=1e-10,
)

A raw Function is accepted directly for constraint/multiplier (wrapped in a MixedConstraint/Multiplier with the OCP's own time/variable dependence, just like Flow(ocp, u::Function) — see Control laws):

julia
f2 = Flows.Flow(ocp, law; constraint=(x, u) -> x, multiplier=(x, p) -> c, reltol=1e-10)
xf, pf = f2(0.0, 1.0, 0.5, 1.0)
xf, pf
(0.7925558476344876, 0.8436563657430431)

A labelled :path constraint from the OCP itself also resolves — pass its Symbol label as constraint. The model stores it as   ; the flow normalises it as   (see the sign convention below) from its bounds:

Bounds stored in the modelconstraint=
lower only,  :label 
upper only,  :label 
equality,  :label 
two-sided,  :label => :lower 
two-sided,  :label => :upper 

So x(t) - lb ≥ 0, x(t) ≥ lb and lb - x(t) ≤ 0 give the same flow. A two-sided constraint has no active side known from the model: :label alone throws an IncorrectArgument, and the side is given as a Pair. For a vector constraint, the side is one Symbol for every component, or a tuple with one Symbol per component. (A Tuple of labels is a set of constraints, see Multiple constraints.)

The same constraint  , declared in the model as x ≥ 0.1 and as 0.1 - x ≤ 0, gives the flow of the function form   :

julia
function ocp_with(f!, lb, ub)
    pre = CTModels.Building.PreModel()
    CTModels.Building.time_dependence!(pre; autonomous=true)
    CTModels.Building.time!(pre; t0=0.0, tf=1.0)
    CTModels.Building.state!(pre, 1)
    CTModels.Building.control!(pre, 1)
    CTModels.Building.dynamics!(
        pre, (r, t, x, u, v) -> (r[1] = -x + u; nothing),
    )
    CTModels.Building.objective!(
        pre, :min; lagrange=(t, x, u, v) -> 0.5 * u^2,
    )
    CTModels.Building.constraint!(
        pre, :path; f=f!, lb=lb, ub=ub, label=:xcon,
    )
    return CTModels.Building.build(pre)
end

ocp_lower = ocp_with((r, t, x, u, v) -> (r[1] = x; nothing), 0.1, Inf)
ocp_upper = ocp_with((r, t, x, u, v) -> (r[1] = 0.1 - x; nothing), -Inf, 0.0)
μ(x, p) = c
f_lower = Flows.Flow(ocp_lower, law; constraint=:xcon, multiplier=μ)
f_upper = Flows.Flow(ocp_upper, law; constraint=:xcon, multiplier=μ)
f_fun = Flows.Flow(ocp, law; constraint=(x, u) -> x - 0.1, multiplier=μ)
f_lower(0.0, 1.0, 0.5, 1.0), f_upper(0.0, 1.0, 0.5, 1.0), f_fun(0.0, 1.0, 0.5, 1.0)
((0.7925558476786874, 0.8436563657634168), (0.7925558476786874, 0.8436563657634168), (0.7925558476786874, 0.8436563657634168))

Sign convention ​

Flow(ocp, law; constraint, multiplier) uses   with    . This is the opposite sign from Pontryagin's Maximum Principle as usually stated with a constraint   and multiplier   (Maurer's formulation): the two are related by

Under g ≥ 0, the multiplier is   on an active (boundary) arc, and   off it.

Multiple constraints ​

A single constraint carrier — a :path label or a raw function — can already be vector-valued (returning a vector, paired with a vector multiplier): the   pairing is dimension-agnostic, so this needs no special support.

To combine separately-defined constraints, pass tuples of matched length:

julia
g1(x, u) = x
g2(x, u) = 1.0 - x
μ1(x, p) = 0.4
μ2(x, p) = 0.1
ft = Flows.Flow(ocp, law; constraint=(g1, g2), multiplier=(μ1, μ2), reltol=1e-10)

Each gᵢ/μᵢ is resolved independently (raw functions, :path labels, or explicit carriers can be mixed within the tuple); the flow sums  . Both tuples must be given together, and both must have the same length — otherwise Flow throws a CTBase.Exceptions.IncorrectArgument.

hamiltonian_type for constrained flows ​

The hamiltonian_type keyword (:total default, :partial) works as for unconstrained flows (see Control laws), extended to the constraint term:

  • :total: compose with the law into a CTBase.Data.ComposedHamiltonian and differentiate through the law and the multiplier (total derivative in , , ).

  • :partial: freeze both the control   and the multiplier   at their feedback values, then differentiate   at those fixed values — is differentiated in , is not.

julia
ft_ = Flows.Flow(
    ocp, law; constraint=Data.StateConstraint(x -> x),
    multiplier=Data.Multiplier((x, p) -> c), hamiltonian_type=:total, reltol=1e-10,
)
fp_ = Flows.Flow(
    ocp, law; constraint=Data.StateConstraint(x -> x),
    multiplier=Data.Multiplier((x, p) -> c), hamiltonian_type=:partial, reltol=1e-10,
)
x0, p0, tf = 1.0, 0.5, 1.0
ft_(0.0, x0, p0, tf), fp_(0.0, x0, p0, tf)   # constant μ, stationary law ⇒ agree
((0.7925558476344876, 0.8436563657430431), (0.7925558476344876, 0.8436563657430431))

Why :total and :partial agree on a constrained arc ​

Pontryagin's Maximum Principle with a state constraint   (Maurer's formulation;    in this package's convention — see above) gives a pseudo-Hamiltonian      . On an active arc the true (closed-loop) Hamiltonian bakes in the optimal control and multiplier ( ); its symplectic gradient decomposes as

:partial is exactly — the genuine PMP field with the control and multiplier frozen at their feedback values — correct everywhere. :total adds the two extra chain-rule terms. Both vanish on the arc:

  •   vanishes because   on the arc (any constraint order).

  •   vanishes either because   (an invariant of the constrained dynamics — order 1: ; order 2: a richer invariant), or because   (order 2: a constant control on the arc).

So :total and :partial coincide on a constrained arc, and a shooting method's entry conditions put the trajectory exactly there — this is why both modes converge to the same solution in practice, as the two worked examples below demonstrate.

Plotting a constrained trajectory ​

A constrained Flow(ocp, law; constraint, multiplier) with a DynClosedLoop law is an OptimalControlFlow, so a trajectory call returns a CTModels.Solution with the control reconstructed from the law — plot it directly once Plots is loaded:

julia
sol = f((0.0, 1.0), 1.0, 0.5)   # trajectory call on the constrained flow f
plot(sol)

Worked example: Goddard problem (order 1) ​

The Goddard rocket problem (maximise final mass under a velocity path constraint, free final time) has four arcs — two bang arcs ( ,  ), a singular arc, and a velocity-boundary arc where the state constraint is active. On the boundary arc,   is a genuine invariant of the order-1 constrained dynamics (): both :total and :partial converge, from the known solution and from a perturbed Newton restart, to the same shooting solution (residual < 1e-7 in both modes). This is also an end-to-end integration test of the whole convenience surface — raw-function control law, raw-function constraint (MixedConstraint,  ), raw-function multiplier, and a NonFixed free final time — built directly from a CTModels.Model and rebuilt across all four arcs with Flow(ocp, law; constraint, multiplier), then reconstructed across phases with the * multi-phase operator (see Multi-phase flows) into a single CTModels.Solution. See test/suite/integration/test_goddard_ocp.jl in the repository.

Worked example: double integrator (order 2) ​

The time-minimal double integrator with a second-order position constraint   ( ) has a boundary arc where the constrained dynamics sit in a different corner of the theory: the boundary control is constant ( ) and the multiplier is identically zero ( ) on the arc — so H_1\,u_c' vanishes because   rather than because  , a distinct reason for the same :total ≡ :partial equivalence. The costate jumps at the two junctions between the interior and boundary arcs (f1 * (t1, jump, f2), see Multi-phase flows), and the reconstructed multi-phase solution's control is exactly on the boundary arc and non-zero on both interior arcs. Ported from the OptimalControl.jl state-constraint example. See test/cron/integration/test_double_integrator_state.jl (part of the weekly test suite).

See also ​