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).

[ocots]

Ensemble flows — multiple shooting, initial-guess sweeps, and ensemble optimal control

Expanding on §5 (GPU support, issue #249) — the local GPU design work (see PR #326 and the probe/gpu/ diagnostic tooling) validates two GPU execution models and explicitly defers the production ensemble/batch API on Flow to its own design note. Recording three concrete use cases here so they don't get lost before that note gets written. They are architecturally different from each other and will likely need different call shapes on Flow.

1. Multiple shooting (time-domain splitting)

Split [t0, tf] into subintervals, integrate each leg in parallel, then a synchronization/communication phase solves the matching (junction) conditions between legs — the classic Newton step on the shooting residual. The parallel-integration phase maps onto ensemble GPU (model 2 of the design report, DiffEqGPU); the synchronization phase is inherently coupled/sequential and stays a small CPU-side (or GPU-resident) linear solve between rounds of parallel integration.

2. Multiple shooting over many initial guesses (globalization)

No time splitting — instead, many independent initial guesses (p0, or the full shooting unknown) integrated in parallel to robustify Newton convergence, keeping whichever guess(es) converge. This is the most direct fit for DiffEqGPU.EnsembleProblem as already probed and validated on H200 (design report §6bis.1, run 11): genuinely independent trajectories, no coupling at all between them.

3. Ensemble optimal control (coupled via control and/or criterion)

A different problem class entirely — not independent trajectory solves, but one optimal control problem whose state is the concatenation of N small, structurally-identical subsystems, coupled through a shared control u(t) and/or a pooled criterion (average or worst-case over the ensemble). This is the classic "ensemble/robust control" setting (e.g. NMR/MRI pulse design robust to a spread of physical parameters).

Concrete example: MRI saturation for a single spin (saturation.html) generalizes to multi-epsilon-saturation — N spins/isochromats with different resonance offsets εᵢ, all driven by the same control, optimized jointly (e.g. average or worst-case saturation error over the ensemble).

Note the different GPU fit here compared to #1/#2: because the RHS is block-diagonal (each spin's dynamics depend only on its own state, the shared control, and the pooled criterion — not on each other), this is naturally a matrix-batch RHS call (u::AbstractMatrix, one column per spin), which HamIpAugRHS already supports today. And because N can be large (hundreds to thousands of isochromats), the concatenated state can genuinely be large enough for array-level GPU (model 1) to pay off — the rare case in optimal control where model 1's "large state" threshold from the design report is actually met, unlike a single small OCP.

---

Summary: #1 and #2 are ensembles of independent trajectories (DiffEqGPU EnsembleProblem, model 2); #3 is a single coupled system best served by the existing matrix-batch RHS mode and array-level GPU (model 1). The future batch/ensemble Flow API design note (§8) should start from these three cases rather than a blank page.

---

Tracking: Captured as issue #341 — Ensemble flows — batch API for multiple shooting, globalization, and ensemble optimal control. Related: #249 (GPU support), #111 (GPU for ODEs).

[ocots]

Roadmap v4 — remaining work (renumbered)

Status update: §5 (GPU support & testing — array-level GPU model validated end-to-end on H200 hardware, PR #326; full device support, tests, and docs shipped through v0.15.0-beta) and the last open case of §8 (variable-costate for the with-control case, issue #231) are now done, on top of the §2/§3/§6/§7/§9 items already marked done above. Re-listing what's left, renumbered:

1. Implicit control (IFT/DAE) — ∂H̃/∂u = 0 (was §1). Issue Implicit control via IFT/DAE — Flow(ocp) with ∂H̃/∂u = 0 #46.
2. Flow derivative via dual numbers (variational equations / IND) — (was §4). Issue Flow derivative via dual numbers — variational equations (IND) #93.
3. Ensemble flows (multiple shooting, initial-guess sweeps, ensemble optimal control) — (was item 4 above). Issue Ensemble flows — batch API for multiple shooting, globalization, and ensemble optimal control #341.
4. Internal ODE solver (optional / strategic) — (was §10). Issue [Dev] Faster ode solver #81.

The "Plus" items from the previous update are now closed: #340 (missing integration tests) and #343 (compatibility table page) are both done. #315 remains the general roadmap tracking issue.

---

Release: v0.16.0-beta is out (bump + CHANGELOG merged in 97b3142) — fixes batched (Matrix)/SVector states on Flow(ocp)/Flow(ocp, law) (#362), and extends the "1-D = scalar" convention to Flow(::VectorField) and Flow(::ODEFunction) (#363 — a breaking output-shape change for those two constructors). Full details in the CHANGELOG.

Bubbling this up into a new OptimalControl.jl beta is the logical next step, but it's currently blocked: CTDirect.jl#610 point 7 ("1-D = scalar") hasn't landed in CTDirect yet, so the convention is inconsistent between the direct (CTDirect) and indirect (CTFlows) legs of the ecosystem. Once that item is done, a new OptimalControl.jl beta becomes reasonable.

Original item numbering (§1–§10) in the opening post stays the reference for full context on each item; this comment is just a shorter pointer to what remains.