[ChrisRackauckas-Claude]

Note

Rebased 2026-09-10 onto bf1efdb992, and the outstanding review is addressed — see the comment below for the per-thread mapping. Versions bumped to ModelingToolkit 11.43.0 / ModelingToolkitBase 1.70.0.

Since the last push the traversals were rewritten onto IRStructure (the homotopy operator's template), the limitnew/limitold sentinels became GlobalScoped — which deletes the name-suffix canonicalization machinery outright — and the corrector now compiles to a single generated function instead of one per limiter.

The first CI run also exposed three plumbing breakages, all fixed here: a missing NonlinearSolveBase test dependency on the root package, postcondition floors that were too low for the downgrade job (so limiting was silently inert there — 179 Newton steps instead of 11), and an ExplicitImports finding on Base.setindex.

Important

This PR should be ignored until reviewed by @ChrisRackauckas.

Summary

Adds the limited(actual, limiter) operator: symbolic, component-level declaration of SPICE-style iterate limiting, lowered automatically to the Predictor/Corrector Newton-Raphson (PCNR) formulation of Aadithya, Keiter & Mei. A device author annotates a sensitive quantity once —

i ~ Is * (exp(limited(v, pnjlim(limitnew, limitold, Vt, vcrit)) / Vt) - 1)

— and every nonlinear solve built from any model containing that component gets predictor/corrector limiting fully automatically: mtkcompile performs the PCNR augmentation and NonlinearProblem construction compiles the limiters into the postcondition corrector hook the problem carries.

The operator follows the architecture of the Modelica homotopy(actual, simplified) operator (#4601): the actual expression first, the helper second; limitnew/limitold are reserved placeholder symbols (fixed sentinels, precompile-safe like HOMOTOPY_LAMBDA) for the proposed and previously-accepted values inside the limiter expression.

What the lowering does

For time-independent systems (in both the ModelingToolkitBase light pipeline and the ModelingToolkit structural pipeline, inserted post-expand_connections):

1. each unique limited(actual, limiter) node gets an auxiliary irreducible unknown limited_k; the node is replaced by it and the consistency equation limited_k ~ actual is appended (EquationSourceInformation is padded accordingly for TearingState);
2. because limited_k is irreducible, alias elimination keeps the limited quantity as the surviving representative of its alias class — the PCNR augmentation is reduced symbolically, at no runtime cost;
3. limited_k receives actual as both default and guess, so no user-provided initial value is needed;
4. the limiter registry is stored in LimitedCtx system metadata; problem construction compiles the whole registry with one generate_custom_function call (parameters, including bound parameters, resolve through standard codegen) and merges it into the problem's keywords as postcondition, so it rides to solve/init like any other solver option and an explicitly passed corrector wins over the model-declared one.

For time-dependent systems the operator is stripped to actual during mtkcompile, so component libraries carrying limiters compile unchanged for transient simulation — but the stripped limiters are remembered in StageLimitedCtx metadata, and an ODEProblem built with nlstep = true re-attaches them to the implicit stage system, conjugated by the affine stage substitution q = a*z + b into φ(znew, zold) = (L(a*znew + b, a*zold + b) - b) / a. Un-lowered limited nodes are rejected with clear errors at NonlinearFunction/ODEFunction construction; nested limited and limiters referencing other unknowns are rejected likewise.

The same conjugation delivers the bounds metadata of an ODE state to the stage solve as a per-stage clamp. That fixes a pre-existing bug in passing: generate_ODENLStepData used to let a state's static lb/ub land on the stage problem, whose unknowns are the increments z — the wrong variable, and enough to put the problem into NonlinearSolve's transformed coordinates that the raw stage values the stepper writes back do not live in.

Changes

• ModelingToolkitBase (1.69.1 → 1.70.0): new src/systems/limited_operator.jl (operator, GlobalScoped sentinels, IR-backed detection/strip/lowering passes, attach_stage_limiters, generate_limited_postcondition); lowering invoked from __mtkcompile; postcondition merged into the problem keywords with guards; DeriveBounds sentinel so lb = nothing, ub = nothing is distinguishable from "not supplied"; nlstep_limit_bounds threaded through ODEFunction; exports limited, limitnew, limitold; tests in test/limited_operator.jl (registered in the InterfaceII group).
• ModelingToolkit (11.42.1 → 11.43.0): lowering invoked from the structural __mtkcompile (with source-info padding); inner_nlsystem attaches the stage limiters; nlstep_scc = true rejects declared limiters rather than silently dropping them.
• Docs: new docs/src/basics/Limiting.md (registered in pages.jl) with a component-based diode-circuit walkthrough, the nlstep story, lowering semantics, and contracts.
• Test plumbing: NonlinearSolveBase added to the root package's test target; postcondition-era floors for NonlinearSolve / NonlinearSolveBase / OrdinaryDiffEqNonlinearSolve in both Project.tomls and docs/Project.toml; Base.setindex added to the all_qualified_accesses_are_public ignore list.

Verification (all local, Julia 1.12.7, registered dependencies — nothing unmerged is needed)

Run	Result
lib/ModelingToolkitBase/test/limited_operator.jl with ModelingToolkit loaded (CI's group_interfaceii.jl path)	96/96, 19 testsets, exit 0
Root GROUP=InterfaceII julia --project -e 'using Pkg; Pkg.test()'	1234 passed, 7 broken, 0 failed, 54m13s, Testing ModelingToolkit tests passed
lib/ModelingToolkitBase GROUP=QA	20 passed, 1 failed, 0 errored — byte-identical to master at the rebase base bf1efdb992 (JET report_package, 234 reports, every one a Moshi @match temporary; tracked in #5063). The extra error the previous run had here (Base.setindex is not public) is gone.
homotopy_lowering.jl (shares the touched lowering pipeline)	43/43
Downgrade floors — env pinned to NonlinearSolve 4.29.1, NonlinearSolveBase 2.45.0, OrdinaryDiffEqNonlinearSolve 2.8.0, MTKBase only	limited_operator.jl 66/66 (the nlstep testsets gate off without ModelingToolkit, as in downgrade-mtkbase). Before the floor bump this run gave 179 Newton steps instead of 11 and calls[] == 0.
The three @example limiting blocks of docs/src/basics/Limiting.md, extracted and executed	all three run clean
runic --check and typos over every changed file	clean

Not verified locally: the full docs build (only the new page's example blocks were executed), the Julia lts/pre jobs, GPU and downstream jobs.

Things a reviewer should push back on

• The stage-limiting path requires each limited quantity to be affine in one stage unknown. A declared limited quantity that is not is an ArgumentError; a bounds-derived clamp that is not is silently dropped. That asymmetry is deliberate — the user asked for the first and not the second — but it is a judgment call.
• bounds clamping on the stage solve is on by default (nlstep_limit_bounds = true). The argument is that clamping an intermediate Newton iterate cannot move the root, only the path, and it keeps log/sqrt/exp right-hand sides in domain. If you disagree, the default should flip.
• NonlinearSolve and OrdinaryDiffEqNonlinearSolve floors were raised (test/docs targets only) rather than the tests being made tolerant of an older stack.
• The fully implicit DAE stage path (inner_dae_nlsystem) deliberately gets neither limiting nor a derived box, preserving its pre-limited behaviour. Extending it is untested work left out on purpose.

Design notes

• Sentinels are parameters (not variables) so the time-dependent System(eqs, t) constructor's function-of-t validation accepts equations containing them, and GlobalScoped so a limiter written inside a component still refers to the two toplevel symbols after flattening; both lowering paths drop them from the discovered unknowns/parameters by identity.
• Every traversal goes through IRStructure — SU.query for detection, search_variables! with an is_atomic predicate for node collection, IRSubstituter for every rewrite — following the template the homotopy operator was refactored onto. This matters because has_any_limited runs on every ODEFunction/NonlinearFunction construction.
• The whole limiter registry compiles to a single generated function: correctors registered against the same unknown are composed symbolically before codegen, so the corrector is concretely typed and its compile time does not grow with the number of limited quantities in a device library.
• The lowering keeps constructor-built caches coherent (var_to_name, and crucially the irreducibles field, which structural simplification consults rather than per-variable metadata — found the hard way).
• Derivative rules register ∂₁ = 1, ∂₂ = 0: Jacobians treat the operator as actual, matching the hook semantics (the limiter is a corrector, not part of the residual).

Closes the ModelingToolkit third of the SciML/NonlinearSolve.jl#351 + PCNR feature arc.

🤖 Generated with Claude Code

https://claude.ai/code/session_01T2Y8c85c7AiQHEjK9KLzdA

[ChrisRackauckas-Claude]

Verification detail (local, Julia 1.12.4). The full GROUP=InterfaceII harness run is currently blocked by registered-ecosystem version skew, not by this PR: dev'd SciMLBase master requires FindFirstFunctions 3 while the MTKBase test env pins DataInterpolations 8 (needs FFF <2), and dev'd NonlinearSolve requires LinearSolve 5 while registered OrdinaryDiffEq pins <5. Any MTK branch tested against SciMLBase master hits the same wall until the compat sweep accompanying those releases.

What was run instead, in a coherent dev environment (SciMLBase#1449 + NonlinearSolve#1084 + this branch):

Testset	Result
test/limited_operator.jl (new, MTKBase light pipeline)	28/28
test/homotopy_lowering.jl (shares the touched pipeline)	43/43
test/homotopy_problem.jl (exercises the edited problems/nonlinearproblem.jl)	29/29
test/nonlinearsystem.jl	23/24 — the 1 error is an using OrdinaryDiffEq precompile failure from the version skew above (old OrdinaryDiffEqDifferentiation vs SciMLBase master), unrelated to this diff
Full ModelingToolkit structural pipeline (root __mtkcompile) end-to-end	flat diode: 179 → 11 steps; hierarchical diode+resistor components (namespaced + bound vcrit parameter): 11 steps, correct v/i via observed; transient System(eqs, t): stripped, compiles, integrates, RHS matches analytic value

Also verified the operator's derivative rules (∂₁=1, ∂₂=0, Jacobians see actual), nested-limited rejection, the un-lowered guard on NonlinearFunction/ODEFunction, and the limiter-references-unknowns error.

🤖 Generated with Claude Code

[ChrisRackauckas-Claude]

CI triage (first full run): Runic Format Check and Spell Check pass. All 83 failing jobs (tests/sublibrary-ci/Downgrade/Documentation/Benchmarks/Catalyst downstream) fail at dependency resolution with Unsatisfiable requirements for SciMLBase — the documented state until SciML/SciMLBase.jl#1449 is merged and 3.37.0 is tagged (this PR's compat requires it). Verified by reading the InterfaceII and Downgrade logs directly; no non-resolution failure modes present.

Release order for green CI here: SciMLBase#1449 → tag 3.37.0 → re-run (tests pass with limiting inert) → NonlinearSolve#1084 + release → hook-gated solver-behavior assertions activate.

🤖 Generated with Claude Code

[AayushSabharwal]

This should use https://github.com/JuliaSymbolics/SymbolicUtils.jl/blob/master/src/irstructure.jl#L659 to have any hope at scaling to Multibody models. Naive recursion on an expression is strictly off the table in such cases.

[AayushSabharwal]

In a similar fashion, this should use IRStructure, though we don't have an inbuilt utility for this. An alternative formulation is to propagate a cache::Base.IdDict{SymbolicT, SymbolicT} throughout the recursion, populate it with the result of maketerm and check it at the beginning to avoid duplicating work.

[AayushSabharwal]

Similar cache required.

[AayushSabharwal]

Isn't this just substitute? It should use IRSubstituter.

[AayushSabharwal]

This is just search_variables!. It should use the IRStructureSearchBuffer version.

[AayushSabharwal]

Can't the lowering use Unicode or a name with # in it?

[AayushSabharwal]

Using IRStructureSearchBuffer and get_variables! here will be faster.

[AayushSabharwal]

Can these closures be made callable structs so we don't accidentally forget to add something to the let block later?

[AayushSabharwal]

Do we want to FunctionWrapper these somehow to cut down on compile times?

[ChrisRackauckas-Claude]

Thanks @AayushSabharwal — all addressed in 204a5ce. Every traversal now goes through the IR machinery rather than naive recursion:

Comment	Change
_has_limited naive recursion (scaling to Multibody)	has_limited is now SU.query(IsLimitedNode(), expr) — DAG-memoized, plus a callable-struct predicate reused as the is_atomic below
_strip_limited needs IRStructure/cache	strip is now node collection + one SU.IRSubstituter pass over equations and observed
_replace_limited_nodes "isn't this just substitute?"	yes — replaced by IRSubstituter{false}(ir, rules); no hand-rolled maketerm rebuild left
_collect_limited_nodes! "just search_variables!"	search_variables! with IRStructureSearchBuffer and is_atomic = IsLimitedNode(), so the limited nodes themselves land in the buffer. Kept an OrderedSet so insertion order stays deterministic — the generated names have to be a pure function of the system for RGF-cache/precompile reasons
sentinel canonicalization recursion	same treatment: search_variables! with an IsSentinel atomic predicate → substitution rules → IRSubstituter
free-variable check → IRStructureSearchBuffer/get_variables!	done, with the buffer reused across limiters
"Can't the lowering use Unicode or a name with #?"	now Symbol("#limited_", k) — unwritable as a Julia identifier, so the name-collision guard is deleted outright
closures → callable structs	LimitedPostcondition{iip, L} with an explicit limiters field

On the last one — FunctionWrapper for the compiled limiters: I left them unwrapped for now, because I don't think there's a signature that survives the use sites. The limiter is called as fn((u_new, u_old), p) where the value type must stay generic (ForwardDiff Duals flow through the same entries when the outer problem is differentiated, e.g. parameter estimation over a limited circuit) and p is an MTKParameters whose type varies per system, so a concrete FunctionWrapper{Float64, Tuple{Tuple{Float64,Float64}, P}} would either break AD or need one wrapper type per system anyway. The tuple-of-limiters does specialize on limiter count, so compile time grows with the number of distinct limited quantities in a model — if that bites on a large device library, the cheap fix is a length threshold above which we switch to a Vector of a single wrapper type with a per-eltype cache. Happy to add that upfront rather than on measurement if you'd prefer.

Re-verified after the rewrite: limited_operator.jl 28/28, and the full ModelingToolkit structural pipeline end-to-end — flat diode 179 → 11 steps, hierarchical component circuit 11 steps with identical root (0.6698509496766559), transient strip path unchanged.

🤖 Generated with Claude Code

[ChrisRackauckas-Claude]

Updated for the SciMLBase/NonlinearSolve redesign (commit febcd0c82): precondition/postcondition are now solve keywords rather than NonlinearFunction fields — see SciML/SciMLBase.jl#1449 for the rationale.

For this PR the change is small and strictly an improvement: the corrector compiled from the limiter registry is merged into the problem's keywords by merge_limited_postcondition instead of being threaded through NonlinearFunction, so it rides to solve/init like any other option. Two consequences worth noting:

• A corrector passed explicitly at solve time now takes precedence over the model-declared one, so a user can override a device library's limiting without rebuilding the model. merge_limited_postcondition deliberately leaves a user-supplied postcondition in place rather than clobbering it.
• NonlinearFunction(sys) no longer needs the expression = Val{true} special case for limited systems at function-construction time; the guard moved to the problem-construction path where the closure is actually created.

Re-verified after the change: limited_operator.jl 29/29 (up from 28, with a new test asserting the solve-time override actually replaces the model's corrector), and the full ModelingToolkit structural pipeline still gives 179 → 11 Newton steps on the flat diode circuit and 11 steps through the hierarchical component circuit.

🤖 Generated with Claude Code

[ChrisRackauckas-Claude]

limited now works inside transient solves, via NLStep

limited was inert in a transient simulation — stripped to actual because there was
nowhere to hang a corrector. There is now: src/systems/solver_nlprob.jl builds a genuine
time-independent System for the implicit stage equations and wraps its NonlinearProblem
in SciMLBase.ODENLStepData, and OrdinaryDiffEqNonlinearSolve solves that through
NonlinearSolve's own init cache path. So an ODEProblem built with nlstep = true now
carries the model's limiters into the stage Newton.

The obstacle that shaped the design

The obvious plan — re-attach the limited nodes to the stage system and let the existing
lower_limited PCNR augmentation run — does not work, for a reason outside this repo.
ODENLStepData requires the stage system's unknowns to be an index subset of the ODE
unknowns:

• MTK: subsetidxs = [findfirst(isequal(y), unknowns(sys)) for y in unknowns(nlsys)]
• OrdinaryDiffEqNonlinearSolve/src/newton.jl:87: nlstep_data.nlprob.u0 .= @view z[nlstep_data.u0perm]
• newton.jl:183-195: atmp_sub = @view(atmp[u0perm]), and the Newton convergence norm is
  taken over that view.

The augmentation introduces an auxiliary irreducible unknown #limited_k per limited
quantity. That is not an ODE state: findfirst returns nothing, the u0 fill and the
residual index map both break, and there is no entry in z to seed it from. Making it work
would mean changing the ODENLStepData contract in SciMLBase and
OrdinaryDiffEqNonlinearSolve.

What is done instead

The augmentation is not needed in the stage system. The stage substitution
v -> γ₂*v + inner_tmp[i] already makes each limited quantity an affine function
q = a*z + b of a stage unknown, with a/b built from parameters (γ₂, inner_tmp[i])
only. Correcting q is correcting z, so the corrector is the model's limiter L
conjugated by that affine map:

φ(znew, zold) = (L(a*znew + b, a*zold + b) - b) / a

φ is registered against z in the same LimitedCtx metadata the time-independent
path uses, so merge_limited_postcondition / generate_limited_postcondition /
LimitedPostcondition are reused verbatim — no new codegen. For the RC-diode below the
registry that lands on the stage system is literally

v(t) => (-inner_tmpₘₜₖ[1] + pnjlim(inner_tmpₘₜₖ[1] + __limitnew_ₘₜₖ*γ₂ₘₜₖ,
                                   inner_tmpₘₜₖ[1] + __limitold_ₘₜₖ*γ₂ₘₜₖ, Vt, 0.71)) / γ₂ₘₜₖ

φ inherits L's fixed-point property, and the corrector runs before the residual is
evaluated at each iterate, so a converged stage still solves the unmodified stage equations
— limiting changes how the step is found, never what it is.

Restrictions, both errors rather than silently dropped limiting:

• the limited quantity must resolve to an affine function of one stage unknown;
• nlstep_scc = true is rejected — the SCC decomposition splits the problem the corrector
  is attached to.

A limited quantity that no stage unknown feeds is skipped: the solver never iterates on it.

Surviving mtkcompile

mtkcompile still strips limited to actual for time-dependent systems, byte for byte as
before; strip_limited_system additionally records actual => limiter in a new
StageLimitedCtx system metadata key.

I picked this over keeping the nodes through structural simplification (or making the strip
lazy) because a limited(v, ...) call node is opaque to the structural pass — carrying
it would change tearing and alias elimination for every time-dependent model bearing the
annotation, including ones that never build NLStep data. Metadata is inert, so nothing about
transient compilation changes; the test asserts the generated ODEFunction right-hand side
is still exactly actual.

The recorded actual never has to survive literally in the simplified equations.
attach_stage_limiters resolves it with substitute_observed(sys, ·) → the stage subrules
→ substitute_observed(stagesys, ·), walking it through both compilations' eliminated
variables into stage coordinates.

What limitold means in a stage solve

NonlinearSolve passes cache.u_cache — the previous accepted iterate — as u_prev to
apply_postcondition!!. So inside a stage solve, limitold is the previous Newton iterate
of the current implicit stage, evaluated as the physical stage value a*z_prev + b. It is
not the previous time step. That is the SPICE reading of limiting: it damps the iteration,
not the trajectory; the first iterate of a stage compares against the step's predictor.
Documented in the limited docstring, the limitold docstring, and a new section of
docs/src/basics/Limiting.md.

Also in this push

• Dropped the NLS_APPLIES_POSTCONDITION test gate now that NonlinearSolveBase v2.41.0 and
  NonlinearSolve v4.25.0 are registered with postcondition; those assertions run
  unconditionally. Compat bumped to match (NonlinearSolveBase = "2.41",
  NonlinearSolve = "4.25"), and NonlinearSolveBase added to the root test target.
• Removed five stray [deps] entries (NonlinearSolve, NonlinearSolveBase,
  NonlinearSolveFirstOrder, NonlinearSolveQuasiNewton, NonlinearSolveSpectralMethods)
  that the first commit accidentally added to ModelingToolkitBase/Project.toml — its source
  never imports them and they were already [extras]. (Aqua would have flagged them.)
• limited_operator.jl is now also run from the root suite (test/group_interfaceii.jl via
  @mtktestset), which is what makes the nlstep testsets execute — they are gated on
  @isdefined(ModelingToolkit) since generate_ODENLStepData lives in ModelingToolkit.jl.
  (This incidentally gives NLStep its first test in this repo; grep -rn nlstep test/ was
  previously empty.)
• NonlinearSolveAlg is reached as OrdinaryDiffEqNonlinearSolve.NonlinearSolveAlg in both
  the test and the docs: OrdinaryDiffEq v7 no longer re-exports it and the owning package
  never exported it, so there is no public path to the one stepper option that consumes
  nlstep_data.

Local test results

Environment: this branch's MTK + MTKBase (dev), SciMLBase 3.43.0 (dev), registry
NonlinearSolve 4.25.0 / NonlinearSolveBase 2.41.0, OrdinaryDiffEqSDIRK +
OrdinaryDiffEqNonlinearSolve. Julia 1.12.4. Run after Runic.

Test Summary:                                    | Pass  Total     Time
with ModelingToolkit loaded (root package suite) |   53     53  1m49.1s
  Limited operator (PCNR iterate limiting)       |   53     53  1m49.0s

Test Summary:                                       | Pass  Total   Time
without ModelingToolkit (ModelingToolkitBase suite) |  127    127  51.9s
  Limited operator (PCNR iterate limiting)          |   42     42   3.1s
  Homotopy lowering                                 |   43     43   9.3s
  Homotopy problem construction & sweep             |   29     29  33.4s
  Homotopy OMC parity                               |   13     13   6.1s

limited_operator.jl went 29 → 42 assertions standalone, 53 with ModelingToolkit loaded
(the two nlstep testsets). New coverage:

• the StageLimitedCtx registry survives mtkcompile and the ODEFunction right-hand side
  is unchanged;
• attach_stage_limiters produces the conjugated limiter, checked numerically: with
  a = 0.5, b = 0.1, the corrector maps a proposed z = 3.0 to
  (pnjlim(0.5*3.0 + 0.1, 0.5*0.0 + 0.1, …) - 0.1) / 0.5, and is the identity at a fixed
  point;
• rejection of a quantity spanning two stage unknowns, and of one nonlinear in its stage
  unknown;
• ODEProblem(..., nlstep = true) puts a postcondition on nlstep_data.nlprob.kwargs,
  u0perm is unchanged ([1], no auxiliary unknown), and — the assertion that actually
  proves the wiring — init(nlstep.nlprob, NewtonRaphson()) gives a cache whose
  NonlinearSolveBase.get_postcondition is the generated LimitedPostcondition;
• an ImplicitEuler(nlsolve = NonlinearSolveAlg()) transient solve succeeds and agrees with
  the unannotated model to rtol = 1e-6;
• nlstep_scc = true throws.

Manual checks beyond the suite, on the RC-diode (Vsrc = 5V, R = 1kΩ, C = 1nF, Is = 1e-14, Vt = 25mV, vcrit = 0.71), ImplicitEuler + NonlinearSolveAlg, adaptive = false:

	limited	unannotated
dt = 1e-6, tspan = 20dt	v = 0.0, ConvergenceFailure	v = -5.1e16, "Success"
adaptive to t = 1e-5	nf = 652, v = 0.6698509496766557	nf = 661, v = 0.6698509496766557

The first row is the limiter doing its job — the unannotated stage Newton walks off into
exp(v/Vt) and returns garbage, the limited one refuses to accept the step. The second is
the honest caveat: on a smoothly stepped RC ramp the limiter almost never binds, so it buys
~1% of f-evaluations, not the 179 → 11 iteration collapse the DC operating point shows. The
value here is robustness on hard turn-ons, not speed on easy ones.

[ChrisRackauckas]

The discussions from before miss that limiting is a function of operators, not variables or even equations. You don't want to limit a full equation, and the choice of limiter can be operator-specific. So doing this based on variable bounds is just a bad approach because you don't necessarily know the function form to do with that. Even if we do something more automatic, I think we want this manual form and the automatic form just lowers to this, so it still needs to be exposed.

[ChrisRackauckas-Claude]

bounds metadata was leaking onto the stage problem as a wrong static box

While extending the nlstep stage-limiting path I found a pre-existing bug worth fixing here, since it is the same machinery.

NonlinearProblem(sys, op) auto-derives lb/ub from the unknowns' bounds metadata whenever they are not passed explicitly. inner_nlsystem builds the stage system with new_dvs = unknowns(sys) — the same symbolic variables, carrying the same metadata — so an ODE state's box silently landed on the stage problem. Measured on D(x) ~ -k*x^2*(1-x) with @variables x(t) = 0.5 [bounds = (0.0, 1.0)] and nlstep = true, before the fix:

stage unknowns : x(t)      stage u0 : [0.5]
stage lb : [0.0]           stage ub : [1.0]
NonlinearSolveBase.needs_bounds_transform(nlprob, NewtonRaphson()) = true
transform maps u0 [0.5] -> [0.0]

That is wrong twice:

1. Wrong variable. The stage unknown is the increment z; the physical quantity is γ₂*z + inner_tmp[i]. lb ≤ u ≤ ub is not lb ≤ z ≤ ub. The correct stage box is (lb - b)/a ≤ z ≤ (ub - b)/a with a = γ₂, b = inner_tmp[i] — and the stepper rewrites a and b every stage, so the correct box is time-varying and cannot be a static vector.
2. Coordinate corruption. Bounds trigger NonlinearSolve's reparametrization, but OrdinaryDiffEq keeps writing raw stage values into the problem each stage (OrdinaryDiffEqNonlinearSolve/src/newton.jl: nlstep_data.nlprob.u0 .= @view z[u0perm], then reinit!). Raw values into a transformed problem is silent corruption.

(A linear ODE cannot show this: mtkcompile solves the stage system exactly and tears it to zero unknowns, so there is no stage unknown for a box to land on. The tests use a nonlinear RHS for this reason.)

The fix

Stop the wrong static box. generate_ODENLStepData now builds the stage NonlinearProblem with lb = nothing, ub = nothing. That needed a sentinel default (DeriveBounds()) on the NonlinearProblem/NonlinearLeastSquaresProblem constructors, because previously an explicit lb = nothing, ub = nothing was indistinguishable from "not supplied" and still triggered derivation. Derivation behaviour is unchanged for every other caller — the existing Bounds metadata is forwarded to NonlinearProblem/NonlinearLeastSquaresProblem testset passes untouched (17/17).

Re-deliver the box correctly, through the machinery that is already there. A clamp is just a limiter, L(new, old) = clamp(new, lo, hi), so bounds on a state becomes a limiter of that state and goes through the existing _conjugate_limiter path, landing on the stage unknown as

φ(znew, zold) = (clamp(a*znew + b, lo, hi) - b) / a

The key property that makes this correct: a and b are symbolic expressions in the stage parameters (γ₂, inner_tmp), and the conjugated corrector is compiled by generate_custom_function against the system's parameters — so they are read live at call time. The box the corrector enforces on the increment moves with the stepper, for free. The tests assert exactly this by driving set_γ_c/set_inner_tmp the way the stepper does and checking that the same proposed z gets projected differently.

Only the finite sides are emitted (max(lo, min(limitnew, hi))), so a one-sided bound costs one comparison and an unbounded state costs nothing.

Composition. LimitedCtx is a Vector{Pair}, and LimitedPostcondition applies entries in order on the running value, so several correctors on one stage unknown compose rather than clobber. Bounds clamps are appended after the model's own limiters, giving clamp ∘ limiter: the limiter damps the proposed move, the clamp then projects whatever it produced, so the corrected iterate always lands in range. The test pins the order — with γ₂ = 0.5, b = 0.1, a halving limiter and a [0, 1] box, z = 6 gives physical 1.0 as clamp ∘ limiter and would give 0.55 the other way round.

Automatic, not opt-in

Bounds-derived clamping is on by default. Clamping an intermediate Newton iterate cannot change the root a stage solve converges to, only the path it takes there, and it keeps a right-hand side containing log/sqrt/exp from being evaluated outside its domain — which would NaN the whole solve rather than merely fail a step. That is a safety property, and safety properties belong on by default.

nlstep_limit_bounds = false remains as an escape hatch, kept for one concrete reason: a clamp is a projection, not just damping, so if a model's stage solutions legitimately leave the declared box (bounds used as an advisory range rather than an invariant) the projection would keep the stage solve from converging. One keyword lets such a model opt out without deleting its metadata.

Two consequences follow from "automatic" that would otherwise be regressions, and are handled explicitly:

• A box that cannot be conjugated (a bounded state the stage system tears away as a nonlinear function of other unknowns) is dropped, not raised. A declared limited(...) quantity still errors in that situation — the user asked for limiting and cannot have it — but a model that integrates today must not stop integrating merely because a state carries a box.
• nlstep_scc = true skips bounds clamping rather than rejecting the model, for the same reason. Declared limited quantities keep erroring there (the SCC decomposition splits the problem the postcondition is attached to).

Tests

lib/ModelingToolkitBase/test/limited_operator.jl, run in both modes it is run in CI.

With ModelingToolkit loaded (as test/group_interfaceii.jl runs it) — 96/96, up from 53:

Test Summary:                                                              | Pass  Total     Time
Limited operator (PCNR iterate limiting)                                   |   96     96  2m19.0s
  runtime numeric fallback is the actual branch                            |    1      1     0.0s
  derivative treats the operator as actual                                 |    3      3     2.3s
  lowering: augmentation, irreducibility, guesses, registry                |    7      7     0.5s
  nested limited operators are rejected                                    |    1      1     0.1s
  compiled system solves with automatic PCNR limiting                      |    6      6    39.5s
  hierarchical components: limiting composes through namespacing           |    5      5     5.0s
  time-dependent systems strip the operator but record it                  |    6      6     5.6s
  stage limiters: the limiter is conjugated onto the stage unknown         |    6      6     3.2s
  bounds become a clamping limiter of the physical quantity                |    8      8     7.0s
  stage limiters: bounds are conjugated into stage coordinates             |    7      7     0.4s
  stage limiters: a user limiter and bounds compose as clamp ∘ limiter     |    5      5     0.4s
  stage limiters: a box that cannot be conjugated is dropped, not an error |    1      1     0.1s
  explicit `lb`/`ub` of `nothing` suppress the metadata-derived box        |    3      3     0.4s
  stage limiters: quantities that cannot be conjugated are rejected        |    4      4     0.7s
  nlstep transient solve limits its stage Newton iterates                  |   10     10    37.5s
  nlstep rejects limiters it cannot attach                                 |    2      2     3.8s
  nlstep does not put the state's box on the stage increments              |    6      6     4.3s
  nlstep delivers the box as a per-stage clamp, automatically              |   12     12    22.9s
  guards                                                                   |    3      3     0.3s

Standalone (MTKBase only, the lib/ModelingToolkitBase/test/runtests.jl path) — 66/66, up from 42.

Existing bounds-derivation coverage, unaffected:

Test Summary:                                                                     | Pass  Total   Time
Bounds metadata is forwarded to `NonlinearProblem`/`NonlinearLeastSquaresProblem` |   17     17  16.2s

And the probe from the top of this comment, after the fix:

stage lb            : nothing
stage ub            : nothing
>>> bounds are NOT attached to the stage problem.

Runic clean on all five changed files. Documented in docs/src/basics/Limiting.md (new section, with a runnable @example block I verified prints (nothing, true)).

[ChrisRackauckas-Claude]

Rebased onto current master (de05533509). The branch was 221 commits behind, and master had refactored the exact constructors this PR touches, so this was not a mechanical replay:

• SciMLFunctionOptions now replaces the GeneratedFunctionOptions pattern, and ODEFunction/NonlinearFunction are each split into a kwargs... wrapper plus an opts-based inner method. The has_any_limited guards moved into the inner methods beside check_complete, and nlstep_limit_bounds is threaded through both signatures and the forwarding call.
• Versions are master's, minor-bumped for the new public API: ModelingToolkit 11.39.0, ModelingToolkitBase 1.62.0.
• The six commits were squashed into one before rebasing. That was deliberate rather than cosmetic: the intermediate commits contained designs that later commits reverse (the hooks began as NonlinearFunction fields and moved to solve keywords), so replaying them through the drift would have meant resolving conflicts in code the next commit deletes. The squashed message records the full design arc.

Dependencies are all released now, so this PR no longer depends on anything unmerged: SciMLBase 3.43 (#1449), NonlinearSolveBase 2.41 (SciML/NonlinearSolve.jl#1084, #1142), and SciML/OrdinaryDiffEq.jl#4158 for the stage-solve forwarding.

Verified locally after the rebase (Julia 1.12.4):

Run	Result
limited_operator.jl with ModelingToolkit visible (CI's group_interfaceii.jl path)	96/96
limited_operator.jl standalone ModelingToolkitBase	66/66
homotopy_lowering.jl (shares the touched lowering pipeline)	43/43

The MTK-visible run is the one that matters for this rebase: the nlstep testsets are gated on @isdefined(ModelingToolkit), and they are what exercise the nlstep_limit_bounds plumbing I re-threaded by hand. Running the file under @safetestset hides that gate — its module isolation makes an outer using ModelingToolkit invisible, silently dropping those 30 assertions — so it was run without the wrapper.

🤖 Generated with Claude Code

[AayushSabharwal]

Also needs a rebase

[AayushSabharwal]

This still doesn't use IRStructure. It will not scale, I had to refactor the homotopy operator to make Multibody models stop stalling. Use the new implementation of the homotopy operator to implement this. It has a very similar template of functions.

[AayushSabharwal]

Do we want limitnew/limitold to be namespaced? If we don't, it can easily be avoided by making them GlobalScoped.

[AayushSabharwal]

This really shouldn't build and call the substituter every time. There should be a method that returns the relevant IRSubstituter and it should then be reused.

[AayushSabharwal]

IRSubstituter should be used here.

[AayushSabharwal]

This should just be a normal function instead of a singleton callable struct.

[AayushSabharwal]

Should use the IRStructure method of query.

[ChrisRackauckas-Claude]

@AayushSabharwal — rebased onto bf1efdb992 and the review is addressed. Per-thread mapping:

Comment	Change
"This still doesn't use IRStructure… Use the new implementation of the homotopy operator"	has_limited now has the two-method shape has_homotopy has: a convenience method that builds an IRStructure, and has_limited(ir, expr) doing SU.query(__has_limited_predicate, ir, expr). has_any_limited walks equations and observed through the system's own get_irstructure(sys), so the DAG is visited once. This one mattered most: has_any_limited runs on every ODEFunction/NonlinearFunction construction, so the naive recursion was on the hot path for every user, not just limited models.
"Should use the IRStructure method of query"	Same change.
"This should just be a normal function instead of a singleton callable struct"	IsLimitedNode is now __has_limited_predicate, a plain function, used both as the query predicate and as the is_atomic for search_variables!. IsSentinel is gone entirely (below).
"Do we want limitnew/limitold to be namespaced? … can easily be avoided by making them GlobalScoped"	Done, and it is a strict deletion: GlobalScope means a limiter written inside a component still refers to the two toplevel symbols after flattening, so _sentinel_kind (name-suffix matching), IsSentinel and _canonicalize_sentinels all disappear and _remove_sentinels becomes an identity test.
"This really shouldn't build and call the substituter every time. There should be a method that returns the relevant IRSubstituter and it should then be reused."	The per-call substituter this pointed at was _canonicalize_sentinels, which no longer exists. In attach_stage_limiters, the stage substitution now builds one IRSubstituter before the loop and reuses it across every limited quantity.
"IRSubstituter should be used here" (attach_stage_limiters)	The stage rewrite goes through that shared substituter; _conjugate_limiter takes the IRStructure and uses IRSubstituter too. The two remaining substitute calls there were the derivative-plus-substitution pair extracting a and b from q = a*z + b — those are now a single cached Symbolics.LinearExpander via get_linear_expander_for!(stagesys, z, true), which also makes the affine check exact (strict, so a z inside an ifelse condition is not called affine).
"Do we want to FunctionWrapper these somehow to cut down on compile times?"	Fixed, but not with a FunctionWrapper — see below.

On the last one. A FunctionWrapper does not fit: the limiter is called as fn((u_new, u_old), p) where the value type has to stay generic (ForwardDiff Duals flow through those entries when the outer problem is differentiated) and p is an MTKParameters whose type varies per system, so a concrete wrapper would either break AD or need one wrapper type per system anyway.

What was actually costing compile time is that there was one generated function per limiter, held in a heterogeneous tuple. generate_limited_postcondition now composes correctors symbolically — several correctors on the same unknown fold into one expression before codegen, which is where the clamp ∘ limiter ordering now lives — and emits one generate_custom_function call for the whole registry. LimitedPostcondition is {iip, F} over a single function plus a Vector{Int} of indices: concretely typed, and compile time independent of how many quantities a device library limits.

---

Separately, the first full CI run exposed three plumbing breakages, all fixed here (details and the failing-before evidence are in the PR body):

• tests / InterfaceII errored with Package NonlinearSolveBase not found — the test file is registered in the root suite too, but the root Project.toml never listed it.
• downgrade-mtkbase (InterfaceII) was resolving NonlinearSolveBase 2.38.1 and NonlinearSolveFirstOrder 2.2.4, neither of which knows about postcondition, so limiting was silently inert: 179 Newton steps instead of 11, and the solve-time override never called. get_postcondition first exists in NonlinearSolveBase 2.41 and supports_postcondition in NonlinearSolveFirstOrder 2.3; the latter is not a declared dependency, so the floor has to be carried by NonlinearSolve 4.29.1, which pins NonlinearSolveBase to ≥ 2.45. Verified by pinning an env to exactly those floors and re-running the file: 66/66.
• Aqua / all_qualified_accesses_are_public flagged Base.setindex; it joins the existing Base entries in the ignore list.

The one QA failure that remains locally (JET report_package, 234 reports, all Moshi @match temporaries) reproduces identically on master at the rebase base and is tracked in #5063 — not from this PR.

🤖 Generated with Claude Code (model: claude-opus-5[1m])

https://claude.ai/code/session_01T2Y8c85c7AiQHEjK9KLzdA