Add the limited(actual, limiter) operator: symbolic PCNR iterate limiting for nonlinear solves

Project.toml

@@ -1,6 +1,6 @@
 name = "ModelingToolkit"
 uuid = "961ee093-0014-501f-94e3-6117800e7a78"
-version = "11.38.2"
+version = "11.39.0"
 authors = ["Yingbo Ma <mayingbo5@gmail.com>", "Chris Rackauckas <accounts@chrisrackauckas.com> and contributors"]
 
 [deps]
@@ -87,11 +87,12 @@ Libdl = "1"
 LinearAlgebra = "1"
 LinearSolve = "3.66, 4, 5"
 Logging = "1"
-ModelingToolkitBase = "1.57"
+ModelingToolkitBase = "1.62"
 ModelingToolkitStandardLibrary = "2.20"
 ModelingToolkitTearing = "1.19.2"
 Moshi = "0.3.6"
-NonlinearSolve = "4.3"
+NonlinearSolve = "4.25"
+NonlinearSolveBase = "2.41"
 OffsetArrays = "1.11"
 Optimization = "5.5"
 OptimizationBase = "5.1.0"
@@ -112,7 +113,7 @@ REPL = "1"
 Reexport = "1"
 RuntimeGeneratedFunctions = "0.5.12"
 SCCNonlinearSolve = "1.13"
-SciMLBase = "3.19"
+SciMLBase = "3.37"
 SciMLPublic = "1.0.0"
 Serialization = "1"
 Setfield = "1"
@@ -169,6 +170,7 @@ LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
 Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
 ModelingToolkitStandardLibrary = "16a59e39-deab-5bd0-87e4-056b12336739"
 NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
+NonlinearSolveBase = "be0214bd-f91f-a760-ac4e-3421ce2b2da0"
 Optimization = "7f7a1694-90dd-40f0-9382-eb1efda571ba"
 OptimizationBase = "bca83a33-5cc9-4baa-983d-23429ab6bcbb"
 OptimizationOptimJL = "36348300-93cb-4f02-beb5-3c3902f8871e"
@@ -199,4 +201,4 @@ TestEnv = "1e6cf692-eddd-4d53-88a5-2d735e33781b"
 URIs = "5c2747f8-b7ea-4ff2-ba2e-563bfd36b1d4"
 
 [targets]
-test = ["BenchmarkTools", "BoundaryValueDiffEqMIRK", "BoundaryValueDiffEqAscher", "ControlSystemsBase", "DataInterpolations", "DelayDiffEq", "NonlinearSolve", "ForwardDiff", "ModelingToolkitStandardLibrary", "Optimization", "OptimizationOptimJL", "OrdinaryDiffEq", "OrdinaryDiffEqCore", "OrdinaryDiffEqDefault", "REPL", "Random", "ReferenceTests", "SafeTestsets", "SciMLTesting", "StableRNGs", "Statistics", "SteadyStateDiffEq", "Test", "StochasticDiffEq", "Sundials", "OrdinaryDiffEqNonlinearSolve", "Logging", "OptimizationBase", "LinearSolve", "Latexify", "Distributed", "DiffEqNoiseProcess", "DynamicQuantities", "DiffEqCallbacks", "URIs", "JumpProcesses", "RecursiveArrayTools", "SpecialFunctions", "SciCompDSL", "OrdinaryDiffEqRosenbrock", "OrdinaryDiffEqSDIRK", "OrdinaryDiffEqBDF", "OrdinaryDiffEqFunctionMap", "OrdinaryDiffEqFIRK"]
+test = ["BenchmarkTools", "BoundaryValueDiffEqMIRK", "BoundaryValueDiffEqAscher", "ControlSystemsBase", "DataInterpolations", "DelayDiffEq", "NonlinearSolve", "NonlinearSolveBase", "ForwardDiff", "ModelingToolkitStandardLibrary", "Optimization", "OptimizationOptimJL", "OrdinaryDiffEq", "OrdinaryDiffEqCore", "OrdinaryDiffEqDefault", "REPL", "Random", "ReferenceTests", "SafeTestsets", "SciMLTesting", "StableRNGs", "Statistics", "SteadyStateDiffEq", "Test", "StochasticDiffEq", "Sundials", "OrdinaryDiffEqNonlinearSolve", "Logging", "OptimizationBase", "LinearSolve", "Latexify", "Distributed", "DiffEqNoiseProcess", "DynamicQuantities", "DiffEqCallbacks", "URIs", "JumpProcesses", "RecursiveArrayTools", "SpecialFunctions", "SciCompDSL", "OrdinaryDiffEqRosenbrock", "OrdinaryDiffEqSDIRK", "OrdinaryDiffEqBDF", "OrdinaryDiffEqFunctionMap", "OrdinaryDiffEqFIRK"]

docs/Project.toml

@@ -26,7 +26,9 @@ Optimization = "7f7a1694-90dd-40f0-9382-eb1efda571ba"
 OptimizationOptimJL = "36348300-93cb-4f02-beb5-3c3902f8871e"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 OrdinaryDiffEqFunctionMap = "d3585ca7-f5d3-4ba6-8057-292ed1abd90f"
+OrdinaryDiffEqNonlinearSolve = "127b3ac7-2247-4354-8eb6-78cf4e7c58e8"
 OrdinaryDiffEqRosenbrock = "43230ef6-c299-4910-a778-202eb28ce4ce"
+OrdinaryDiffEqSDIRK = "2d112036-d095-4a1e-ab9a-08536f3ecdbf"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 PreallocationTools = "d236fae5-4411-538c-8e31-a6e3d9e00b46"
 SciCompDSL = "91a8cdf1-4ca6-467b-a780-87fda3fff15e"
@@ -63,7 +65,9 @@ Optimization = "5"
 OptimizationOptimJL = "0.4"
 OrdinaryDiffEq = "6.31, 7"
 OrdinaryDiffEqFunctionMap = "1, 2"
+OrdinaryDiffEqNonlinearSolve = "1.5.0, 2"
 OrdinaryDiffEqRosenbrock = "1, 2"
+OrdinaryDiffEqSDIRK = "1, 2"
 Plots = "1.36"
 PreallocationTools = "0.4, 1"
 SciCompDSL = "1"

docs/pages.jl

@@ -48,6 +48,7 @@ pages = [
         "basics/Composition.md",
         "basics/Events.md",
         "basics/Homotopy.md",
+        "basics/Limiting.md",
         "basics/Linearization.md",
         "basics/InputOutput.md",
         "basics/MTKLanguage.md",

docs/src/basics/Limiting.md

@@ -0,0 +1,207 @@
+# [Iterate Limiting: the `limited` operator](@id limiting)
+
+Nonlinear solves over device equations with exponential I-V characteristics — diodes,
+BJTs, MOSFETs — are the classic hard case for Newton's method: a volt-sized overshoot of
+a junction voltage puts `exp(v/Vt)` astronomically far from the linearization. SPICE-family
+simulators handle this with *limiting*: each Newton update of a sensitive quantity is
+clipped to a trusted move relative to its previous value. The Predictor/Corrector
+Newton-Raphson (PCNR) method of Aadithya, Keiter & Mei recasts limiting consistently by
+making each limited quantity an explicit unknown, applying the limiter as a corrector
+between Newton steps, and always evaluating residuals at the corrected iterates.
+
+ModelingToolkit exposes this symbolically through the [`limited`](@ref) operator, in the
+spirit of the Modelica-style [`homotopy`](@ref homotopy) operator: a component author
+annotates the model once, and every nonlinear solve built from the model gets
+predictor/corrector limiting automatically.
+
+```julia
+limited(actual, limiter)
+```
+
+  - `actual` — the expression being limited (e.g. a junction voltage).
+  - `limiter` — the correction rule, written in terms of the reserved placeholders
+    [`limitnew`](@ref) (the proposed value) and [`limitold`](@ref) (the previously
+    accepted value), plus parameters.
+
+## Example: a diode component with `pnjlim`
+
+Register the SPICE3 junction limiting function as an opaque symbolic function (keeping
+Julia's short-circuit branch semantics), and annotate the diode's junction voltage:
+
+```@example limiting
+using ModelingToolkit, NonlinearSolve
+
+function pnjlim(vnew, vold, vt, vcrit)
+    if vnew > vcrit && abs(vnew - vold) > 2vt
+        if vold > 0
+            arg = 1 + (vnew - vold) / vt
+            vnew = arg > 0 ? vold + vt * log(arg) : vcrit
+        else
+            vnew = vt * log(vnew / vt)
+        end
+    end
+    return vnew
+end
+@register_symbolic pnjlim(vnew, vold, vt, vcrit)
+
+function DCDiode(; name, Is = 1.0e-14, Vt = 0.025)
+    @variables v i
+    ps = @parameters begin
+        (Is::Float64 = Is)
+        (Vt::Float64 = Vt)
+        (vcrit::Float64 = Vt * log(Vt / (sqrt(2) * Is)))
+    end
+    eqs = [i ~ Is * (exp(limited(v, pnjlim(limitnew, limitold, Vt, vcrit)) / Vt) - 1)]
+    return System(eqs, [v, i], ps; name)
+end
+
+function DCResistor(; name, R = 1.0e3)
+    @variables v i
+    @parameters R = R
+    return System([v ~ i * R], [v, i], [R]; name)
+end
+
+@named diode = DCDiode()
+@named res = DCResistor()
+@parameters Vs = 5.0
+connections = [res.i ~ diode.i, Vs ~ res.v + diode.v]
+@named circuit = System(connections, [], [Vs]; systems = [diode, res])
+csys = mtkcompile(circuit)
+
+prob = NonlinearProblem(csys, [diode.v => 0.0, res.i => 0.0])
+sol = solve(prob, NewtonRaphson())
+sol[diode.v], sol[diode.i], sol.stats.nsteps
+```
+
+Without the annotation, plain Newton needs a couple hundred millivolt-creep iterations on
+this circuit; with it, the solve converges in about a dozen. Nothing about the *usage*
+changed — the limiting behavior travels with the component definition.
+
+## What `mtkcompile` does with `limited`
+
+For **time-independent** systems, each unique `limited(actual, limiter)` node is lowered
+in the PCNR augmented form:
+
+ 1. an auxiliary irreducible unknown `limited_k` is introduced and the node is replaced
+    by it (irreducible, so structural simplification keeps the limited quantity as the
+    surviving representative of its alias class — the reduction of the augmented system
+    happens symbolically);
+ 2. the consistency equation `limited_k ~ actual` is appended, and `limited_k` receives
+    the symbolic guess `actual`;
+ 3. the limiters are compiled into a `postcondition` corrector hook on the generated
+    `SciMLBase.NonlinearFunction`, which NonlinearSolve.jl's native solvers apply to every
+    accepted iterate before evaluating the residual there.
+
+Solving therefore requires a solver that supports `postcondition` (e.g. `NewtonRaphson`
+and the other native NonlinearSolve.jl methods); unsupported solvers throw instead of
+silently ignoring the limiter.
+
+For **time-dependent** systems the operator is stripped to `actual` during `mtkcompile`,
+so the same component library compiles unchanged for transient simulation. The stripped
+limiters are recorded, and come back in the stage solves of an implicit step — see below.
+
+## Limiting inside a transient solve
+
+The right-hand side of an ODE is not where transient limiting belongs: it is the *nonlinear
+solve of an implicit step* whose Newton iterates need damping. ModelingToolkit can hand an
+implicit solver a symbolically built nonlinear system for exactly those stage equations —
+`M*z = outer_tmp + γ₁*f(γ₂*z + inner_tmp, p, c)` — by constructing the problem with
+`nlstep = true`. When the model declares limited quantities, the limiters are re-attached to
+that stage system as its `postcondition`:
+
+```@example limiting
+using ModelingToolkit: t_nounits as t, D_nounits as D
+using OrdinaryDiffEqSDIRK
+import OrdinaryDiffEqNonlinearSolve
+
+@variables v(t) = 0.0
+@parameters Is = 1e-14 Vt = 0.025 C = 1e-6 R = 1e3 Vsrc = 5.0
+vcrit = 0.71
+rc = [D(v) ~ ((Vsrc - v) / R -
+              Is * (exp(limited(v, pnjlim(limitnew, limitold, Vt, vcrit)) / Vt) - 1)) / C]
+@mtkcompile rcsys = System(rc, t)
+
+prob = ODEProblem(rcsys, [], (0.0, 1e-4); nlstep = true)
+alg = ImplicitEuler(nlsolve = OrdinaryDiffEqNonlinearSolve.NonlinearSolveAlg())
+sol = solve(prob, alg; dt = 1e-6, adaptive = false)
+sol[v, end]
+```
+
+Two things are worth being explicit about.
+
+**`limitold` is the previous Newton iterate, not the previous time step.** The corrector
+runs inside the stage solve, so it compares each proposed value of a limited quantity
+against the previous iterate of that same solve (which starts from the step's predictor).
+That is exactly what SPICE-style limiting means: it damps the iteration, not the
+trajectory. The converged step is unaffected — a limiter is the identity at a fixed point,
+and the residual is always evaluated after the correction — so limiting changes how the
+step is found, never what it is.
+
+**The limited quantity must be an affine function of one stage unknown.** The stage system's
+unknowns have to remain the ODE unknowns (the solver maps between them by index), so unlike
+the standalone nonlinear case there is no room for an auxiliary unknown per limited
+quantity. Instead the correction on the quantity `q = a*z + b` is applied to its stage
+unknown `z` as the conjugated limiter `(L(a*znew + b, a*zold + b) - b) / a`. A limited
+quantity that resolves to several stage unknowns, or nonlinearly to one, is an error at
+problem construction rather than a silently dropped limiter; build with `nlstep = false` to
+compile such a model without limiting. `nlstep_scc = true` is likewise rejected, since the
+SCC decomposition splits the stage problem the corrector is attached to.
+
+## `bounds` metadata limits the stage iterates too
+
+A stage unknown is the Newton *increment* `z`, not the state: the physical value is
+`γ₂*z + inner_tmp`, and `γ₂` and `inner_tmp` are rewritten by the stepper every stage. A
+state's `bounds` metadata is therefore not a box on anything the stage solver iterates on,
+and `nlstep` never hands the stage problem a static `lb`/`ub`.
+
+What it does instead is deliver the box the same way a limiter is delivered. `bounds` on a
+state becomes the limiter `clamp(limitnew, lo, hi)` of that state, which goes through the
+same conjugation as everything else and lands on the stage unknown as
+`(clamp(a*znew + b, lo, hi) - b) / a`. Because `a` and `b` are expressions in the stage
+parameters, the compiled corrector reads them live and enforces the *correct* box on the
+increment at every stage — which a static vector could not do.
+
+```@example limiting
+@variables c(t) = 0.5 [bounds = (0.0, 1.0)]
+@parameters k = 1.0
+@mtkcompile csys2 = System([D(c) ~ -k * c^2 * (1 - c)], t)
+
+bprob = ODEProblem(csys2, [], (0.0, 1.0); nlstep = true)
+bprob.f.nlstep_data.nlprob.lb, haskey(bprob.f.nlstep_data.nlprob.kwargs, :postcondition)
+```
+
+This is on by default, because it is a safety property rather than a modelling choice:
+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` or
+`exp` from being evaluated outside its domain — which would otherwise `NaN` the entire
+solve rather than merely fail a step. Pass `nlstep_limit_bounds = false` to opt out, for a
+model whose `bounds` are an advisory range that its stage solutions legitimately leave.
+
+When a state carries both a `limited(...)` annotation and `bounds`, the two compose in that
+order: the model's limiter damps the proposed move, the clamp then projects whatever it
+produced, so the corrected iterate always ends up in range. Unlike a declared limiter, a
+box that cannot be conjugated — a bounded state that the stage system tears away as a
+nonlinear function of other unknowns — is dropped rather than raising, since it is a safety
+net and not a request; the same is true under `nlstep_scc = true`.
+
+## Contracts
+
+  - `limiter` must satisfy `limiter == limitnew` when `limitnew == limitold`, so
+    solutions are fixed points of the correction.
+  - `limiter` may reference `limitnew`, `limitold`, and parameters (including bound
+    parameters); referencing other unknowns is an error.
+  - `limited` nodes may not be nested.
+  - Jacobians treat `limited(actual, limiter)` as `actual`: the limiter is a corrector
+    between iterations, not part of the residual.
+
+## Docstrings
+
+```@docs
+ModelingToolkitBase.limited
+ModelingToolkitBase.limitnew
+ModelingToolkitBase.limitold
+```
+
+**Reference:** K. V. Aadithya, E. R. Keiter, T. Mei, *Predictor/Corrector Newton-Raphson
+(PCNR): A Simple, Flexible, Scalable, Modular, and Consistent Replacement for Limiting in
+Circuit Simulation*, Scientific Computing in Electrical Engineering, 2020.

lib/ModelingToolkitBase/Project.toml

@@ -1,7 +1,8 @@
 name = "ModelingToolkitBase"
 uuid = "7771a370-6774-4173-bd38-47e70ca0b839"
+version = "1.55.0"
 authors = ["Yingbo Ma <mayingbo5@gmail.com>", "Chris Rackauckas <accounts@chrisrackauckas.com> and contributors"]
-version = "1.61.0"
+version = "1.62.0"
 
 [deps]
 ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b"
@@ -102,6 +103,8 @@ BenchmarkTools = "1"
 BifurcationKit = "0.4, 0.5, 0.8"
 BipartiteGraphs = "0.1.3"
 BlockArrays = "1.1"
+BoundaryValueDiffEqAscher = "1.6.0"
+BoundaryValueDiffEqMIRK = "1.7.0"
 CasADi = "1.0.7"
 ChainRulesCore = "1"
 Combinatorics = "1"

lib/ModelingToolkitBase/src/ModelingToolkitBase.jl

@@ -213,6 +213,9 @@ include("systems/problem_utils.jl")
 # Operator + lowering layer; must load before the problem constructors that
 # consume it (problems/nonlinearproblem.jl selector + problems/homotopyproblem.jl).
 include("systems/homotopy_operator.jl")
+# The `limited` operator: symbolic SPICE-style iterate limiting (PCNR), lowered during
+# `mtkcompile` (`apply_limited_lowering`) and consumed by `problems/nonlinearproblem.jl`.
+include("systems/limited_operator.jl")
 
 include("problems/compatibility.jl")
 include("problems/odeproblem.jl")
@@ -357,6 +360,7 @@ export alg_equations, diff_equations, has_alg_equations, has_diff_equations
 export get_alg_eqs, get_diff_eqs, has_alg_eqs, has_diff_eqs
 
 export homotopy
+export limited, limitnew, limitold
 
 export @variables, @parameters, @independent_variables, @constants, @brownians, @brownian,
     @poissonians, @discretes