[Problems] Proof of concept: a ModelingToolkit modeling, solved with Corleone, with our own bryson_denham

Agents: read the guide for agents in #310 first (order, difficulty, how to claim an issue, working method). The specification is in this description and in the comments below.

Updated on 2026-11-06 after #360 (convergence is reported, not tested) and #375 (reference objectives). Where an older comment contradicts this description, the description wins; what changed is listed in the update comment.

Goal

A proof of concept of a new modeling, ModelingToolkit: write our own bryson_denham as a ModelingToolkit system in this package, as an extension, and show how to solve it with Corleone (SciML, direct shooting with DiffEq), with retrieval, tests and documentation. This shows that every part of the infrastructure (registry, statuses, tests, crons, results pipeline, documentation) handles a modeling whose method differs. Then decide whether to write more problems this way.

Modeling and solving are separate, as for OptimalControl (problem written with @def, method chosen in solve): the modeling is the language the problem is written in, ModelingToolkit; Corleone is one way to solve a ModelingToolkit system, as Ipopt is for JuMP. Other methods could solve the same system (ModelingToolkit's own dynamic optimization interfaces, e.g. JuMPDynamicOptProblem, CasADiDynamicOptProblem, to check).

Part of #310, step 7 (last). Depends on #303, #304, #306, #314, #307.

What it is not

Not an import of SciML's OptimalControlBenchmarks catalogue. The problem is ours, rewritten: our dynamics, parameters, cost, initial guess; only the method changes (shooting with piecewise-constant controls instead of collocation). Their catalogue is a source of inspiration (MTK and Corleone usage, prepare_benchmark, solve_with_corleone, plotting) and of problems to port later.

Problem

bryson_denham as in ext/OptimalControlModels/bryson_denham.jl: ẋ₁ = x₂, ẋ₂ = u, x(0) = (0, 1), x(1) = (0, -1), x₁(t) ≤ 1/9, cost ∫ ½ u², parameters from ext/MetaData/bryson_denham.jl. Analytic optimum (ℓ = 1/9 ≤ 1/6): J = 4/(9ℓ) = 4, recorded as the analytic objective with its source (#379). The reference_objective of #375 is the objective of the trapezoidal discretization at the default grid (4.000143…): it does not apply to a shooting solve, whose discretization differs.

Inherent to shooting with Corleone, to state on the problem page: the state constraint is imposed at the points of a constraint grid only, and controls are piecewise constant on a control grid.

The ModelingToolkit system must be continuous. In SciML's catalogue, the constructor receives the grids and samples the state constraint itself (x(tᵢ) ≲ 1/9 for tᵢ in constraint_grid). Here:

• (a) write a continuous path constraint x₁(t) ≲ 1/9 in the system, if ModelingToolkit can represent it;
• (b) sample it on Corleone's constraint grid when preparing the solve.

First task: check (a). If ModelingToolkit cannot represent a continuous path constraint, fall back to passing the constraint grid at retrieval, as grid_size for JuMP, and say so here.

Scope

1. Modeling: entry ModelingToolkitBackend in the modeling table (#303); extension triggered by ModelingToolkit only; file ext/ModelingToolkitModels/bryson_denham.jl; status :draft first.
2. Retrieval: bryson_denham(ModelingToolkitBackend(); parameters=…) returns our generic Problem (#283); model(pb) is the continuous ModelingToolkit system (with its time span).
3. Solving with Corleone (not part of the extension; in tests, documentation and user code): build the control, shooting and constraint grids from a grid size (single shooting by default, multiple shooting as an option), prepare the Corleone model (as prepare_benchmark in SciML's catalogue), solve with Ipopt.
4. Initial guess: our shared :init (#306) is used: evaluated on the control grid for the controls, and at the shooting nodes for the states in multiple shooting, written into the parameter vector of the prepared Corleone model (model.ps, through pack / reconstruct as in Corleone's examples). Check first that this is possible.
5. Tests (test/ModelingToolkit environment, which also loads Corleone; label test ModelingToolkit backend; Saturday cron of #305):
   Three kinds of tests, as for the other modelings (#360):
   • does not break (every run): build the system; solve with Corleone on a small grid with a few iterations (Ipopt); no exception, finite objective;
   • model equivalence (every run): the cross-evaluation: project the optimal control of OptimalControl onto the piecewise-constant control grid, evaluate the Corleone model with it, compare cost and constraint violation with OptimalControl's result. This checks that the mathematical problem is the same although the methods differ; the comparison needs a tolerance derived from the grid (piecewise-constant controls, sampled constraint). Check first what Corleone lets us evaluate with imposed parameters;
   • convergence report (cron, never fails, no tolerance): solve on the default grid; report the solver status and the gaps to the analytic objective and to the OptimalControl objective.
6. Documentation: the ModelingToolkit documentation adapter of #307 / #314: the tab of the problem page shows how to get the system and solve it with Corleone; result in result.json; plot. Nothing else changes, which checks that the documentation handles every modeling generically.

Points to check first

• Environments: with an extension, users who do not load ModelingToolkit install nothing; but the test/ModelingToolkit and documentation environments must resolve the OptimalControl v2 stack and the SciML stack together (shared dependencies such as ForwardDiff or Symbolics). If the documentation environment cannot, compute these results in their own job (#314) and only read them in the documentation.
• Compat bounds for ModelingToolkit (extension) and Corleone (test and documentation environments), which move fast.

Noted while reading OptimalControlBenchmarks (for later)

Not needed here, since we write our own problems; to report to the authors if useful:

• robbins: the state constraint is x₃(tᵢ) ≥ 0 in a variable named grid_cons_x₁; the Robbins problem constrains x₁.
• double_oscillator: m₁ = 200 (ours 100), and the cost integrated over (0, 1) while the horizon is [0, 2π].
• bryson_denham: cost ∫ w² (no ½), state constraint sampled.

Done when

• Continuous path constraint in the ModelingToolkit system checked (a), or fallback documented.
• ModelingToolkitBackend in the modeling table; bryson_denham available in it, with our parameters; solved with Corleone using our initial guess.
• Tests pass in test/ModelingToolkit (does not break, cross-evaluation); the convergence report shows the gaps; status switched to :valid once the model is equivalent.
• Problem page shows the ModelingToolkit tab (get, then solve with Corleone) and results without specific code in the generators.
• Decision in this issue: write more problems with ModelingToolkit or not.