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

[BatchZero]

Summary

Implements Modelica's homotopy(actual, simplified) operator (spec 3.7.4) for initialization. A system that annotates an equation with homotopy(...) can now be initialized either by directly solving the actual equation (the spec's trivial form) or by continuation that sweeps a parameter λ from 0 (the easy-to-solve simplified equation) to 1 (the actual equation) — the spec's transformation form.

The default user experience mirrors OpenModelica: try the trivial solve first, and automatically fall back to the parameter sweep if it fails. No configuration is required. Because the trivial attempt solves at λ=1 (≡ actual), the injected default is a no-op on the happy path — numerically identical to ordinary initialization — and the sweep engages only when the plain solve fails.

m_flow = homotopy(
    actual     = sign(dp) * sqrt(abs(dp)) * k,        # true turbulent law (hard at dp≈0)
    simplified = m_flow_nominal * dp / dp_nominal,    # linear nominal (always easy)
)
# default: trivial-then-sweep, zero config
sol = solve(prob, Rodas5P())

This unblocks downstream component libraries from adding homotopy() annotations to pressure-driven-flow / nonlinear-equation components and getting robust initialization out of the box.

No code changes to SciMLBase, NonlinearSolve, or OrdinaryDiffEq — the entire feature lives in ModelingToolkit, and the sweep algorithms plug into the existing OverrideInit.nlsolve slot. The only dependency-table change is internal to ModelingToolkitBase: NonlinearSolve moves from [extras] to [weakdeps] and gains a small MTKNonlinearSolveExt extension that supplies the default inner Newton solver. When NonlinearSolve is absent the feature degrades gracefully to SimpleNonlinearSolve (already a runtime dependency), so a homotopy(...) system stays buildable and solvable with no extra packages.

Progresses #3988 and #4039 from "no homotopy support" toward "both spec forms supported." RFC #1385 stays open as the design thread.

User API

Three nlsolve algorithms, all <: SciMLBase.AbstractNonlinearAlgorithm, mapping directly onto OpenModelica's flags:

MTK	OpenModelica equivalent	behavior
(default, no initializealg)	default	TrivialThenSweep: trivial first, sweep on failure
OverrideInit(nlsolve = HomotopySweep(...))	-homotopyOnFirstTry	sweep from the start
OverrideInit(nlsolve = TrivialHomotopy())	-noHomotopyOnFirstTry	trivial only, no fallback

# force sweep from the start
sol = solve(prob, Rodas5P();
    initializealg = OverrideInit(nlsolve = HomotopySweep(set_λ! = ...)))

# disable sweep entirely
sol = solve(prob, Rodas5P();
    initializealg = OverrideInit(nlsolve = TrivialHomotopy()))

An explicit initializealg passed at solve time always overrides the injected default.

Implementation

Operator + lowering — lib/ModelingToolkitBase/src/systems/homotopy.jl:

• @register_symbolic homotopy(actual, simplified) — symbolic primitive (scalar; Julia broadcast covers the vectorized case). Runtime fallback homotopy(a::Real, s::Real) = actual.
• rewrite_trivial / rewrite_with_lambda — hand-written walks over the TermInterface protocol (NOT @rule, per @AayushSabharwal's Modelica's homotopy operator #1385 feedback that per-node matcher overhead is unacceptable at scale). rewrite_with_lambda replaces each homotopy(a, s) with (1 - λ)*s + λ*a using one shared __homotopy_λ parameter (default 1.0, so λ=1 ≡ actual).
• add_homotopy_parameter(sys) — lowers homotopy(...) nodes in both equations and observed, injects __homotopy_λ into the parent system's parameter list.

Lowering happens in complete() (abstractsystem.jl), at the same lifecycle point as MTK's own add_initialization_parameters. This is the key design decision: the parent system and the derived initialization system end up sharing one identity for __homotopy_λ, so MTKParametersReconstructor can resolve λ when it appears in init-system observed expressions (e.g. when an algebraic variable carrying a homotopy() annotation is eliminated into an observed equation — the Buildings PressureDrop pattern).

Sweep algorithms — lib/ModelingToolkitBase/src/systems/homotopy_sweep.jl:

• HomotopySweep{Inner, Sched, Setter} — walks λ ∈ schedule (default 0.0:0.1:1.0), warm-starting u0 from the previous step, writing λ via a SymbolicIndexingInterface.setp handle. On a non-converging step it returns that step's unsuccessful retcode (no silent fallback); before returning it resets λ to 1.0 (actual form) so the surfaced solution's parameter vector never exposes a half-homotopied system.
• TrivialHomotopy{Inner} — single solve at λ=1.
• TrivialThenSweep{T, S} — composite: trivial first, sweep on failure, annotating sol.original.path with :trivial / :sweep_fallback.

Default inner solver via a package extension. The algorithms' default inner is NonlinearSolve.NewtonRaphson, supplied through a new MTKNonlinearSolveExt extension keyed on a [weakdeps] dependency on NonlinearSolve. The extension populates a factory Ref at load time; when NonlinearSolve is not loaded the default falls back to SimpleNonlinearSolve.SimpleNewtonRaphson (a runtime [deps]). This replaces an earlier construction-time Base.require(<hardcoded UUID>) with the idiomatic weak-dependency pattern.

Default-alg injection — when lowering fired, InitializationMetadata carries the setp handle and a ready OverrideInit(nlsolve = TrivialThenSweep(...)). process_kwargs injects it as the default initializealg into prob.kwargs, but only when the caller didn't supply one. This is wired for ODEProblem, NonlinearProblem, SteadyStateProblem, DAEProblem, and SDEProblem (each forwards f.initialization_data to process_kwargs).

Tests

GROUP=Initialization, all green (101 assertions across 5 files):

• Homotopy Operator L0 Trivial (24): operator registration, rewrite_trivial on nested / Base.ifelse-branched / broadcast expressions, ODESystem init integration, Buildings PressureDrop fixture.
• Homotopy lowering (30): rewrite_with_lambda algebra (λ=0 → simplified, λ=1 → actual, shared λ), OverrideInitData metadata round-trip, default-initializealg injection + explicit-override + non-homotopy-no-injection, and a regression guard that observed equations are homotopy-free after lowering (an eliminated variable's definition living in observed must be lowered, not left as an opaque homotopy node).
• Homotopy nlsolve algorithms (17): HomotopySweep convergence + failure retcode + λ-reset-on-failure, TrivialHomotopy, TrivialThenSweep trivial-success and sweep-fallback paths, and that the default inner is supplied by the NonlinearSolve extension.
• Homotopy integration (23): default / explicit-sweep / explicit-trivial end-to-end through ODEProblem; an out-of-basin guess that only the sweep continuation can rescue (atan(y-3) basin escape); a clean unsuccessful-retcode (not an internal error) when the trivial path genuinely cannot converge; and that the default injection reaches NonlinearProblem and SteadyStateProblem.
• Homotopy OMC parity (7): MTK init vs. OpenModelica reference on the Buildings PressureDrop fixture.

OMC numerical parity

The PressureDrop fixture (actual = sign(dp)·√|dp|·k, simplified = m_flow_nominal·dp/dp_nominal) was run through OpenModelica v1.27.0-dev. OMC's init gives m_flow(0) = √5 ≈ 2.2360679774997896. All three MTK paths — default TrivialThenSweep, explicit HomotopySweep, explicit TrivialHomotopy — match to within 1e-6. The OMC values are frozen as constants in the test (MTK CI does not depend on an OMC install).

Other local verification

• Graceful degradation without NonlinearSolve: in a clean environment with NonlinearSolve absent, a homotopy(...) ODEProblem still constructs (no construction-time throw), the default inner falls back to SimpleNewtonRaphson, and the injected TrivialThenSweep's trivial+sweep inners are both SimpleNewtonRaphson.
• Docs: the new docs/src/API/homotopy.md page builds under Documenter — all four @docs entries (homotopy, TrivialThenSweep, HomotopySweep, TrivialHomotopy) resolve and render.

Scope notes for reviewers

• End-to-end solve coverage is exercised through the ODESystem InitializationProblem path (including the PressureDrop and out-of-basin fixtures). The default-alg injection is now also wired for NonlinearProblem / SteadyStateProblem / DAEProblem / SDEProblem and asserted structurally, but the deep numerical exercise lives on the ODE path.
• Adaptive λ schedule and local-tearing sweep (OpenModelica's adaptiveGlobal, equidistantLocal, adaptiveLocal modes) are out of scope; this PR ships the equidistantGlobal equivalent only.

Test plan

• Operator + both lowering passes (unit)
• Three nlsolve algorithms (unit, incl. sweep fallback + λ-reset-on-failure)
• Default inner solver via NonlinearSolve weak-dependency extension; SimpleNonlinearSolve fallback
• End-to-end integration through ODEProblem (default + both explicit modes + out-of-basin sweep rescue)
• Default-alg injection reaches NonlinearProblem / SteadyStateProblem / DAEProblem / SDEProblem
• OMC numerical parity on the Buildings PressureDrop fixture
• Graceful build/solve in an environment without NonlinearSolve
• docs/src/API/homotopy.md builds under Documenter (all @docs resolve)
• GROUP=Initialization regression: no new failures vs. master baseline
• No SciMLBase / NonlinearSolve / OrdinaryDiffEq changes (NonlinearSolve added only as [weakdeps])
• No @rule, no IfElse.ifelse, scalar @register_symbolic only (per Modelica's homotopy operator #1385 feedback)

(Diff: 21 files, +1267 / −7. New files: homotopy.jl, homotopy_sweep.jl, ext/MTKNonlinearSolveExt.jl, five test/homotopy_*.jl, docs/src/API/homotopy.md.)

Refs

• Modelica spec 3.7.4: https://specification.modelica.org/master/operators-and-expressions.html#homotopy
• OpenModelica homotopy flags (homotopyOnFirstTry, noHomotopyOnFirstTry, --homotopyApproach): OpenModelica User's Guide, Simulation Runtime Flags
• Modelica Buildings: Fluid/FixedResistances/PressureDrop.mo (from_dp=true branch)
• RFC thread + Aayush feedback turn: Modelica's homotopy operator #1385 (comment)
• Init use case: Homotopy Initialization #3988 — Related cyclic-guess pain: Initialization Failure when Providing Certain State Guesses (Cyclic Guesses Error) #4039

[BatchZero]

Hi @AayushSabharwal, the CI workflows seem to be waiting on a maintainer to approve the run. Could you kick it off when you get a chance?

[ChrisRackauckas]

This doesn't get the structure right. The homotopy solver should be separately designed with an MTK independent API in NonlinearSolve that is then targeted.

[BatchZero]

@ChrisRackauckas Understood — the continuation solver really doesn't belong in MTK. Apologies for the original structure: I went with the single-repo approach because it was the fastest way to get something landed, and I hope it didn't cause any frustration. I'll restructure it as an MTK-independent NonlinearSolve algorithm over a plain NonlinearProblem, and reduce MTK to lowering homotopy(actual, simplified) and targeting it. Below is the concrete plan; I'd like your read on it before I start, since it spans three repos.

Can the existing HomotopyContinuationJL / HomotopyNonlinearFunction be reused? I don't think so, and want to confirm the reasoning. That path traces solution paths through ℂⁿ and AD-differentiates the residual over the complex field (the wrapper reinterpret(Complex, ...)s the state before calling f, and solve filters only_real results). The systems this targets are non-holomorphic and real-only — homotopy(actual = sign(dp)*sqrt(abs(dp))*k, simplified = m_flow_nominal*dp/dp_nominal) and friends (abs, sign) can't be evaluated or differentiated on a complex tracker. The single-path GuessHomotopy ({false}) is structurally closer (one path from the initial guess) but still evaluates f in ℂ, so it hits the same wall; polynomialize/unpolynomialize doesn't help either, since it's for systems polynomial in a nonlinear basis with a closed-form inverse. So what I'd add is a sibling, not a reuse: a real-arithmetic convex embedding homotopy H(u, λ) = (1-λ)·simplified + λ·actual, λ swept 0→1 with warm-started inner Newton and step-size bisection on a non-converging step. This is the Modelica homotopy() semantics and a genuinely different algorithm class from polynomial homotopy continuation.

Planned changes, per repo:

1. SciMLBase — a small interface function type, EmbeddingHomotopyFunction <: AbstractNonlinearFunction, describing the homotopy family. Why here: it mirrors how HomotopyNonlinearFunction lives in SciMLBase so the algorithm and the problem-construction decouple across the package boundary; and I'd avoid touching HomotopyNonlinearFunction itself, whose polynomialize/unpolynomialize/denominator are polynomial-specific and shouldn't be generalized. (There's an open question on exactly what this type carries — see question 1 below.)

2. NonlinearSolve — the algorithm HomotopySweep <: AbstractNonlinearSolveAlgorithm plus CommonSolve.solve(prob::NonlinearProblem, alg::HomotopySweep). The λ-sweep / warm-start / λ schedule / step adaptivity live here; the inner solver is composable (defaults to the standard polyalgorithm, not hardcoded). Why here: this is the MTK-independent surface you asked for — a non-MTK user with a hand-built NonlinearProblem can solve(prob, HomotopySweep()). I'd ship a standalone test in NonlinearSolve's own suite (no MTK in deps) on a cold-Newton-divergent real system like sign(x)*sqrt(abs(x)) - c that the sweep rescues.

3. ModelingToolkit — lowering only: rewrite homotopy(actual, simplified) into an EmbeddingHomotopyFunction-backed NonlinearProblem, then target HomotopySweep (it can still plug into the existing OverrideInit.nlsolve slot for init). Add Modelica homotopy(actual, simplified) operator for initialization (spec 3.7.4) #4576 shrinks down to this. Why: symbolic lowering is MTK's job; the solver is not.

Merge order would be SciMLBase → NonlinearSolve → MTK, since each depends on the previous.

A few questions before I build it where you want it:

1. What the interface type carries — two reasonable encodings: (a) the two endpoints f (actual) + simplified, with the solver forming H(u, λ) = (1-λ)·simplified + λ·actual internally — keeps (simplified, actual) first-class and stays close to the Modelica operator; or (b) a single parametrized residual H(u, p, λ) plus a locator for which knob is λ — more general (any embedding, not just the convex one), with the convex combination built in the MTK lowering. (b) feels more in line with a generic continuation solver; (a) is more legible for this use case. Which do you prefer?
2. Home for the interface type — SciMLBase (mirroring HomotopyNonlinearFunction), or would you rather it start in NonlinearSolve to avoid a SciMLBase change? This is the difference between a 3-repo and a 2-repo chain.
3. Home for the algorithm — NonlinearSolve core, or a sublib like NonlinearSolveHomotopyContinuation? It has no heavy deps (just the inner solver), so core seems right to me, but your call.
4. Naming — I went with EmbeddingHomotopyFunction + HomotopySweep to keep the homotopy lineage and the Modelica connection while deliberately not colliding with the polynomial complex-tracking HomotopyContinuationJL. Any preference?

Does this structure look right to you, and is there anything you'd change before I start?

[ChrisRackauckas]

Can the existing HomotopyContinuationJL / HomotopyNonlinearFunction be reused? I don't think so, and want to confirm the reasoning.

Yeah that's pretty clear. Because it's only polynomial systems. It should probably be renamed PolynomialHomotopyProblem or something.

I'll restructure it as an MTK-independent NonlinearSolve algorithm over a plain NonlinearProblem, and reduce MTK to lowering homotopy(actual, simplified) and targeting it. Below is the concrete plan; I'd like your read on it before I start, since it spans three repos.

It probably needs a new problem type, since it needs more information than just a normal nonlinear problem. Or a function type, like you outline, that works fine too.

Home for the interface type — SciMLBase (mirroring HomotopyNonlinearFunction), or would you rather it start in NonlinearSolve to avoid a SciMLBase change? This is the difference between a 3-repo and a 2-repo chain.

Yes, the 3 repo chain is required.

[BatchZero]

Closing in favor of #4601: per the review here, the solver was redesigned as an MTK-independent API — SciMLBase.HomotopyProblem (SciML/SciMLBase.jl#1375, registered in SciMLBase v3.19.0) swept by NonlinearSolveBase.HomotopySweep (SciML/NonlinearSolve.jl#962, registered in NonlinearSolveBase v2.31.0), both now merged. #4601 carries the remaining MTK side: lowering homotopy(actual, simplified) and targeting that API.