Modelica homotopy(actual, simplified) operator for initialization (spec 3.7.4)

NEWS.md

@@ -1,3 +1,22 @@
+# Unreleased
+
+## New features
+
+  - Added support for Modelica's `homotopy(actual, simplified)` operator for
+    initialization (Modelica spec 3.7.4.2). The operator stays opaque through
+    `mtkcompile` and evaluates to `actual` at runtime; when the initialization
+    problem is built for a system containing `homotopy` calls, it is lowered to
+    a `SciMLBase.HomotopyProblem` whose residual `f(u, p, λ)` blends
+    `(1 - λ)*simplified + λ*actual`, solved by the
+    `NonlinearSolveBase.HomotopySweep` continuation solver — a
+    natural-parameter sweep from `λ=0` (the easy-to-solve `simplified` form) to
+    `λ=1` (`actual`), injected as the default `initializealg` (an explicit
+    user-passed `initializealg` takes precedence). This makes initialization of
+    pressure-driven-flow and similar nonlinear systems robust to poor guesses
+    by starting from the simplified form. `NonlinearSolveBase` is now a direct
+    dependency of `ModelingToolkitBase`. Systems that do not use the operator
+    are unaffected.
+
 # ModelingToolkit v11 Release Notes
 
 ## Symbolics@7 and SymbolicUtils@4 compatibility

docs/pages.jl

@@ -47,6 +47,7 @@ pages = [
     "Basics" => Any[
         "basics/Composition.md",
         "basics/Events.md",
+        "basics/Homotopy.md",
         "basics/Linearization.md",
         "basics/InputOutput.md",
         "basics/MTKLanguage.md",

docs/src/basics/Homotopy.md

@@ -0,0 +1,167 @@
+# [Homotopy Initialization](@id homotopy)
+
+ModelingToolkit implements the Modelica `homotopy(actual, simplified)` operator
+([Modelica Specification 3.7.4.2](https://specification.modelica.org/master/operators-and-expressions.html#homotopy))
+as an initialization aid for systems with hard-to-initialize algebraic equations.
+
+## What Is Homotopy Initialization?
+
+Many physical models (pressure-driven pipe networks, power-flow equations, chemical
+equilibria) have algebraic constraints that are easy to solve under simplified
+assumptions but hard to solve from a cold start. The `homotopy` operator encodes
+both the actual expression and a simplified approximation in a single annotation:
+
+```julia
+homotopy(actual, simplified)
+```
+
+During initialization the solver may start from the easy `simplified` form and
+continuously deform it into `actual`; outside initialization the operator
+evaluates to `actual`, as the Modelica specification prescribes.
+
+## Runtime Semantics
+
+The operator stays an opaque symbolic function through `System` construction,
+`mtkcompile`, and runtime code generation — no continuation parameter is added to
+the system. At runtime the generated code calls the numeric fallback
+
+```julia
+homotopy(actual::Real, simplified::Real) = actual
+```
+
+so the operator evaluates to `actual`. Note the honest cost: the `simplified`
+argument expression is still evaluated and its value discarded (arguments are
+evaluated before the call). This small overhead is borne only by systems that use
+`homotopy` — systems without the operator go through a byte-identical pipeline and
+are completely unaffected.
+
+## Initialization Semantics
+
+When the initialization problem is built (after `mtkcompile`), the compiled
+initialization system is scanned for `homotopy` nodes in its equations and
+observed equations. If any are found and the initialization problem is a square
+`NonlinearProblem` with unknowns, the initialization residual is regenerated with
+every `homotopy(actual, simplified)` replaced by the blend
+
+```math
+(1 - \lambda)\,\text{simplified} + \lambda\,\text{actual}
+```
+
+and compiled as `f(u, p, λ)`, where `λ` is an explicit trailing argument of the
+residual — it is never added to the system's parameters, and your parameter object
+`p` passes through untouched. All `homotopy` calls in a system share the single
+`λ`, per the Modelica spec's recommendation of (conceptually) one homotopy
+iteration over the whole model; this includes nested `homotopy` calls. The result is wrapped as a `SciMLBase.HomotopyProblem`
+(its `λspan` defaults to `(0.0, 1.0)`).
+
+Problems built through the default path then receive
+
+```julia
+initializealg = SciMLBase.OverrideInit(nlsolve = NonlinearSolveBase.HomotopySweep())
+```
+
+as their default initialization algorithm. `HomotopySweep` performs a
+natural-parameter continuation: it sweeps `λ` from `0` (the easy `simplified`
+equations) to `1` (`actual`), solving a standard nonlinear problem at each step and
+warm-starting from the previous step's solution. `NonlinearSolveBase` is a hard
+dependency of `ModelingToolkitBase`, so no extra package load is needed to
+activate the sweep.
+
+Symbolic differentiation works through the operator: nodewise derivative rules
+keep symbolic jacobians, `tgrad`, and index reduction consistent. At runtime,
+differentiated equations reproduce `actual`'s derivative; in the swept residual
+they follow the derivative of the blended expression above.
+
+## When Does the Sweep Activate?
+
+  - A system whose compiled initialization system contains `homotopy` nodes (in
+    equations or observed equations) gets a swept initialization, as described
+    above.
+  - A system without `homotopy` is untouched: no detection cost beyond the scan,
+    no problem rewrapping, no `initializealg` injection.
+  - An explicit user-passed `initializealg` always wins — the default is only
+    injected when you did not pass one, and a solve-time `initializealg` overrides
+    the problem default as usual.
+  - If simplification eliminates every `homotopy` node before the initialization
+    system is built, initialization proceeds as standard (detection happens on the
+    compiled initialization system).
+  - Multiple (or nested) `homotopy` calls in one system are blended with the one
+    shared `λ` and swept together in a single continuation.
+  - If the initialization would otherwise target an `SCCNonlinearProblem`, a
+    homotopy-using system is forced to the monolithic `NonlinearProblem` target,
+    since the continuation deforms the coupled system along one `λ`.
+
+## Example: Out-of-Basin Rescue
+
+The equation `0 = atan(y - 3)` has a root at `y = 3`, but a Newton solver starting
+from `y = 12` diverges because `atan` saturates. Using `homotopy` with
+`simplified = y` (whose root is `y = 0`) lets the continuation walk from the easy
+root to the true one:
+
+```julia
+using ModelingToolkit, NonlinearSolve, OrdinaryDiffEqRosenbrock
+using ModelingToolkit: t_nounits as t, D_nounits as D
+
+@variables x(t) y(t)
+@named sys = System([D(x) ~ -x, 0 ~ homotopy(atan(y - 3), y)], t; guesses = [y => 12.0])
+sys = mtkcompile(sys)
+prob = ODEProblem(sys, [x => 1.0], (0.0, 1.0))
+integ = init(prob, Rodas5P())
+integ[y] # ≈ 3.0 — the sweep rescued the out-of-basin guess
+```
+
+The `guesses` keyword provides the starting point of the continuation (here
+`y = 12`); the sweep deforms the equations so the solver reaches `y ≈ 3` at
+`λ = 1`.
+
+## Customizing the Continuation Solver
+
+To tune the sweep, pass your own `initializealg` (which always takes precedence
+over the injected default):
+
+```julia
+integ = init(prob, Rodas5P();
+    initializealg = SciMLBase.OverrideInit(
+        nlsolve = NonlinearSolveBase.HomotopySweep(nsteps = 30)))
+```
+
+`HomotopySweep` accepts the keyword arguments `inner` (the nonlinear algorithm
+used at each step), `nsteps`, `adaptive`, `initial_step_factor`, and `min_dλ`; see
+the `NonlinearSolveBase.HomotopySweep` docstring for their meanings and defaults.
+Tolerances are inherited from the integrator's `abstol`/`reltol`.
+
+## Limitations
+
+  - **Non-square initialization cannot be swept.** If the initialization problem
+    is not a square `NonlinearProblem` (e.g. over-determined initialization, which
+    falls back to least squares), a one-time warning is emitted and initialization
+    solves the `actual` (`λ = 1`) form instead. Least-squares and per-SCC sweeps
+    are future work.
+  - **Trivial initialization is silently correct.** If the initialization system
+    has no unknowns, there is nothing to solve and no sweep happens; observed
+    evaluation folds `homotopy` to `actual`, the correct
+    outside-initialization semantics.
+  - **`expression = Val{true}` does not auto-inject the `initializealg`.** When
+    generating problem-construction expressions, pass
+    `initializealg = SciMLBase.OverrideInit(nlsolve = NonlinearSolveBase.HomotopySweep())`
+    explicitly.
+  - **Parameter bindings and defaults are not rewritten.** A `homotopy` call
+    inside a parameter binding or default value evaluates as `actual` at all `λ`.
+  - **Scalar `Real` expressions only**, matching Modelica's restriction; the
+    numeric fallback method is `homotopy(::Real, ::Real)`.
+  - **Failure surfaces as a failed initialization.** If the continuation cannot
+    reach `λ = 1`, the sweep returns a failure return code and initialization
+    fails (see `NonlinearSolveBase.HomotopySweep` for the failure contract).
+
+## API Reference
+
+```@docs
+ModelingToolkit.homotopy
+```
+
+## See Also
+
+  - [Initialization of Systems](@ref initialization) — general DAE initialization,
+    bindings, guesses, and `OverrideInit`.
+  - [Modelica Specification 3.7.4.2](https://specification.modelica.org/master/operators-and-expressions.html#homotopy)
+    — the upstream specification this operator implements.

lib/ModelingToolkitBase/Project.toml

@@ -36,6 +36,7 @@ Libdl = "8f399da3-3557-5675-b5ff-fb832c97cbdb"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Moshi = "2e0e35c7-a2e4-4343-998d-7ef72827ed2d"
 NaNMath = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
+NonlinearSolveBase = "be0214bd-f91f-a760-ac4e-3421ce2b2da0"
 OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
 OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
 PreallocationTools = "d236fae5-4411-538c-8e31-a6e3d9e00b46"
@@ -143,7 +144,8 @@ ModelingToolkitStandardLibrary = "2.20"
 Mooncake = "0.5"
 Moshi = "0.3.6"
 NaNMath = "0.3, 1"
-NonlinearSolve = "4.3"
+NonlinearSolve = "4.19"
+NonlinearSolveBase = "2.31"
 OffsetArrays = "1"
 OptimizationBase = "5"
 OptimizationIpopt = "1"
@@ -169,12 +171,12 @@ ReferenceTests = "0.10"
 RuntimeGeneratedFunctions = "0.5.12"
 SCCNonlinearSolve = "1.8.1"
 SafeTestsets = "0.1"
-SciMLBase = "3.18"
+SciMLBase = "3.19"
 SciMLPublic = "1.0.0"
 SciMLStructures = "1.7"
 Serialization = "1"
 Setfield = "0.7, 0.8, 1"
-SimpleNonlinearSolve = "0.1.0, 1, 2"
+SimpleNonlinearSolve = "2.11.1"
 SparseArrays = "1"
 SpecialFunctions = "1, 2"
 StableRNGs = "1"

lib/ModelingToolkitBase/src/ModelingToolkitBase.jl

@@ -101,6 +101,8 @@ import Random: AbstractRNG
 import DomainSets
 # For `LinearInitializationProblem`
 import SCCNonlinearSolve
+# For the default `homotopy(...)` continuation solver (`HomotopySweep`)
+import NonlinearSolveBase
 using TaskLocalValues: TaskLocalValue
 
 export @derivatives
@@ -214,6 +216,7 @@ include("modelingtoolkitize/sdeproblem.jl")
 include("modelingtoolkitize/optimizationproblem.jl")
 include("modelingtoolkitize/nonlinearproblem.jl")
 
+include("systems/homotopy.jl")
 include("systems/nonlinear/homotopy_continuation.jl")
 include("systems/nonlinear/initializesystem.jl")
 include("systems/diffeqs/basic_transformations.jl")
@@ -300,6 +303,8 @@ export generate_initializesystem, Initial, isinitial, InitializationProblem
 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 @variables, @parameters, @independent_variables, @constants, @brownians, @brownian,
     @poissonians, @discretes
 export @named, @nonamespace, @namespace, extend, compose, complete, toggle_namespacing

lib/ModelingToolkitBase/src/problems/initializationproblem.jl

@@ -123,7 +123,98 @@ All other keyword arguments are forwarded to the wrapped problem constructor.
         sys, isys; warn_initialize_determined,
         kwargs...
     )
-    TProb{_iip}(isys, op; kwargs..., build_initializeprob = false, is_initializeprob = true)
+    # Modelica homotopy(actual, simplified) (spec 3.7.4.2): scan the compiled init system
+    # (equations ∪ observed) for homotopy nodes. If any are present, the square init
+    # problem is rebuilt by `_maybe_sweep_homotopy_initprob` below as a
+    # `SciMLBase.HomotopyProblem` whose residual `f(u, p, λ)` sweeps
+    # `(1 - λ)*simplified + λ*actual`, so a continuation solver can walk λ from 0 to 1.
+    homotopic = has_any_homotopy(isys)
+    if homotopic && TProb === SCCNonlinearProblem
+        # A continuation sweep deforms the COUPLED system along one λ; per-SCC
+        # sweeping is a separate future improvement. Root ModelingToolkit's
+        # `get_initialization_problem_type` prefers `SCCNonlinearProblem` for split
+        # fully-determined systems — force the monolithic target here.
+        TProb = NonlinearProblem
+    end
+    initprob = TProb{_iip}(isys, op; kwargs..., build_initializeprob = false, is_initializeprob = true)
+    if homotopic
+        initprob = _maybe_sweep_homotopy_initprob(initprob, isys, _iip, kwargs)
+    end
+    return initprob
+end
+
+# Codegen keyword arguments which `_maybe_sweep_homotopy_initprob` threads from
+# `InitializationProblem`'s kwargs into `generate_homotopy_residual`, so the swept
+# residual is compiled with the same options as the standard one. Keep in sync with
+# the codegen-relevant keywords of `SciMLBase.NonlinearFunction{iip, spec}(::System)`
+# (problems/nonlinearproblem.jl), which the standard init residual build receives.
+const HOMOTOPY_INIT_CODEGEN_KWARGS = (
+    :eval_expression, :eval_module, :checkbounds, :cse, :compiler_options,
+)
+
+"""
+    $TYPEDSIGNATURES
+
+Independent post-build pass for Modelica `homotopy(actual, simplified)` (spec 3.7.4.2),
+called from `InitializationProblem` when the compiled init system `isys` contains
+homotopy nodes. If `initprob` is a square `NonlinearProblem` with unknowns, rebuild
+it as a `SciMLBase.HomotopyProblem` whose residual `f(u, p, λ)` sweeps
+`(1 - λ)*simplified + λ*actual`, so a continuation solver can walk λ from 0 to 1.
+Otherwise return `initprob` unchanged — silently when initialization is trivial
+(no solve happens, so nothing is lost), with a one-time warning when a real init
+solve will run on the `actual` (λ = 1) form instead of being swept.
+
+`iip` is the resolved in-place flag of `initprob`. `kwargs` are
+`InitializationProblem`'s keyword arguments, from which the codegen options listed
+in `HOMOTOPY_INIT_CODEGEN_KWARGS` are threaded into the swept-residual build.
+"""
+function _maybe_sweep_homotopy_initprob(initprob, isys::AbstractSystem, iip::Bool, kwargs)
+    if initprob isa NonlinearProblem && state_values(initprob) !== nothing
+        shadow, λ = lower_homotopy_for_init(isys)
+        codegen_kwargs = (;
+            (
+                kw => kwargs[kw] for kw in HOMOTOPY_INIT_CODEGEN_KWARGS
+                    if haskey(kwargs, kw)
+            )...,
+        )
+        hf = generate_homotopy_residual(shadow, λ; codegen_kwargs...)
+        nf = initprob.f
+        swept_f = SciMLBase.NonlinearFunction{iip}(
+            hf; sys = nf.sys, observed = nf.observed,
+            resid_prototype = nf.resid_prototype
+        )
+        # Deliberately NOT carried over: jac / jac_prototype / sparsity — they
+        # encode the λ = 1 (opaque-actual) system and would be wrong mid-sweep.
+        # (The standard build did generate that λ = 1 residual — and a jac, if the
+        # user requested one — only for them to be replaced here: a deliberate
+        # layering cost on the homotopy path; the un-swept path is unaffected.)
+        # λspan is deliberately not passed: the `SciMLBase.HomotopyProblem`
+        # default is the single source of truth.
+        return SciMLBase.HomotopyProblem{iip}(
+            swept_f, initprob.u0, initprob.p; initprob.kwargs...
+        )
+    elseif initprob isa NonlinearProblem
+        # No unknowns (`state_values === nothing`): initialization is trivial —
+        # the init solve is skipped entirely, so there is no system to sweep and
+        # no homotopy intent is discarded. Downstream observed evaluation folds
+        # `homotopy(actual, simplified)` to `actual`, which is exactly the
+        # Modelica outside-initialization semantics. Correct result; no warning.
+        @debug "`homotopy(actual, simplified)` present, but initialization is " *
+            "trivial (no unknowns); nothing to sweep."
+    else
+        # Remaining legs: the NLLS fallback for non-square/singular init systems,
+        # the linear-init path (homotopy nodes without unknowns; its build product
+        # is SCC-wrapped), or an `Expr` from `expression = Val{true}`. Any future
+        # target type lands here too — warn + λ = 1 solve is the safe default. In
+        # all of them the generated code's numeric fallback folds homotopy to
+        # `actual`.
+        @warn "This system contains `homotopy(actual, simplified)`, but its " *
+            "initialization problem (a `$(nameof(typeof(initprob)))`) is not a " *
+            "square `NonlinearProblem` and cannot be swept by continuation; " *
+            "initialization proceeds with the `actual` (λ = 1) form. A " *
+            "least-squares/SCC homotopy sweep is planned as a follow-up." maxlog = 1
+    end
+    return initprob
 end
 
 function overdetermined_initialization_message(neqs::Integer, nunknown::Integer, extra::AbstractString)

lib/ModelingToolkitBase/src/systems/codegen.jl

@@ -81,6 +81,7 @@ function generate_rhs(
             rhss = [_iszero(eq.lhs) ? eq.rhs : eq.rhs - eq.lhs for eq in eqs]
         end
     else
+        # NOTE: generate_homotopy_residual in homotopy.jl mirrors the time-independent path here; keep in sync.
         if !override_discrete && !is_discrete_system(sys)
             check_operator_variables(eqs, Differential)
             check_lhs(eqs, Differential, Set(dvs))
@@ -273,7 +274,7 @@ function generate_jacobian(
 end
 
 function assert_jac_length(arr::SparseMatrixCSC, I::Vector{<:Integer}, J::Vector{<:Integer})
-    @assert findnz(arr)[1:2] == (I, J)
+    return @assert findnz(arr)[1:2] == (I, J)
 end
 
 function assert_jac_length_header(sys)
@@ -1103,8 +1104,10 @@ end
 function assemble_maj(majv::Vector{U}, unknowntoid, pmapper) where {U <: MassActionJump}
     rs = [numericrstoich(maj.reactant_stoch, unknowntoid) for maj in majv]
     ns = [numericnstoich(maj.net_stoch, unknowntoid) for maj in majv]
-    return MassActionJump(rs, ns; param_mapper = pmapper, nocopy = true,
-        scale_rates = false, rescale_rates_on_update = false)
+    return MassActionJump(
+        rs, ns; param_mapper = pmapper, nocopy = true,
+        scale_rates = false, rescale_rates_on_update = false
+    )
 end
 
 function numericrstoich(mtrs::Vector{Pair{V, W}}, unknowntoid) where {V, W}
@@ -1389,8 +1392,8 @@ Base.@nospecializeinfer function build_explicit_observed_function(
             p_start + wrap_delays, length(args) - length(rps) + 1 + wrap_delays, is_split(sys),
         ),
     }(
-            oop, iip
-        )
+        oop, iip
+    )
     return (return_inplace isa Val{true} || return_inplace isa Bool && return_inplace) ? (f, f) : f
 end