CTFlows.jl

Appearance

Flow(ocp, law) compatibility ​

This page is a living compatibility reference for Flow(ocp, law), built from an OCP with control and a control law: which state/costate types and call styles each of the three feedback kinds (DynClosedLoop, OpenLoop, ClosedLoop) accepts, each shown with a minimal, executable example. Every ✓ / ✗ is demonstrated by a code block re-run on every documentation build.

Scope: CPU, default AD backend. The tables are generated from probe/cpu/probe_cpu.jl. GPU compatibility is a separate effort — see probe/gpu.

As on Flow(ocp)'s page, the state dimension is fixed at OCP construction, so this page uses two problems built once — scalar (n=1) and vector (n=2) — and every table tests the scalar containers against the first, the vector-family ones against the second.

Last probed: 2026-07-24

This page's table shape — which containers and axes are tested — reflects probe/cpu/probe_cpu.jl as of this date. Every ✓/✗ cell below is still re-executed on every documentation build regardless (see Compatibility overview) — only the scope of what's tested can go stale, not the results shown. Until 2026-07-24, Matrix-family states failed on every path here and SVector failed on the non-AD (OpenLoop/ClosedLoop) path — the same fixed-size internal buffer as Flow(ocp), tracked as #358. Both are now fixed; the tables below reflect current (fixed) behaviour.

As documented in Control laws, the constructor dispatches on the law's feedback trait into two structurally different flow types:

Law typeFlow typeDynamicsInternal AD
DynClosedLoopOptimalControlFlow (Hamiltonian)state + costateyes (AutoForwardDiff)
OpenLoop / ClosedLoopControlledFlow (state flow)state onlyno

The OCP used on this page (the same one as Control laws):   ,  .

julia
function build_ocp_ctrl(n)
    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, n)
    CTModels.Building.control!(pre, n)
    CTModels.Building.dynamics!(pre, (r, t, x, u, v) -> (r .= u .- x; nothing))
    CTModels.Building.objective!(pre, :min; lagrange=(t, x, u, v) -> 0.5 * sum(abs2, u))
    return CTModels.Building.build(pre)
end

ocp1 = build_ocp_ctrl(1)
ocp2 = build_ocp_ctrl(2)
The (autonomous) optimal control problem is of the form:

    minimize  J(x, u) = ∫ f⁰(x(t), u(t)) dt, over [0.0, 1.0]

    subject to

        ẋ(t) = f(x(t), u(t)), t in [0.0, 1.0] a.e.,

    where x(t) ∈ R² and u(t) ∈ R².

DynClosedLoop — Hamiltonian flow, :total vs :partial ​

The law u = p is the stationary point of the pseudo-Hamiltonian     (   ), so the composed Hamiltonian is     , giving

julia
law = Data.DynClosedLoop((x, p) -> p)

f1_total = Flows.Flow(ocp1, law; reltol=1e-8)                                # :total (default)
f2_total = Flows.Flow(ocp2, law; reltol=1e-8)
f1_partial = Flows.Flow(ocp1, law; hamiltonian_type=:partial, reltol=1e-8)
f2_partial = Flows.Flow(ocp2, law; hamiltonian_type=:partial, reltol=1e-8)
OptimalControlFlow
├─ system: PseudoHamiltonianSystem
│  ├─ time_dependence: Autonomous
│  ├─ variable_dependence: Fixed
│  ├─ PseudoHamiltonian: autonomous, fixed (no variable)
│  │    natural call: h̃(x, p, u)
│  │    uniform call: h̃(t, x, p, u, v)
│  ├─ ControlLaw: dyn-closed-loop, autonomous, fixed (no variable)
│  │    natural call: u(x, p)
│  │    uniform call: u(t, x, p, v)
│  └─ backend: DifferentiationInterface{CPU}(ad_backend=AutoForwardDiff)
└─ integrator: SciML{CPU} (instance, id=:sciml)
   ├─ internalnorm = real_norm  [default]
   ├─ alg = Tsit5  [default]
   ├─ reltol = 1.0e-8  [user]
   ├─ save_everystep = auto  [default]
   ├─ abstol = 1.0e-8  [default]
   ├─ save_start = auto  [default]
   └─ dense = auto  [default]
   Tip: use describe(SciML{CPU}) to see all available options.

Because the feedback is stationary, :total (AD through the law) and :partial (AD at fixed u) agree exactly — the two columns below are a genuine consistency cross-check, not two independent measurements:

State type:total point:total traj:partial point:partial traj
Scalar Real✓✓✓✓
Vector Real✓✓✓✓
MVector Real✓✓✓✓
SVector Real✓✓✓✓
Matrix/MMatrix/SMatrix Real✓✓✓✓
Complex (any container)✗ (a)✗ (a)✗ (a)✗ (a)
ForwardDiff.Dual (any container)✗ (c)✗ (c)✗ (c)✗ (c)

Notes (a)/(c) are exactly Flow(ocp)'s notes (a)/(c): Complex fails the same AutoForwardDiff way, and a hand-built Dual collides with the flow's own AD the same way — neither is specific to having a control law. Matrix/MMatrix/SMatrix states work here too, fixed as of 2026-07-24 (#358).

julia
julia> x0, p0 = 1.0, 0.5;

julia> xf, pf = f1_total(0.0, x0, p0, 1.0);

julia> xf ≈ x0/exp(1.0) + (p0/2)*(exp(1.0) - 1/exp(1.0)), pf ≈ p0*exp(1.0)
(true, true)
julia
julia> f1_partial(0.0, x0, p0, 1.0)   # matches f1_total exactly — stationary feedback
(0.9554800381294778, 1.3591409143554003)
julia
julia> f2_total(0.0, SA[1.0, 2.0], SA[0.5, 0.3], 1.0)
([0.9554800381223815, 1.0883192406416302], [1.3591409143502615, 0.8154845486101568])

julia> f2_total(0.0, [1.0 2.0; 3.0 4.0], [0.5 0.3; 0.2 0.1], 1.0)
([0.9554800381100731 1.0883192406213191; 1.3386785625038267 1.589037884386333], [1.3591409143394941 0.8154845486036962; 0.5436563657357976 0.2718281828678988])

DynClosedLoop, constrained — does adding constraint=/multiplier= keep working? ​

Using a constant constraint=(x,u)->1.0 and multiplier=(x,p)->0.3: the product    is a pure constant added to , so it contributes nothing to or — the dynamics, and hence the analytic reference above, are unchanged. This isolates "does the constrained path still run" from "is the constrained physics right" (already covered by the constrained-flow tests), and directly answers the question this table exists to answer:

julia
f1c_total = Flows.Flow(
    ocp1, law; constraint=(x, u) -> 1.0, multiplier=(x, p) -> 0.3, reltol=1e-8
)
f2c_total = Flows.Flow(
    ocp2, law; constraint=(x, u) -> 1.0, multiplier=(x, p) -> 0.3, reltol=1e-8
)
OptimalControlFlow
├─ system: HamiltonianSystem
│  ├─ time_dependence: Autonomous
│  ├─ variable_dependence: Fixed
│  ├─ ComposedHamiltonian: autonomous, fixed (no variable)
│  │    natural call: h(x, p)
│  │    uniform call: h(t, x, p, v)
│  └─ backend: DifferentiationInterface{CPU}(ad_backend=AutoForwardDiff)
└─ integrator: SciML{CPU} (instance, id=:sciml)
   ├─ internalnorm = real_norm  [default]
   ├─ alg = Tsit5  [default]
   ├─ reltol = 1.0e-8  [user]
   ├─ save_everystep = auto  [default]
   ├─ abstol = 1.0e-8  [default]
   ├─ save_start = auto  [default]
   └─ dense = auto  [default]
   Tip: use describe(SciML{CPU}) to see all available options.
State type:total point:total traj:partial point:partial traj
Scalar Real✓✓✓✓
Vector/MVector/SVector Real✓✓✓✓
Matrix-family Real✓✓✓✓
Complex (any container)✗ (a)✗ (a)✗ (a)✗ (a)
ForwardDiff.Dual (any container)✗ (c)✗ (c)✗ (c)✗ (c)

Measured: the ✓/✗ footprint is identical to the unconstrained table above, cell for cell. Adding a path constraint does not change what state types the flow accepts.

julia
julia> f1c_total(0.0, x0, p0, 1.0) == f1_total(0.0, x0, p0, 1.0)
true

OpenLoop / ClosedLoop — state flow, no internal AD ​

Flow(ocp, ClosedLoop(x -> -x)) eliminates the control (  ), giving      . This path has no AD anywhere — the control law and the dynamics are both evaluated directly:

julia
law_cl = Data.ClosedLoop(x -> -x)
fcl1 = Flows.Flow(ocp1, law_cl; reltol=1e-8)
fcl2 = Flows.Flow(ocp2, law_cl; reltol=1e-8)
ControlledFlow
├─ system: VectorFieldSystem
│  ├─ time_dependence: Autonomous
│  ├─ variable_dependence: Fixed
│  └─ ComposedVectorField: autonomous, fixed (no variable), out-of-place
│       natural call: f(x)
│       uniform call: f(t, x, v)
└─ integrator: SciML{CPU} (instance, id=:sciml)
   ├─ internalnorm = real_norm  [default]
   ├─ alg = Tsit5  [default]
   ├─ reltol = 1.0e-8  [user]
   ├─ save_everystep = auto  [default]
   ├─ abstol = 1.0e-8  [default]
   ├─ save_start = auto  [default]
   └─ dense = auto  [default]
   Tip: use describe(SciML{CPU}) to see all available options.
State typepointtraj
Scalar Real✓✓
Vector/MVector/SVector Real✓✓
Matrix-family Real✓✓
Scalar/Vector/MVector/SVector Complex✓✓
Matrix-family Complex✓✓
ForwardDiff.Dual scalar/Vector/MVector/SVector✓✗ (e)

Complex and Dual work here (no AD in the RHS), like Flow(ocp)'s basic call — SVector/Matrix-family now work too, fixed as of 2026-07-24 (#358) — one asymmetry remains:

(e) Dual works at a point, but not in a trajectory

ForwardDiff.Dual states integrate correctly for the point call, but the trajectory call raises MethodError: no method matching Float64(::ForwardDiff.Dual{...}). The trajectory path additionally computes the OCP objective (Mayer + Lagrange) from the reconstructed control and state, and that computation is not Dual-transparent — unlike the point call, which only integrates the state. If you need a trajectory with a Dual state, this is the concrete blocker to work around.

julia
julia> xf = fcl1(0.0, x0, 1.0)
0.13533528340900094

julia> xf ≈ x0 * exp(-2.0)
true
julia
julia> fcl1(0.0, ForwardDiff.Dual(1.0, 1.0), 1.0)   # point call: works
Dual{Nothing}(0.13533528340900094,0.13533528340900094)
julia
julia> try
           fcl1((0.0, 1.0), ForwardDiff.Dual(1.0, 1.0))   # trajectory call: fails (note e)
       catch e
           showerror(stdout, e)
       end
MethodError: no method matching Float64(::ForwardDiff.Dual{Nothing, Float64, 1})
The type `Float64` exists, but no method is defined for this combination of argument types when trying to construct it.

Closest candidates are:
  (::Type{T})(::Real, !Matched::RoundingMode) where T<:AbstractFloat
   @ Base rounding.jl:265
  (::Type{T})(::T) where T<:Number
   @ Core boot.jl:965
  Float64(!Matched::IrrationalConstants.Fourinvπ)
   @ IrrationalConstants ~/.julia/packages/IrrationalConstants/RokwY/src/macro.jl:111
  ...
julia
julia> sol = fcl2((0.0, 1.0), [1.0, 2.0]);

julia> Trajectories.state(sol)(1.0), Trajectories.control(sol)(1.0)
([0.1353352833695603, 0.2706705667391206], [-0.1353352833695603, -0.2706705667391206])

SVector and Matrix (batch) both work, point and trajectory:

julia
julia> fcl2(0.0, SA[1.0, 2.0], 1.0)
2-element SVector{2, Float64} with indices SOneTo(2):
 0.1353352833695603
 0.2706705667391206

julia> fcl2(0.0, [1.0 2.0; 3.0 4.0], 1.0)
2×2 Matrix{Float64}:
 0.135335  0.270671
 0.406006  0.541341

See also ​