Skip to content

Flow(ocp) compatibility (control-free)

This page is a living compatibility reference for Flow(ocp) built from a control-free optimal control problem (ocp::CTModels.Models.Model with no control! declared): which state/costate types and call styles it accepts, each shown with a minimal, executable example. Every ✓ / ✗ in the tables below is demonstrated by a code block on this page and is re-run on every documentation build, so the page cannot drift from the code.

Scope: CPU, default AD backend (AutoForwardDiff). The tables are generated from probe/cpu/probe_cpu.jl (run it locally with julia --project=probe/cpu probe/cpu/probe_cpu.jl). GPU compatibility is a separate effort — see probe/gpu.

Unlike every other flow on this site, the state dimension of an OCP is fixed at construction (CTModels.Building.state!(pre, n)) — one Flow(ocp) does not accept every container size the way Flow(VectorField)/Flow(HamiltonianVectorField)/ Flow(Hamiltonian) do. This page therefore uses two control-free problems, built once: a scalar one (n=1) and a vector one (n=2), and every table below tests the scalar containers against the first and the vector-family containers 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 and SVector failed on the non-AD (basic) path — a fixed-size internal buffer, tracked as #358. Both are now fixed; the tables below reflect current (fixed) behaviour.

Flow(ocp) dispatches on the OCP's CTBase.Traits.ControlDependence trait: a control-free problem builds an OptimalControlFlow directly from the OCP's own dynamics and cost —      . With no Lagrange cost and dynamics   , this reduces to    , giving

julia
function build_ocp_free(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.dynamics!(pre, (r, t, x, u, v) -> (r .= -x; nothing))
    CTModels.Building.objective!(pre, :min; mayer=(x0, xf, v) -> sum(xf))
    return CTModels.Building.build(pre)
end

ocp1 = build_ocp_free(1)   # scalar
ocp2 = build_ocp_free(2)   # vector
f1 = Flows.Flow(ocp1; reltol=1e-8)
f2 = Flows.Flow(ocp2; reltol=1e-8)
OptimalControlFlow
├─ system: HamiltonianSystem
  ├─ time_dependence: Autonomous
  ├─ variable_dependence: Fixed
  ├─ Hamiltonian: 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.

Hamiltonian call — f(t0, x0, p0, tf)

The Hamiltonian point call f(t0, x0, p0, tf) returns (xf, pf); the trajectory call f((t0, tf), x0, p0) returns a CTModels.Solutions.Solution.

State/costate typepointtraj
Scalar Real
Vector Real
MVector Real
SVector Real
Matrix Real (batch)
MMatrix Real
SMatrix Real
Complex (any container)✗ (a)✗ (a)
ForwardDiff.Dual (any container)✗ (c)✗ (c)

Legend

  • — works; the result matches and is executed on this page.

  • — raises an error; see the corresponding note.

(a) Complex is not supported — same cause as Flow(Hamiltonian)

Flow(ocp) is AutoForwardDiff-backed internally, exactly like Flow(Hamiltonian): ArgumentError: Cannot create a dual over scalar type ComplexF64. See Flow(Hamiltonian)'s note (a) for the full explanation and the real/imaginary-splitting workaround — it applies verbatim here.

(c) A hand-built Dual collides with the flow's own internal AD

Same nested-AD tag collision as Flow(Hamiltonian): DualMismatchError. Differentiate the flow call with an outer ForwardDiff.jacobian/gradient instead of passing a hand-built Dual as x0/p0 — see that page's pattern, which applies unchanged here.

Examples

julia
julia> x0, p0 = 1.0, 0.5;

julia> f1(0.0, x0, p0, 1.0)   # (xf, pf) ≈ (x0·e⁻¹, p0·e)
(0.3678794412539543, 1.3591409143402517)
julia
julia> sol = f1((0.0, 1.0), x0, p0);

julia> CTModels.Components.state(sol)(1.0), CTModels.Components.costate(sol)(1.0)
(0.3678794412539543, 1.3591409143402515)

Vector, MVector, and SVector states (ocp2, n=2):

julia
julia> f2(0.0, [1.0, 2.0], [0.5, 0.3], 1.0)
([0.36787944124749616, 0.7357588824949923], [1.3591409143329156, 0.8154845485997495])

julia> f2(0.0, MVector{2}(1.0, 2.0), MVector{2}(0.5, 0.3), 1.0)
([0.36787944124749616, 0.7357588824949923], [1.3591409143329156, 0.8154845485997495])

julia> f2(0.0, SA[1.0, 2.0], SA[0.5, 0.3], 1.0)
([0.36787944124749616, 0.7357588824949923], [1.3591409143329156, 0.8154845485997495])

Matrix (batch) state — each column an independent trajectory, same convention as Flow(VectorField)/Flow(Hamiltonian):

julia
julia> f2(0.0, [1.0 2.0; 3.0 4.0], [0.5 0.3; 0.2 0.1], 1.0)
([0.3678794412394858 0.7357588824789716; 1.103638323718457 1.4715177649579432], [1.359140914324038 0.8154845485944227; 0.5436563657296151 0.27182818286480753])

Not supported: Complex, hand-built Dual

julia
julia> try
           f1(0.0, 1.0 + 2.0im, 0.5 + 0.1im, 1.0)
       catch e
           showerror(stdout, e)
       end
ArgumentError: Cannot create a dual over scalar type ComplexF64. If the type behaves as a scalar, define ForwardDiff.can_dual(::Type{ComplexF64}) = true.

Basic (state-only) call — f(t0, x0, tf)

For a control-free OCP, Flow(ocp) also exposes a state-only call with no costate (the direct-shooting use case, #230), dispatched by arity on the same f object. This path does not use automatic differentiation at all — the state equation is computed exactly — so its compatibility profile is structurally different from the Hamiltonian table above:

State typepointtraj
Scalar Real
Vector/MVector/SVector Real
Matrix/MMatrix/SMatrix Real
Scalar/Vector/MVector/SVector Complex
Matrix/SMatrix Complex
ForwardDiff.Dual scalar/Vector/MVector/SVector

Note the reversal from the Hamiltonian table: Complex and Dual work here (no AD on this path, exactly like Flow(VectorField)) — the Hamiltonian table's own Complex/Dual restrictions come from its internal AD backend, which this path never invokes. This path is fully green: every container tested works.

Examples

julia
julia> xf_basic = f1(0.0, x0, 1.0)
0.36787944125650673

julia> xf_basic  x0 * exp(-1.0)
true
julia
julia> sol_basic = f1((0.0, 1.0), x0);

julia> Trajectories.state(sol_basic)(1.0)
0.3678794412565068

Complex and Dual work directly — no AD on this path:

julia
julia> f1(0.0, 1.0 + 2.0im, 1.0)
0.36787944123358063 + 0.7357588824671613im

julia> f1(0.0, ForwardDiff.Dual(1.0, 1.0), 1.0)
Dual{Nothing}(0.36787944125650673,0.36787944125650673)

SVector and Matrix (batch) both work too — fixed as of 2026-07-24 (#358):

julia
julia> f2(0.0, SA[1.0, 2.0], 1.0)
2-element SVector{2, Float64} with indices SOneTo(2):
 0.3678794412375269
 0.7357588824750538

julia> f2(0.0, [1.0 2.0; 3.0 4.0], 1.0)
2×2 Matrix{Float64}:
 0.367879  0.735759
 1.10364   1.47152

variable= / variable_costate=

Flow(ocp) also accepts a NonFixed (variable-dependent) control-free problem. With    (autonomous, variable!(pre, 1)), the Hamiltonian is    , giving

and — since   is constant along the trajectory (the two exponentials cancel) — a closed form for the augmented variable-costate variable_costate=true (see Optimal control § Free times for what this mechanism is for):

julia
function build_ocp_free_variable(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.variable!(pre, 1)
    CTModels.Building.dynamics!(pre, (r, t, x, u, v) -> (r .= v[1] .* x; nothing))
    CTModels.Building.objective!(pre, :min; mayer=(x0, xf, v) -> sum(xf))
    return CTModels.Building.build(pre)
end

ocpv1 = build_ocp_free_variable(1)
ocpv2 = build_ocp_free_variable(2)
fv1 = Flows.Flow(ocpv1; reltol=1e-8)
fv2 = Flows.Flow(ocpv2; reltol=1e-8)
OptimalControlFlow
├─ system: HamiltonianSystem
  ├─ time_dependence: Autonomous
  ├─ variable_dependence: NonFixed
  ├─ Hamiltonian: autonomous, variable
    natural call: h(x, p, v)
    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/costate typevariable=v+ variable_costate=true
Scalar Real
Vector/SVector Real
Scalar Complex✗ (a)✗ (a)
Dual scalar✗ (c)✗ (c)

This axis was flagged as a possible gap before measuring it — it turns out not to be irregular: it fails exactly where the Hamiltonian table above fails (Complex/Dual, same AD cause), and works everywhere Real does, including variable_costate=true.

julia
julia> v = 0.7;

julia> xf, pf = fv1(0.0, x0, p0, 1.0; variable=v);

julia> [xf, pf]  [x0 * exp(v), p0 * exp(-v)]
true
julia
julia> xf, pf, pvf = fv1(0.0, x0, p0, 1.0; variable=v, variable_costate=true);

julia> pvf  -(x0 * p0)
true

See also