Roadmap

[ocots]

Roadmap v4

Last updated: July 7, 2026 Supersedes: roadmap-v3.md

Context

Between the v3 roadmap (July 3) and now, a burst of work (PRs #303, #304, #307, #310 and the "1-D = scalar" refactor) landed most of v3 §1, §2, §6 and the §7 tests. This v4 records what actually shipped — often in a different shape than v3 predicted — and re-scopes the remaining work with firmer design decisions.

Key differences from what v3 anticipated:

• The control-law wrapper types (OpenLoop, ClosedLoop, DynClosedLoop) live in CTBase.Data (a single parametric ControlLaw{F,FB,TD,VD} with feedback traits), not in CTFlows as three separate structs. Dispatch is on the CTBase.Traits.feedback trait.
• There is no arity-based inference of the wrapper type. Flow(ocp, u::Function) always builds a DynClosedLoop; OpenLoop/ClosedLoop must be constructed explicitly.
• The pseudo-Hamiltonian path is real: PseudoHamiltonianSystem, ComposedHamiltonian, ComposedVectorField, ControlledVectorField (all in CTBase.Data), with a hamiltonian_type = :total | :partial action option.
• The variable_costate free-time semantics were decided: the flow always integrates the naive augmented adjoint ṗv = -∂H/∂v with pv(t0) = 0; the free-time transversality condition is a mitigated residual written in the shooting method, not computed by the flow.

---

Completed since v3 (for the record)

v3 section	Status	Notes
§1 OpenLoop / ClosedLoop / DynClosedLoop	✅ done	Flow(ocp, law), Flow(h̃, law), Flow(fc, law); :total / :partial. PRs #304, #307
§1 control-free Flow(ocp)	✅ done	PR #263
§2 Hamiltonian / pseudo-Hamiltonian / gradient getters	✅ done	hamiltonian, pseudo_hamiltonian, control_law, *_gradient functors. PR #310
§6 scalar / vector ("1-D = scalar")	✅ done	_asvec dropped; _coerce_state / _dim_coerce centralise coercion at the CTModels boundary
§7 variable-costate free-time tests	✅ partial	tests for v = tf; semantics decided (see §7 below). v = t0 and v = (t0, tf) to confirm
v2 carry-over: closures → callable structs	✅ done	#245, #256
v2 carry-over: traits → CTBase	✅ done	#246
v2 carry-over: differential geometry → CTLie	✅ done	#248

The remaining sections below are the v4 work items.

---

1. Implicit control (IFT / DAE) — ∂H̃/∂u = 0

Issue: #46. Status: open, not started. Reference implementation exists in the old gallery.

Goal: When no control law is provided, compute the control automatically from the pseudo-Hamiltonian stationarity condition ∂H̃/∂u = 0, and integrate the resulting flow. This is the missing branch of the with-control OCP story: we now have Flow(ocp, law) (control eliminated by an explicit feedback), but not Flow(ocp) with control implicit.

Two strategies

• IFT (Implicit Function Theorem) — Differentiate ∂H̃/∂u = 0 along the flow to get an explicit ODE for u:

u̇ = -(∂²H̃/∂u²)⁻¹ · (∂²H̃/∂z∂u)ᵀ · ẋ_H, with z = (x, p), ẋ_H = (∂H̃/∂p, -∂H̃/∂x).

The augmented state (x, p, u) is integrated as an ODE. Requires an initial control guess u(t0). The old-gallery RHS computes ∂h/∂z, ∂h/∂u, ∂²h/∂u², ∂²h/∂z∂u by AD, sets dw[1:2n] = hv and dw[2n+1:end] = -(∂²h/∂u²)⁻¹(∂²h/∂z∂u)ᵀ·hv.

• DAE — Keep ∂H̃/∂u = 0 as an algebraic constraint and integrate (x, p, u) as a DAE (DAEProblem + IDA/Sundials or a mass-matrix / DFBDF formulation).

Design decision to make: option vs. strategy

The choice of implicit-control method (and IFT-vs-DAE within it) can be surfaced either as:

• a keyword option on the flow call / constructor (e.g. control = :implicit, implicit_method = :ift | :dae), routed through the existing action-option machinery (like hamiltonian_type), or
• a CTBase.Strategies strategy — a first-class strategy family the same way the backend (:di) and integrator (:sciml) strategies are, so the implicit-control solver is a pluggable, dispatchable component.

➡️ Decision pending. Lean toward a strategy if the IFT/DAE solver needs its own options and extension boundary (e.g. DAE pulls in Sundials as a weak dep); lean toward a simple option if it is just a branch in the RHS builder. Prototype IFT first (pure ODE, no new heavy dep), then decide whether DAE justifies the strategy machinery.

Interaction with the pseudo-Hamiltonian infrastructure

We already have PseudoHamiltonianSystem (AD at fixed u) and the pseudo_*_gradient getters. The IFT RHS needs the Hessian ∂²H̃/∂u² and the mixed Jacobian ∂²H̃/∂z∂u — new AD primitives to add to CTBase.Differentiation, alongside the existing pseudo_hamiltonian_gradient / pseudo_variable_gradient.

Control passed as a keyword (not positional)

For implicit-control flows the initial control u0 is a keyword argument, consistent with variable / variable_costate:

xf, pf, uf = flow(t0, x0, p0, tf; control = u0, variable = v)   # point eval
sol        = flow((t0, tf), x0, p0; control = u0, variable = v)  # trajectory

control is required for implicit-control flows and rejected otherwise (trait-based validation, same pattern as variable for NonFixed/Fixed).

Open questions

• IFT Hessian via AD (ForwardDiff) or user-supplied? Default AD.
• Which DAE solvers to support, and behind which extension?
• Does the augmented (x,p,u) state interact with the existing variable_costate augmentation (i.e. (x,p,u,pv))? Order and coercion of the augmented block.

---

2. Flow on OCP + control + constraints

Status: OCP + control ✅ done; + constraints not started (blocked on §3 design).

Goal: Extend Flow(ocp, law) (and the implicit-control flow of §1) to OCPs that carry path constraints — state g(t,x) ≤ 0, control g(t,u) ≤ 0, mixed g(t,x,u) ≤ 0. The flow must integrate the constrained Hamiltonian/state system.

This is deliberately split from §3 (dual variables): here the concern is only that the constrained dynamics are integrated correctly. Whether/how multipliers are exposed is the subject of §3 and must be settled first, because it dictates the augmented-state layout.

Concretely this needs:

• constraint carriers (function g + its kind + multiplier dimension) — likely in CTBase.Data next to ControlLaw, for symmetry;
• a way for the flow to read the OCP's declared constraints (coordinate with CTModels);
• the constrained RHS (junction/activation handling for state constraints is the hard part — order of the constraint, boundary/tangency conditions).

➡️ Do not start implementation until §3's storage question is decided.

---

3. Constraints & multipliers — design/reflection task (do not return duals)

Issue: #103 (closed as "captured in roadmap"). Status: reopen as a design/reflection item. This reverses the v3 §3 plan.

Revised position (v4): returning the dual variables from the flow is a bad idea and is dropped as the default deliverable. The real problem is not "how to return duals" but how the CTModels.Solutions.Solution type should represent the absence/presence of multipliers. That needs thought before any code.

Why returning duals cleanly is hard

• A problem can have several constraints, but at a given time only some are active. The multiplier only "exists" for the active ones.
• Padding the inactive ones with zeros is possible, but we do not have a clean link from "this stored multiplier" back to "which OCP constraint it belongs to and whether it is active" — so the stored object is ambiguous and hard to interpret downstream.
• Accessors like dual(sol, ...) therefore cannot return something well-defined in the general (multi-constraint, partially-active) case.
• Display/printing of a Solution must degrade gracefully when duals are absent.

The actual questions to resolve

1. Presence/absence of duals in Solution: trait or type?

• A trait on Solution (HasDuals / NoDuals) keeps one type but branches every accessor/printer on the trait.
• A separate type (e.g. ConstrainedSolution <: AbstractSolution) keeps the dual-free Solution clean and pushes dual handling to a subtype.
• Coordinate with CTModels, since Solution lives there.

2. Accessor contract. If dual(sol, ...) cannot be well-defined, what is the supported query? Per-named-constraint lookup? Return nothing / throw for inactive?
3. Active-set bookkeeping. Is the flow even the right place to know the active set, or is that a solver-side concept? Maybe the flow only integrates the constrained dynamics (§2) and multipliers are never a flow output at all.
4. Display. How a Solution prints with / without duals, with / without an active set.

Deliverable for v4

A written design note answering (1)–(4), agreed with CTModels, before any constraint-multiplier code. Implementation (if any) is scoped afterward. The constrained dynamics (§2) can proceed independently of duals once the storage question is closed.

---

4. Flow derivative via dual numbers (variational equations / IND)

Issue: #93 (refs CTBase #25). Status: open, not started. Confirmed worth doing.

Goal: Differentiate through the flow correctly. When the flow is evaluated on ForwardDiff.Dual inputs, integrate the variational equations (Internal Numerical Differentiation, IND) instead of naively pushing duals through the black-box integrator, so that ∂φ/∂x₀, ∂φ/∂p₀, ∂φ/∂t₀, ∂φ/∂t_f are obtained from the sensitivity ODE.

Design

• Code the variational system. Augment (x, p) with (δx, δp) propagated by the Jacobian of the vector field: δż = J·δz. For Hamiltonian flows this is the Hessian of H (Jacobian of X_H). This is the IND part.
• Dispatch on dual numbers + a tag. Overload the flow call for Dual-typed inputs and use a dedicated tag so we can recognise "our" duals and route to the variational path (rather than colliding with an outer AD pass). The tag disambiguates nested/foreign differentiation.
• Activation: option vs. strategy — same open choice as §1. Either a keyword (variational = true / derivative = true) or a CTBase.Strategies strategy that selects the derivative-propagation method (IND/variational vs. plain dual push-forward vs. future adjoint/reverse). ➡️ Decision pending, ideally made together with §1 so implicit-control and derivative both use the same option-or-strategy convention.

Implementation steps

1. Variational RHS functor augmenting the system with J·δz (Hessian of H via AD).
2. Tagged Dual dispatch on the flow call → variational path.
3. Extract derivatives from the dual components of the output.
4. Cross-check against finite differences and against a plain (black-box) dual push-forward.

---

5. GPU support & testing

Issues: #249, #111. Status: open, not started. Testing on GPU is called out as important.

Goal: Flows run on GPU arrays (CuArray, and ideally AMDGPU/Metal) with GPU-compatible integrators, and we have actual GPU tests in CI or a guarded test env.

Steps

1. Audit array allocations in Systems, Configs, Trajectories, and the RHS functors — infer the array type from the input instead of hard-coding Vector{Float64}.
2. Ensure Integrators.build_* preserves the input array type through ODEProblem.
3. Add a CPU / GPU strategy (singleton types) — this dovetails with the option-vs-strategy discussions in §1/§4; GPU is a natural strategy candidate.
4. GPU tests (guarded by a CUDA/extension test environment) — the important part.
5. Document GPU usage.

---

6. Non-Hamiltonian flow without a costate (Flow(ocp) for direct shooting)

Issue: #230. Status: open; the Hamiltonian branch exists, the basic branch does not.

Goal: For an OCP without control but with variables, support a basic (non-Hamiltonian) state flow where the user does not pass a costate p0:

xf = f(t0, x0, tf; variable = v)          # basic flow, no p0 — for direct shooting

in addition to the existing Hamiltonian flow

xf, pf = f(t0, x0, p0, tf; variable = v)  # Hamiltonian flow (already works)

This is the direct-shooting use case (parameter estimation / optimal design): integrate ẋ = f(t, x, ∅, v) as a plain ODE with variables, no adjoint.

Note: the v3-era augment=true keyword mentioned in issue #230 no longer exists — it has been replaced by variable_costate. So issue #230's item 1 (Hamiltonian + augment) maps to today's variable_costate; only item 2 (basic non-Hamiltonian flow with variables, no p0) is the actual remaining work here.

Steps

1. Detect the "no p0, control-free OCP" call shape and route to a VectorFieldSystem built from the OCP dynamics f(t,x,∅,v) (not a Hamiltonian system).
2. Enforce the compulsory variable keyword for NonFixed (same rule as elsewhere).
3. Return the final state (and a state trajectory on a tspan call) — no costate, no Solution-with-costate. Decide the trajectory return type (bare state trajectory vs. a trimmed Solution).

---

7. Getter naming & vector-field getter

Issue: #185. Status: partially addressed by PR #310; naming reconciliation pending.

We shipped hamiltonian, pseudo_hamiltonian, control_law, and the *_gradient functors (§2 above). Issue #185 asks specifically for, from a flow:

f  = Flow(ocp, u)
h  = hamiltonian(f)              # ✅ have it (delegates through the system)
hv = vector_field(f)             # ❓ or symplectic_gradient(h) — not clearly exposed

Work item

• Add convenience getters and/or reconcile names:
• vector_field(flow) / hamiltonian_vector_field(flow) — expose the (symplectic) vector field X_H of a Hamiltonian flow directly (we have hamiltonian_vector_field on systems; make sure it is reachable from the flow and named intuitively).
• symplectic_gradient(h) — decide whether to offer this as an alias/name for X_H of a Hamiltonian.
• Audit the getter surface for consistency (flow → system → data delegation) and document the canonical names, keeping backward-friendly aliases where cheap.

---

8. Variable costate for variable times — confirm remaining cases

Issue: #231. Status: mechanism ✅; v = tf ✅ tested; v = t0 and v = (t0, tf) to confirm.

The semantics are settled (this closes the main v3 §7 open question): the flow always integrates the naive ṗv = -∂H/∂v with pv(t0) = 0, uniformly, for any variable — time or not. When the variable is a free time, the boundary term is a mitigated transversality condition written in the shooting method, not by the flow:

p_{t0}(tf) = -H(t0, x0, p0, v),    p_{tf}(tf) = H(tf, xf, pf, v).

Remaining work

1. Confirm/extend tests for v = t0 (free initial time) and v = (t0, tf) (both free) — the test_variable_costate_free_time.jl suite covers v = tf and a 2-D variable; make the t0-only and combined (t0, tf) cases explicit with analytic references.
2. Document the eval-time ≠ variable-value distinction ([Bug] Variable times #183) and the shooting-side transversality residual in the guide.
3. Extend to the with-control case (§1/§2) once implicit control lands — variable-costate for free times must work alongside DynClosedLoop and implicit-control flows.

---

9. Cache trajectory projections at construction (quality / perf)

Status: open, small, low-risk cleanup. Suggested from controlled_trajectory.jl.

Observation: the trajectory accessors rebuild an immutable projection functor on every call instead of returning a stored one. For example control(sol::ControlledTrajectory) constructs a fresh ControlProjection(sol.traj, sol.law, sol.variable, sol.state_coerce) each time, although all four fields are known at construction. Same for state(sol::ControlledTrajectory) (ControlledStateProjection).

Idea: build the projections once at trajectory construction and store them as fields; the accessors then just return the stored functor. This trades a little more memory in the struct (a couple of small immutable functors) for zero reallocation on repeated state(sol) / control(sol) calls (which happen a lot — plotting, objective quadrature, shooting).

Applicability — two cases, one subtlety

• ControlledTrajectory (ControlProjection, ControlledStateProjection) — the projections wrap sol.traj (the inner trajectory) + law + variable + coerce, all available before the ControlledTrajectory exists. ✅ Directly precomputable and storable.
• HamiltonianVectorFieldTrajectory (StateProjection, CostateProjection) — these wrap sol itself (self-reference), so they cannot be built before sol exists. Options: (a) re-target them to wrap the inner integration result instead of sol, then precompute; (b) a mutable/Ref field set right after construction; (c) leave as-is — they wrap a single field and are the cheapest case. ⚠️ Decide per-case; (a) is cleanest if the inner result is enough.

Steps

1. Audit all trajectory accessors that return a freshly-built projection (ControlledTrajectory, HamiltonianVectorFieldTrajectory, VectorFieldTrajectory, and the CTModels.Solution builder in OptimalControlFlow if it rebuilds closures).
2. For projections depending only on pre-existing data (the ControlledTrajectory case), add fields and store them at construction; accessors return the stored functor.
3. For self-referential projections, pick (a)/(b)/(c) above; prefer re-targeting to the inner result.
4. Keep the accessor API unchanged (state, control, costate still return callables).

---

10. Internal ODE solver (optional / strategic)

Issue: #81. Status: open, exploratory. Not committed for v4 — recorded here because it interacts with the integrator/strategy architecture (§1, §4, §5).

Idea: a standalone in-house RKF 5(4) adaptive solver to reduce the OrdinaryDiffEq.jl dependency (install/precompile weight), keep things simple, and get easier access to tailored AD (which helps §4's variational integration). Strategic; revisit if precompile cost or AD-through-the-integrator becomes a real pain point.

---

Priority Summary

High priority:

1. Implicit control (IFT/DAE) (§1) — the last missing branch of the with-control OCP story; unblocks the classic indirect-method use case. Start with IFT.
2. Getter naming / vector_field(flow) (§7) — small, closes [Dev] Add getter for hamiltonian vector field #185, unblocks downstream.
3. Non-Hamiltonian Flow(ocp) without p0 (§6) — direct-shooting use case, small.

Medium priority:

4. Flow derivative via dual numbers (§4) — variational/IND; decide option-vs-strategy jointly with §1.
5. Constraints design note (§3) — reflection deliverable; must precede any constrained flow. Then constrained dynamics (§2).
6. Variable-costate free-time t0 / (t0,tf) tests (§8) — finish [Dev] Augmented adjoints for variable times #231.

Lower priority:

7. GPU support + tests (§5) — important but independent; benefits from a clean array-type audit.
8. Cache trajectory projections (§9) — small, low-risk quality/perf cleanup.
9. Internal ODE solver (§10) — strategic, optional.

---

Cross-cutting decision: option vs. CTBase.Strategies strategy

Three features (implicit control §1, flow derivative §4, GPU §5) each face the same architectural choice: expose the behaviour as a keyword option (routed through the action-option machinery, like hamiltonian_type) or as a first-class CTBase.Strategies strategy (like the :di backend and :sciml integrator strategies).

➡️ Make this decision once, up front, and apply it consistently. A shared convention avoids three divergent mechanisms. Rule of thumb: strategy when the feature needs its own options and/or an extension boundary (weak deps: Sundials for DAE, CUDA for GPU); plain option when it is a local RHS branch with no new dependency.

---

Dependency Graph

Cross-cutting: option-vs-strategy convention ─┬─→ §1 (implicit control)
                                              ├─→ §4 (flow derivative)
                                              └─→ §5 (GPU)

§3 (constraints design note) ──→ §2 (constrained dynamics) ──→ constraint multipliers (if any)

§1 (implicit control) ──→ §8 (variable costate must work with control flows)

§7 (getter naming), §6 (basic Flow(ocp)), §9 (projection caching) ── independent, small

§10 (internal ODE solver) ── strategic, touches integrator architecture (§1/§4/§5)

---

Carried Over / still relevant

• Flows on generic objects (ODEProblem, ODEFunction) — phase 4 of the umbrella issue Complete Flow Construction features #247; partially present via the SciMLBase extension, audit completeness.
• Input validation at flow call — check user functions define the expected methods.
• Documentation & test revamp — follow the CTBase/Handbook guides; add guides for implicit control (§1), constrained flows (§2/§3 once decided), and flow derivatives (§4).