JuliaCon 2026 · Optimal Control
Optimizing race car track times in Dyad.
Problem
Drive one lap, as fast as possible
- Nonlinear optimal control with free final time:
minimize lap time
T - Two controls: thrust and steering angle
- Subject to
- vehicle dynamics
- tire force saturation
- engine power limit
- periodic boundary: end the lap the way you started
- Reference implementation: the Dymos (NASA / OpenMDAO) racecar example
Model
A symbolic DAE, authored in Dyad
- Dyad components compile to ModelingToolkit systems — the analysis below is plain SciML all the way down
- Longitudinal + lateral dynamics, load transfer, nonlinear tire saturation
- Track geometry: curvature
κ(s)is an interpolation object stored as a parameter mtkcompile(sys; inputs, outputs)— controls stay symbolic inputs, constrained quantities stay addressable outputs
After compilation
| states | 9 |
| control inputs (thrust, steering) | 2 |
| bounded outputs (tires ×4, power) | 5 |
| collocation nodes | 51 |
Costs and constraints are symbolic too:
costs = [-EvalAt(te)(…)], point constraints via
EvalAt, path bounds from variable metadata.
Attempt 1
The direct formulation, time as the independent variable
- Free
tf⇒ normalize to τ ∈ [0, 1], dynamics become∂U/∂τ = tf · f(U, V) tfmultiplies every defect constraint, so one free variable scales the whole NLP- Warm-starting from the exact converged racing line still fails: the ill-conditioning is intrinsic to the formulation
97
iterations, JuMP + RK4
(converges, barely conditioned)
(converges, barely conditioned)
Fail
InfiniteOpt finite differences
(oscillates, dual-inf → 3·10⁷)
(oscillates, dual-inf → 3·10⁷)
Fail
InfiniteOpt orthogonal collocation
(never reaches feasibility)
(never reaches feasibility)
Reparameterize
Drive on arc length, not on the clock
- One call:
change_independent_variable(sys, s), the s-parametrization used in everything that follows - Grid becomes fixed: s ∈ [0, track length]
- Time turns into a passive lap-timer state,
the objective is just
t(L) - No free multiplier on the dynamics ⇒ well-scaled NLP
- Same trick Dymos uses, but here it is a symbolic transformation of the same model, not a rewrite
Interface
From a System to a solver, four ways
System
mtkcompile, inputs + outputs
→
DynamicOptProblem
costs, constraints, bounds, guesses
→
JuMP
InfiniteOpt
CasADi
Pyomo
- JuMPTableaus → generated function → JuMP variables The transcription's Runge–Kutta tableau drives a RuntimeGeneratedFunction for the dynamics; JuMP decision variables are propagated straight through it to build the defect constraints.
- InfOptSubstitute Symbolics variables for InfiniteOpt ones Symbolic states, controls and registered operators are substituted with their InfiniteOpt counterparts — the transcription (orthogonal collocation, finite differences) is then InfiniteOpt's job.
- BVPMTK codegen into the SciML solvers
Not an NLP transcription at all:
BVProblemuses MTK's own codegen and hands the collocation to BoundaryValueDiffEq's MIRK/FIRK solvers.
Results · Accuracy
Direct collocation, three ways into Ipopt
| Transcription | Stack | Lap time | Iterations |
|---|---|---|---|
| Gauss–Lobatto (50, 3), L-BFGS | Dymos (Python) | 22.27 s | 79 |
| RK4 defects, analytic Hessian | JuMPDynamicOptProblem |
22.83 s | 39 |
| Orthogonal collocation (4) | InfiniteOptDynamicOptProblem |
22.00 s | 37 |
- All find the same racing line (V ≈ 71–78 m/s, hugging the inside at n ≈ −2.8 m)
- The spread is discretization error at 51 nodes, not different optima
- Symbolic model ⇒ exact sparse Jacobians and Hessians for free — half of Dymos' iteration count
- OC(4) caveat: δ spikes to 0.06 rad at corner exits — path bounds hold only at the 51 supports
One symbolic model, one
solve per backend — the transcription is a
keyword, not a rewrite.
Results · Runtime
Runtime performance benchmarks
1.43 s
Dyad via MTK & InfiniteOpt wall time
(37 iterations)
(37 iterations)
6.3 s
Dyad via MTK & JuMP wall time
(39 iterations)
(39 iterations)
242.6 s
Dymos wall time
(79 iterations)
(79 iterations)
- All variants solve in Ipopt — cards report Ipopt solver wall time;
same machine, sequential runs, matched settings
(
mu_init = 1e-3, monotone barrier, MUMPS) - 95% of Dymos' time is KKT solves, not Python model evaluation: pyoptsparse cannot pass a Hessian callback ⇒ forced onto L-BFGS, ~2m+2 KKT solves per iteration — the symbolic path hands Ipopt an exact sparse Hessian: half the iterations, each ~20× cheaper
- Why OC(4) wins: each defect row is a single RHS evaluation at a collocation point, where RK4 chains four stage evaluations into every constraint ⇒ flatter expressions for the solver's AD
The speedup is structural: it's what an analytic sparse Hessian
buys when the modeling layer can generate one.
Results · Validation
Same racing line, two independent stacks
Dymos GL(50,3)
MTK · JuMP RK4
MTK · InfOpt OC(4)
Identical apexes (n saturating at ±4 m), braking points and speed envelope (39–80 m/s). The only divergence is on the straights, where lateral position is nearly cost-free.
Results · Racing line
One lap, drawn from the solution data
Dymos · 22.27 s
MTK · InfOpt 22.00 s
Exact segment geometry · lateral offset n to scale (|n| ≤ 4 m) · draw-in paced by the optimal speed profile
Journey · 01
What it actually took, pt. 1: scaling is not optional
- The power constraint is 5 orders of magnitude off.
One O(10⁵) row next to O(1) tire rows ⇒ Ipopt disappears into
restoration phase and never comes back.
- Fix: per-variable
scales— dynamics as scaled residuals(x′ − f)/scale ~ 0, bounds on outputs lifted into auxiliary variables with scaled equalities.
- Fix: per-variable
- Where you start matters twice. Starting a lifted variable exactly on its bound ⇒ “locally infeasible”. Interior start ⇒ optimal in 37 iterations.
- One barrier knob. Default μ₀ = 0.1 parks the solve in a
corner-limited local optimum (46.7 s).
mu_init = 1e-3finds the racing line. Same problem, same start.
These are exactly the things a modeling language can automate:
nominal-value metadata → automatic residual & constraint scaling.
Journey · 02
What it actually took, pt. 2: the track is a parameter you call
- Curvature
κ(s)is aDataInterpolationsobject stored as a callable parameter — evaluated at a decision variable - Symbolic → JuMP tracing calls it with solver types it can't accept:
MethodError: *(::Symbolic, ::GeneralVariableRef) - Now on master: every parameter that is called in the equations is
auto-registered as a JuMP nonlinear operator
- arity read off the call site; works for bare callables and
FunctionWrappers - symbolic derivative rules picked up from Symbolics' registry ⇒ exact gradients for the operator
- a bare
sinsideκ(s)lowers onto the backend's own time/arc variable
- arity read off the call site; works for bare callables and
Interpolated lookup data in optimal control, with no manual operator
plumbing. PR #4909 · merged
Results · Summary
The journey, in one table
| Formulation | Backend | Lap time | Iters | Ipopt wall | Verdict |
|---|---|---|---|---|---|
| time domain, free tf | JuMP RK4 | 22.41 s | 97 | — | converges, ill-conditioned |
| time domain, free tf | InfiniteOpt (FD / OC4) | — | — | — | fails, intrinsically |
| s-parametrization | JuMP RK4 | 22.83 s | 39 | 6.3 s | robust |
| s-parametrization | InfiniteOpt OC(4) | 22.00 s | 37 | 1.43 s | best honest lap, fastest solve |
| s-parametrization, lifted DAE-BVP | MIRK4 + Ipopt | 21.84 s | 119 | 302 s | leakage-flattered; honest 22.02 s |
| reference | Dymos GL(50,3) | 22.27 s | 79 | 242.6 s | brackets the same optimum |
- Formulation dominated everything: the reparameterization mattered more than any solver option
- The BVP's 21.84 s is boundary-enforcement leakage, not a better line — with control bounds + warm start it lands at 22.02 s, right among the honest transcriptions
- Every row is the same symbolic model — the journey is a sequence of keyword arguments, not rewrites
Upstream
All of this has landed upstream
merged ModelingToolkit.jl
- #4399 parameter bindings & observed substitution
- #4406
inputs(sys)for optimal control - #4407 variable bounds → backend bounds
- #4425 nominal-value metadata
- #4702 input-aware codegen on compiled systems
- #4839 functional initial guesses (trajectories)
- #4840 observed-variable bounds lifting
- #4841 residual scaling of dynamics
- #4909 callable-parameter registration
- #4930 InfiniteOpt verbose-flag fix
- #4931 control-aware BVProblem dynamics
open in review now
- BVDE #521 mass-matrix MIRK/FIRK
- BVDE #589 OC constraint Jacobian
- PreallocationTools #198 typed tracer workspaces
All split from the #4394 umbrella — reviewable pieces, each validated against this benchmark.
ModelingToolkitBase 1.65.0 ships all of them — the collocation
workflow needs no patch branches.
Roadmap
Where this is going
- Dyad analyses: this workflow packaged as a DyadModelOptimizer analysis — model in, racing line out
- Control-aware BVProblems & the OC constraint Jacobian — the remaining pieces of the BVP path, next in the upstream queue
- Solver breadth: Optimization.jl decouples transcription from solver — Ipopt today, MadNLP next
- Automatic scaling from nominal metadata — turn pt. 1 of the war stories into a non-event
- Beyond racing: same machinery = trajectory optimization, parameter estimation with path constraints, free-time maneuvers
The gap between “a model you simulate” and “a model you optimize over”
is closing — symbolically.
Wrap-up
One model, three transcriptions, a racing line.
- Written once in MTK, reparameterized on arc length with one symbolic call, then solved through JuMP, InfiniteOpt or BoundaryValueDiffEq
- Matches Dymos on the optimum (22.27 s vs 22.00 s) and beats it on wall-clock, with the whole MTK side shipped in ModelingToolkitBase 1.65.0
- A future version of DyadModelOptimizer will provide an analysis for dynamic optimization, so this workflow becomes a single declarative analysis: model in, racing line out
DYNAMIC OPTIMIZATION TUTORIAL
docs.sciml.ai/ModelingToolkit/
stable/tutorials/dynamic_optimization
stable/tutorials/dynamic_optimization
DYAD DOCUMENTATION
help.juliahub.com/dyad/stable
Backup
Backup: the BVP's 21.84 s was a leak
MIRK bounds live on mesh nodes — so at s ≈ 0 the lifted BVP fires thrust ≈ 40 (normal ±2), hits 115 m/s and swings power to −4000 kW.
Dymos 22.27 s
BVP unbounded 21.84 s
BVP honest 22.02 s
- 01Finer mesh halves the spike, lands in a worse basin — 31.34 s
- 02Control bounds close the exploit, break the cold start — 39.5 s, max_iter
- 03Bounds + warm start: 22.024 s, Optimal, 24 iterations — leakage ≈ 0.19 s