From 36a6a95797f28d4ea6c782103a27c4a0b7df2e12 Mon Sep 17 00:00:00 2001 From: lf Date: Thu, 21 May 2026 14:14:05 +0800 Subject: [PATCH 01/28] test: failing Q1 basic test for homotopy rewrite_trivial --- lib/ModelingToolkitBase/test/homotopy_init.jl | 15 +++++++++++++++ 1 file changed, 15 insertions(+) create mode 100644 lib/ModelingToolkitBase/test/homotopy_init.jl diff --git a/lib/ModelingToolkitBase/test/homotopy_init.jl b/lib/ModelingToolkitBase/test/homotopy_init.jl new file mode 100644 index 0000000000..76a6f74bce --- /dev/null +++ b/lib/ModelingToolkitBase/test/homotopy_init.jl @@ -0,0 +1,15 @@ +using Test +using ModelingToolkitBase +using ModelingToolkitBase: rewrite_trivial +using Symbolics + +@testset "homotopy operator — L0 trivial rewrite" begin + @testset "Q1 basic" begin + @variables x + @parameters p + expr = homotopy(x^2 - p, x - sqrt(p)) + rewritten = rewrite_trivial(expr) + @test isequal(Symbolics.unwrap(rewritten), Symbolics.unwrap(x^2 - p)) + @test !occursin("homotopy", repr(rewritten)) + end +end From 8c28ed9f265b62e921c9a793de0e803cfbede587 Mon Sep 17 00:00:00 2001 From: lf Date: Thu, 21 May 2026 14:19:47 +0800 Subject: [PATCH 02/28] feat: register homotopy operator + rewrite_trivial (L0 trivial form) Implements Modelica spec 3.7.4 trivial form: homotopy(a, s) is a registered symbolic primitive that rewrites to actual=a. Hand-written TermInterface walk (not @rule) per #1385 feedback. Metadata + symtype propagated on rebuilt nodes. Runtime fallback methods for Real and AbstractArray inputs. --- .../src/ModelingToolkitBase.jl | 2 + .../src/systems/homotopy.jl | 71 +++++++++++++++++++ 2 files changed, 73 insertions(+) create mode 100644 lib/ModelingToolkitBase/src/systems/homotopy.jl diff --git a/lib/ModelingToolkitBase/src/ModelingToolkitBase.jl b/lib/ModelingToolkitBase/src/ModelingToolkitBase.jl index da845565bf..3a71cb2124 100644 --- a/lib/ModelingToolkitBase/src/ModelingToolkitBase.jl +++ b/lib/ModelingToolkitBase/src/ModelingToolkitBase.jl @@ -194,6 +194,7 @@ include("systems/ir_info.jl") include("systems/codegen_utils.jl") include("problems/docs.jl") include("systems/codegen.jl") +include("systems/homotopy.jl") include("systems/problem_utils.jl") include("problems/compatibility.jl") @@ -270,6 +271,7 @@ export liouville_transform, change_independent_variable, export respecialize export PDESystem export Differential, expand_derivatives, @derivatives +export homotopy export Equation export Term export SymScope, LocalScope, ParentScope, GlobalScope diff --git a/lib/ModelingToolkitBase/src/systems/homotopy.jl b/lib/ModelingToolkitBase/src/systems/homotopy.jl new file mode 100644 index 0000000000..57e528fee1 --- /dev/null +++ b/lib/ModelingToolkitBase/src/systems/homotopy.jl @@ -0,0 +1,71 @@ +# Modelica homotopy operator — L0 trivial form (spec 3.7.4). +# `homotopy(actual, simplified)` is a symbolic primitive; at runtime and in +# the init pipeline it evaluates to `actual`. The `simplified` argument is +# carried symbolically so future PRs (L1 parameter sweep) can use it for +# continuation without changing call sites. + +using Symbolics +using SymbolicUtils + +@register_symbolic homotopy(actual, simplified) + +""" + homotopy(actual, simplified) + +Modelica homotopy operator (spec 3.7.4). At runtime and at L0 trivial rewrite +time, equivalent to `actual`. Future PRs may use `simplified` as a starting +point for parameter-sweep continuation during initialization. +""" +homotopy(actual::Real, simplified::Real) = actual +homotopy(actual::AbstractArray, simplified::AbstractArray) = actual + +""" + rewrite_trivial(expr) + +Recursively replace every `homotopy(a, s)` node in `expr` with `a`. Preserves +`metadata` and `symtype` on every rebuilt node. Implemented as a hand-written +walk over `TermInterface` — does NOT use `@rule` (per AayushSabharwal feedback +in #1385: per-node matcher overhead is unacceptable at scale). +""" +function rewrite_trivial(expr) + x = Symbolics.unwrap(expr) + return _rewrite_trivial(x) +end + +function _rewrite_trivial(x) + if !SymbolicUtils.iscall(x) + return x + end + op = SymbolicUtils.operation(x) + args = SymbolicUtils.arguments(x) + new_args = map(_rewrite_trivial, args) + if op === homotopy + length(new_args) == 2 || throw(ArgumentError( + "homotopy() expects exactly 2 arguments, got $(length(new_args))")) + return new_args[1] + end + return SymbolicUtils.maketerm( + typeof(x), op, new_args, + SymbolicUtils.metadata(x), + ) +end + +""" + has_homotopy(expr) + +Return `true` iff `expr` contains at least one `homotopy(...)` node. +""" +function has_homotopy(expr) + x = Symbolics.unwrap(expr) + return _has_homotopy(x) +end + +function _has_homotopy(x) + if !SymbolicUtils.iscall(x) + return false + end + if SymbolicUtils.operation(x) === homotopy + return true + end + return any(_has_homotopy, SymbolicUtils.arguments(x)) +end From 2b7b770ad609867abe31ca1fee7358c1c1369779 Mon Sep 17 00:00:00 2001 From: lf Date: Thu, 21 May 2026 14:32:01 +0800 Subject: [PATCH 03/28] refactor: align homotopy.jl with MTKBase style + tighten test MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Per code review: - Drop top-level 'using Symbolics' / 'using SymbolicUtils' inside the included file — MTKBase parent module already imports iscall, operation, arguments, maketerm, metadata bare (L21-23 of ModelingToolkitBase.jl) and re-exports Symbolics, so qualified prefixes are redundant noise versus surrounding files (utils.jl, codegen_utils.jl, if_lifting.jl, etc.). - Drop unreachable AbstractArray runtime fallback (@register_symbolic is scalar-only; users go through broadcast). - Drop arity ArgumentError guard in _rewrite_trivial (@register_symbolic enforces 2-arity at parse time; guard is dead code). - Soften docstring: metadata is preserved explicitly via maketerm arg; symtype is inferred by maketerm itself, not preserved by us. - Replace brittle !occursin("homotopy", repr(rewritten)) assertion with !has_homotopy(rewritten), which uses the public predicate and gives has_homotopy test coverage as a bonus. --- .../src/systems/homotopy.jl | 40 ++++++------------- lib/ModelingToolkitBase/test/homotopy_init.jl | 4 +- 2 files changed, 15 insertions(+), 29 deletions(-) diff --git a/lib/ModelingToolkitBase/src/systems/homotopy.jl b/lib/ModelingToolkitBase/src/systems/homotopy.jl index 57e528fee1..e01b720b6a 100644 --- a/lib/ModelingToolkitBase/src/systems/homotopy.jl +++ b/lib/ModelingToolkitBase/src/systems/homotopy.jl @@ -1,12 +1,3 @@ -# Modelica homotopy operator — L0 trivial form (spec 3.7.4). -# `homotopy(actual, simplified)` is a symbolic primitive; at runtime and in -# the init pipeline it evaluates to `actual`. The `simplified` argument is -# carried symbolically so future PRs (L1 parameter sweep) can use it for -# continuation without changing call sites. - -using Symbolics -using SymbolicUtils - @register_symbolic homotopy(actual, simplified) """ @@ -17,37 +8,32 @@ time, equivalent to `actual`. Future PRs may use `simplified` as a starting point for parameter-sweep continuation during initialization. """ homotopy(actual::Real, simplified::Real) = actual -homotopy(actual::AbstractArray, simplified::AbstractArray) = actual """ rewrite_trivial(expr) Recursively replace every `homotopy(a, s)` node in `expr` with `a`. Preserves -`metadata` and `symtype` on every rebuilt node. Implemented as a hand-written -walk over `TermInterface` — does NOT use `@rule` (per AayushSabharwal feedback -in #1385: per-node matcher overhead is unacceptable at scale). +`metadata` on every rebuilt node; `symtype` is inferred by `maketerm`. +Implemented as a hand-written walk over the TermInterface protocol — does NOT +use `@rule` (per AayushSabharwal feedback in #1385: per-node matcher overhead +is unacceptable at scale). """ function rewrite_trivial(expr) - x = Symbolics.unwrap(expr) + x = unwrap(expr) return _rewrite_trivial(x) end function _rewrite_trivial(x) - if !SymbolicUtils.iscall(x) + if !iscall(x) return x end - op = SymbolicUtils.operation(x) - args = SymbolicUtils.arguments(x) + op = operation(x) + args = arguments(x) new_args = map(_rewrite_trivial, args) if op === homotopy - length(new_args) == 2 || throw(ArgumentError( - "homotopy() expects exactly 2 arguments, got $(length(new_args))")) return new_args[1] end - return SymbolicUtils.maketerm( - typeof(x), op, new_args, - SymbolicUtils.metadata(x), - ) + return maketerm(typeof(x), op, new_args, metadata(x)) end """ @@ -56,16 +42,16 @@ end Return `true` iff `expr` contains at least one `homotopy(...)` node. """ function has_homotopy(expr) - x = Symbolics.unwrap(expr) + x = unwrap(expr) return _has_homotopy(x) end function _has_homotopy(x) - if !SymbolicUtils.iscall(x) + if !iscall(x) return false end - if SymbolicUtils.operation(x) === homotopy + if operation(x) === homotopy return true end - return any(_has_homotopy, SymbolicUtils.arguments(x)) + return any(_has_homotopy, arguments(x)) end diff --git a/lib/ModelingToolkitBase/test/homotopy_init.jl b/lib/ModelingToolkitBase/test/homotopy_init.jl index 76a6f74bce..4c9a59a3a6 100644 --- a/lib/ModelingToolkitBase/test/homotopy_init.jl +++ b/lib/ModelingToolkitBase/test/homotopy_init.jl @@ -1,6 +1,6 @@ using Test using ModelingToolkitBase -using ModelingToolkitBase: rewrite_trivial +using ModelingToolkitBase: rewrite_trivial, has_homotopy using Symbolics @testset "homotopy operator — L0 trivial rewrite" begin @@ -10,6 +10,6 @@ using Symbolics expr = homotopy(x^2 - p, x - sqrt(p)) rewritten = rewrite_trivial(expr) @test isequal(Symbolics.unwrap(rewritten), Symbolics.unwrap(x^2 - p)) - @test !occursin("homotopy", repr(rewritten)) + @test !has_homotopy(rewritten) end end From 86b613ed9ea98d6273ca4bc7cac01694e5d5bef5 Mon Sep 17 00:00:00 2001 From: lf Date: Thu, 21 May 2026 14:35:18 +0800 Subject: [PATCH 04/28] test: Q2 nested homotopy collapses to innermost actual --- lib/ModelingToolkitBase/test/homotopy_init.jl | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/lib/ModelingToolkitBase/test/homotopy_init.jl b/lib/ModelingToolkitBase/test/homotopy_init.jl index 4c9a59a3a6..85d7d65144 100644 --- a/lib/ModelingToolkitBase/test/homotopy_init.jl +++ b/lib/ModelingToolkitBase/test/homotopy_init.jl @@ -12,4 +12,20 @@ using Symbolics @test isequal(Symbolics.unwrap(rewritten), Symbolics.unwrap(x^2 - p)) @test !has_homotopy(rewritten) end + + @testset "Q2 nested" begin + @variables x + @parameters p + inner = homotopy(x^2 - p, x - sqrt(p)) + outer = homotopy(inner, x - 1) + @test isequal(Symbolics.unwrap(rewrite_trivial(outer)), + Symbolics.unwrap(x^2 - p)) + + triple = homotopy(homotopy(homotopy(x, x + 1), x + 2), x + 3) + @test isequal(Symbolics.unwrap(rewrite_trivial(triple)), + Symbolics.unwrap(x)) + + @test !has_homotopy(rewrite_trivial(outer)) + @test !has_homotopy(rewrite_trivial(triple)) + end end From 414261c0701daa4d2b90efda3b0a0e14a176bbd9 Mon Sep 17 00:00:00 2001 From: lf Date: Thu, 21 May 2026 14:37:27 +0800 Subject: [PATCH 05/28] test: Q3 homotopy preserved inside Base.ifelse branch structure --- lib/ModelingToolkitBase/test/homotopy_init.jl | 23 +++++++++++++++++++ 1 file changed, 23 insertions(+) diff --git a/lib/ModelingToolkitBase/test/homotopy_init.jl b/lib/ModelingToolkitBase/test/homotopy_init.jl index 85d7d65144..6977dbf08b 100644 --- a/lib/ModelingToolkitBase/test/homotopy_init.jl +++ b/lib/ModelingToolkitBase/test/homotopy_init.jl @@ -28,4 +28,27 @@ using Symbolics @test !has_homotopy(rewrite_trivial(outer)) @test !has_homotopy(rewrite_trivial(triple)) end + + @testset "Q3 Base.ifelse branches" begin + @variables x + @parameters p + cond = p > 0 + branch_expr = Base.ifelse( + cond, + homotopy(x^2 - p, x - sqrt(abs(p))), + homotopy(-(x^2) - p, x + sqrt(abs(p))), + ) + rewritten = rewrite_trivial(branch_expr) + + # Outer ifelse structure preserved + @test occursin("ifelse", repr(rewritten)) + # Both branches' homotopy nodes are gone + @test !has_homotopy(rewritten) + + # Folding the cond to true/false picks the right actual + true_folded = Symbolics.simplify(Symbolics.substitute(rewritten, Dict(cond => true))) + false_folded = Symbolics.simplify(Symbolics.substitute(rewritten, Dict(cond => false))) + @test isequal(Symbolics.unwrap(true_folded), Symbolics.unwrap(x^2 - p)) + @test isequal(Symbolics.unwrap(false_folded), Symbolics.unwrap(-(x^2) - p)) + end end From 20be7b0ba4b4cee72ab7c011c318c284278fa848 Mon Sep 17 00:00:00 2001 From: lf Date: Thu, 21 May 2026 14:38:00 +0800 Subject: [PATCH 06/28] test: Q4 vectorized homotopy via broadcast --- lib/ModelingToolkitBase/test/homotopy_init.jl | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/lib/ModelingToolkitBase/test/homotopy_init.jl b/lib/ModelingToolkitBase/test/homotopy_init.jl index 6977dbf08b..92be761255 100644 --- a/lib/ModelingToolkitBase/test/homotopy_init.jl +++ b/lib/ModelingToolkitBase/test/homotopy_init.jl @@ -51,4 +51,18 @@ using Symbolics @test isequal(Symbolics.unwrap(true_folded), Symbolics.unwrap(x^2 - p)) @test isequal(Symbolics.unwrap(false_folded), Symbolics.unwrap(-(x^2) - p)) end + + @testset "Q4 broadcast" begin + @variables x[1:3] p[1:3] + actuals = [x[i]^2 - p[i] for i in 1:3] + simplifieds = [x[i] - p[i] for i in 1:3] + vec_expr = homotopy.(actuals, simplifieds) + rewritten = rewrite_trivial.(vec_expr) + @test length(rewritten) == 3 + for i in 1:3 + @test isequal(Symbolics.unwrap(rewritten[i]), + Symbolics.unwrap(actuals[i])) + @test !has_homotopy(rewritten[i]) + end + end end From 4bcc771fe81dd77ad89b9d9b5772fee66a685a0a Mon Sep 17 00:00:00 2001 From: lf Date: Thu, 21 May 2026 14:50:29 +0800 Subject: [PATCH 07/28] test: failing integration test for homotopy in InitializationProblem Adds has_homotopy_in_equations and rewrite_trivial_in_equations helpers to homotopy.jl (internal use; not exported). Adds RED integration test asserting that an InitializationProblem built from a NonlinearSystem containing homotopy(x^2 - p, x - 1) = 0 solves the actual branch (x^2 = p) rather than the simplified branch (x = 1). Currently fails because the init pipeline does not yet apply rewrite_trivial. --- .../src/systems/homotopy.jl | 27 +++++++++++++++++++ lib/ModelingToolkitBase/test/homotopy_init.jl | 24 +++++++++++++++++ 2 files changed, 51 insertions(+) diff --git a/lib/ModelingToolkitBase/src/systems/homotopy.jl b/lib/ModelingToolkitBase/src/systems/homotopy.jl index e01b720b6a..5c0f0abfd0 100644 --- a/lib/ModelingToolkitBase/src/systems/homotopy.jl +++ b/lib/ModelingToolkitBase/src/systems/homotopy.jl @@ -55,3 +55,30 @@ function _has_homotopy(x) end return any(_has_homotopy, arguments(x)) end + +""" + has_homotopy_in_equations(eqs) + +Return `true` iff any equation in `eqs` (lhs or rhs) contains a `homotopy(...)` +node. `eqs` is an iterable of `Equation`. +""" +function has_homotopy_in_equations(eqs) + for eq in eqs + if has_homotopy(eq.lhs) || has_homotopy(eq.rhs) + return true + end + end + return false +end + +""" + rewrite_trivial_in_equations(eqs) + +Return a new vector of `Equation`s with every `homotopy(a, s)` replaced by `a` +on both lhs and rhs. Original `eqs` not mutated; the system caller is +responsible for swapping the equation vector into the system. +""" +function rewrite_trivial_in_equations(eqs) + return [Equation(rewrite_trivial(eq.lhs), rewrite_trivial(eq.rhs)) + for eq in eqs] +end diff --git a/lib/ModelingToolkitBase/test/homotopy_init.jl b/lib/ModelingToolkitBase/test/homotopy_init.jl index 92be761255..5f79052e5d 100644 --- a/lib/ModelingToolkitBase/test/homotopy_init.jl +++ b/lib/ModelingToolkitBase/test/homotopy_init.jl @@ -65,4 +65,28 @@ using Symbolics @test !has_homotopy(rewritten[i]) end end + + @testset "Integration: homotopy in NonlinearSystem init" begin + using ModelingToolkitBase: System, mtkcompile + using ModelingToolkitBase: InitializationProblem + using NonlinearSolve: NewtonRaphson, solve + using SciMLBase + + @variables x + @parameters p + # Equation: homotopy(actual = x^2 - p, simplified = x - 1) = 0 + # L0 trivial must solve actual: x^2 = p ⇒ x = ±√p. + # If rewrite is broken and simplified got selected: x = 1. + # At p = 9: actual root x = 3 (or -3); simplified root x = 1. + # The x^2 ≈ p assertion distinguishes the two outcomes. + eqs = [0 ~ homotopy(x^2 - p, x - 1)] + @named sys = System(eqs, [x], [p]) + sys = mtkcompile(sys) + + prob = InitializationProblem{false}(sys, nothing, Dict(p => 9.0, x => 2.5)) + sol = solve(prob, NewtonRaphson(); abstol = 1e-10, reltol = 1e-10) + @test SciMLBase.successful_retcode(sol) + @test abs(sol.u[1]^2 - 9.0) < 1e-6 # actual equation solved (NOT simplified, which would give x=1) + @test abs(sol.u[1] - 1.0) > 0.5 # not the simplified root + end end From a47c6fefc3f941084255b59cfa2854ae1a2cd800 Mon Sep 17 00:00:00 2001 From: lf Date: Thu, 21 May 2026 15:06:15 +0800 Subject: [PATCH 08/28] test: refine RED integration test to ODESystem with algebraic homotopy The previous NonlinearSystem variant never exercised the rewrite hook because generate_initializesystem_timeindependent does not migrate the parent system's algebraic equations into the initialization system, so the failure mode was an empty u0 unrelated to homotopy. The new design uses an ODESystem with an algebraic constraint `0 ~ homotopy(y^2 - p, y - 1)` and a free init guess for `y`. This routes through generate_initializesystem_timevarying, which preserves the homotopy-bearing equation in the init system. The structural assertion `!has_homotopy_in_equations(equations(prob.f.sys))` distinguishes hooked from unhooked: without the rewrite, the init system carries the opaque homotopy() symbolic call even though runtime dispatch happens to give a numerically correct answer via the homotopy(::Real, ::Real) fallback. --- lib/ModelingToolkitBase/test/homotopy_init.jl | 33 ++++++++++++------- 1 file changed, 21 insertions(+), 12 deletions(-) diff --git a/lib/ModelingToolkitBase/test/homotopy_init.jl b/lib/ModelingToolkitBase/test/homotopy_init.jl index 5f79052e5d..f2f28dbd5a 100644 --- a/lib/ModelingToolkitBase/test/homotopy_init.jl +++ b/lib/ModelingToolkitBase/test/homotopy_init.jl @@ -66,27 +66,36 @@ using Symbolics end end - @testset "Integration: homotopy in NonlinearSystem init" begin + @testset "Integration: homotopy in ODESystem init" begin using ModelingToolkitBase: System, mtkcompile using ModelingToolkitBase: InitializationProblem + using ModelingToolkitBase: t_nounits as t, D_nounits as D + using ModelingToolkitBase: has_homotopy_in_equations, equations using NonlinearSolve: NewtonRaphson, solve using SciMLBase - @variables x + @variables x(t) y(t) @parameters p - # Equation: homotopy(actual = x^2 - p, simplified = x - 1) = 0 - # L0 trivial must solve actual: x^2 = p ⇒ x = ±√p. - # If rewrite is broken and simplified got selected: x = 1. - # At p = 9: actual root x = 3 (or -3); simplified root x = 1. - # The x^2 ≈ p assertion distinguishes the two outcomes. - eqs = [0 ~ homotopy(x^2 - p, x - 1)] - @named sys = System(eqs, [x], [p]) + # `y` is algebraic, constrained by homotopy(actual = y^2 - p, simplified = y - 1). + # L0 trivial rewrite must replace the constraint with `y^2 - p = 0`, so init solves + # to y = ±√p. At p = 9 with guess y = 2.5, init converges to y ≈ 3. If the rewrite + # never fired, the init system would carry the opaque `homotopy(...)` symbolic call + # — which is the structural assertion below. + eqs = [D(x) ~ -x, + 0 ~ homotopy(y^2 - p, y - 1)] + @named sys = System(eqs, t; guesses = [y => 2.5]) sys = mtkcompile(sys) - prob = InitializationProblem{false}(sys, nothing, Dict(p => 9.0, x => 2.5)) + prob = InitializationProblem{false}(sys, 0.0, Dict(x => 1.0, p => 9.0)) + + # Structural: after the init pipeline applies rewrite_trivial, no equation in + # the wrapped initialization system should still contain a `homotopy(...)` node. + @test !has_homotopy_in_equations(equations(prob.f.sys)) + + # Numerical: init must converge to the actual root. sol = solve(prob, NewtonRaphson(); abstol = 1e-10, reltol = 1e-10) @test SciMLBase.successful_retcode(sol) - @test abs(sol.u[1]^2 - 9.0) < 1e-6 # actual equation solved (NOT simplified, which would give x=1) - @test abs(sol.u[1] - 1.0) > 0.5 # not the simplified root + @test abs(sol.u[1]^2 - 9.0) < 1e-6 # actual equation y^2 = p satisfied + @test abs(sol.u[1] - 1.0) > 0.5 # not the simplified root (y = 1) end end From 3cca63680ee93126119d3b3831e91508482152b7 Mon Sep 17 00:00:00 2001 From: lf Date: Thu, 21 May 2026 15:06:31 +0800 Subject: [PATCH 09/28] feat: rewrite homotopy to actual inside InitializationProblem Detects homotopy() in initialization-system equations and applies rewrite_trivial. Gated by has_homotopy_in_equations so non-homotopy systems pay zero traversal cost. Solves the L0 trivial requirement of Modelica spec 3.7.4. The hook lands before mtkcompile(isys) so that downstream structural transforms operate on the actual equations rather than opaque homotopy() calls; future PRs that lower into a parameter sweep will replace this branch when a HomotopySweep algorithm is in use. --- .../src/problems/initializationproblem.jl | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/lib/ModelingToolkitBase/src/problems/initializationproblem.jl b/lib/ModelingToolkitBase/src/problems/initializationproblem.jl index e24346cfea..7056591350 100644 --- a/lib/ModelingToolkitBase/src/problems/initializationproblem.jl +++ b/lib/ModelingToolkitBase/src/problems/initializationproblem.jl @@ -82,6 +82,17 @@ All other keyword arguments are forwarded to the wrapped problem constructor. binds[get_iv(sys)::SymbolicT] = Symbolics.COMMON_ZERO @set! isys.bindings = ROSymmapT(binds) end + # L0 trivial homotopy rewrite (Modelica spec 3.7.4 trivial form). + # Only runs when `homotopy(...)` is actually present, so non-homotopy + # systems pay zero overhead. PR2 (L1 parameter sweep) will replace this + # branch with λ-lowering when an explicit HomotopySweep algorithm is in use. + # Applied before `mtkcompile` so that downstream structural transforms + # operate on the actual equations, not on opaque `homotopy(...)` calls. + if has_homotopy_in_equations(equations(isys)) + new_eqs = rewrite_trivial_in_equations(equations(isys)) + @set! isys.eqs = new_eqs + end + if simplify_system isys = mtkcompile(isys; fully_determined, split = is_split(sys), initsys_mtkcompile_kwargs...) end From 743fd9e254ad7cc613e04b159c5668da47eabc5a Mon Sep 17 00:00:00 2001 From: lf Date: Thu, 21 May 2026 15:06:47 +0800 Subject: [PATCH 10/28] test: register homotopy_init.jl in Initialization group The new homotopy_init.jl covers rewrite_trivial (unit) plus an ODESystem InitializationProblem integration check. It belongs alongside the other Initialization-group entries (guess_propagation, initializationsystem, initial_values) rather than the broader InterfaceI bucket. --- lib/ModelingToolkitBase/test/runtests.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/lib/ModelingToolkitBase/test/runtests.jl b/lib/ModelingToolkitBase/test/runtests.jl index e024cc0200..732d61dfe5 100644 --- a/lib/ModelingToolkitBase/test/runtests.jl +++ b/lib/ModelingToolkitBase/test/runtests.jl @@ -63,6 +63,7 @@ end if GROUP == "All" || GROUP == "Initialization" @safetestset "Guess Propagation" include("guess_propagation.jl") @safetestset "InitializationSystem Test" include("initializationsystem.jl") + @safetestset "Homotopy Operator L0 Trivial" include("homotopy_init.jl") @safetestset "Initial Values Test" include("initial_values.jl") end From e2b82a34730d28972508404ce876d2e191924015 Mon Sep 17 00:00:00 2001 From: lf Date: Thu, 21 May 2026 15:45:27 +0800 Subject: [PATCH 11/28] test: Buildings PressureDrop-style fixture exercises homotopy init Mirrors the structure of Buildings/Fluid/FixedResistances/PressureDrop.mo (from_dp=true branch): nonlinear actual (sqrt-based turbulent law) + linear simplified (nominal-point scaling). PR1 asserts (a) hook fired (no residual homotopy in init equations), (b) init converged, (c) m_flow numerically matches the actual root sqrt(5), not the simplified root m_flow_nominal*dp/dp_nominal = 1.0. Numerical OMC parity deferred to PR2. --- lib/ModelingToolkitBase/test/homotopy_init.jl | 58 +++++++++++++++++++ 1 file changed, 58 insertions(+) diff --git a/lib/ModelingToolkitBase/test/homotopy_init.jl b/lib/ModelingToolkitBase/test/homotopy_init.jl index f2f28dbd5a..3f1bb048cb 100644 --- a/lib/ModelingToolkitBase/test/homotopy_init.jl +++ b/lib/ModelingToolkitBase/test/homotopy_init.jl @@ -98,4 +98,62 @@ using Symbolics @test abs(sol.u[1]^2 - 9.0) < 1e-6 # actual equation y^2 = p satisfied @test abs(sol.u[1] - 1.0) > 0.5 # not the simplified root (y = 1) end + + @testset "Buildings PressureDrop fixture (PR1: init doesn't break)" begin + # Modelica Buildings Fluid/FixedResistances/PressureDrop.mo from_dp=true branch: + # m_flow = homotopy( + # actual = basicFlowFunction_dp(dp, k, m_flow_turbulent), + # simplified = m_flow_nominal_pos * dp / dp_nominal_pos) + # PR1 scope = "init runs and converges to the *actual* equation, not + # the simplified one". Numerical OMC parity is deferred to PR2. + # Uses the ODESystem-with-algebraic-constraint pattern (same as the + # Integration testset) so the hook in InitializationProblem fires + # symbolically; NonlinearSystem direct construction does not migrate + # parent algebraic eqs into init system. + using ModelingToolkitBase: System, mtkcompile + using ModelingToolkitBase: t_nounits as t, D_nounits as D + using ModelingToolkitBase: has_homotopy_in_equations, equations + using OrdinaryDiffEqRosenbrock: Rodas5P + using OrdinaryDiffEqNonlinearSolve # needed for DAE init nlsolve + using SciMLBase + + @variables x(t) m_flow(t) + @parameters dp k m_flow_nominal dp_nominal + + # Turbulent sqrt-law actual, linear-nominal simplified. + basic_flow(dp_val, k_val) = sign(dp_val) * sqrt(abs(dp_val)) * k_val + + # `dp` is a parameter (boundary fixed externally); `m_flow` is + # algebraic, pinned by the homotopy constraint. `x` is a dummy + # differential state so the system is a well-posed DAE (mirrors + # the Integration testset's structure: 1 ODE + 1 algebraic eq). + eqs = [ + D(x) ~ -x, + 0 ~ m_flow - homotopy( + basic_flow(dp, k), + m_flow_nominal * dp / dp_nominal, + ), + ] + @named sys = System(eqs, t; guesses = [m_flow => 2.0]) + sys = mtkcompile(sys) + + prob = ODEProblem(sys, + Dict(x => 1.0, dp => 5.0, k => 1.0, + m_flow_nominal => 1.0, dp_nominal => 5.0), + (0.0, 1.0)) + sol = solve(prob, Rodas5P(); abstol = 1e-10, reltol = 1e-10) + + # Hook fired: no residual homotopy in init equations + @test !has_homotopy_in_equations(equations(prob.f.initialization_data.initializeprob.f.sys)) + # Solver succeeded + @test SciMLBase.successful_retcode(sol) + # m_flow at t=0 should solve the *actual* equation: + # m_flow = sign(dp) * sqrt(abs(dp)) * k at dp=5, k=1 + # m_flow = sqrt(5) ≈ 2.236 + # If the rewrite had failed and simplified were active: + # m_flow = m_flow_nominal * dp / dp_nominal = 1.0 * 5/5 = 1.0 + # The sqrt(5) result distinguishes the two outcomes. + m_flow_at_t0 = sol[m_flow][1] + @test abs(m_flow_at_t0 - sqrt(5.0)) < 1e-6 # actual root + end end From fed7d600ddf5926872bbf2f4c749f9179764c691 Mon Sep 17 00:00:00 2001 From: lf Date: Mon, 25 May 2026 00:55:26 +0800 Subject: [PATCH 12/28] feat: rewrite_with_lambda lowering for L1 sweep path --- .../src/systems/homotopy.jl | 50 ++++++++++++++++ .../test/homotopy_lowering.jl | 60 +++++++++++++++++++ lib/ModelingToolkitBase/test/runtests.jl | 1 + 3 files changed, 111 insertions(+) create mode 100644 lib/ModelingToolkitBase/test/homotopy_lowering.jl diff --git a/lib/ModelingToolkitBase/src/systems/homotopy.jl b/lib/ModelingToolkitBase/src/systems/homotopy.jl index 5c0f0abfd0..3ce95d5f2d 100644 --- a/lib/ModelingToolkitBase/src/systems/homotopy.jl +++ b/lib/ModelingToolkitBase/src/systems/homotopy.jl @@ -82,3 +82,53 @@ function rewrite_trivial_in_equations(eqs) return [Equation(rewrite_trivial(eq.lhs), rewrite_trivial(eq.rhs)) for eq in eqs] end + +""" + rewrite_with_lambda(expr, λ = nothing) + +Recursively replace every `homotopy(a, s)` node in `expr` with `(1 - λ)*s + λ*a`, +where `λ` is a single shared parameter for the whole expression (allocated lazily +if not supplied, default name `__homotopy_λ`, default value `1.0`). +Returns `(new_expr, λ)`. + +At λ=1 the lowered expression reduces numerically to `actual` (trivial form); +at λ=0 it reduces to `simplified`. `HomotopySweep` walks λ from 0 → 1. + +Hand-written TermInterface walk; no `@rule` (per AayushSabharwal #1385 feedback). +""" +function rewrite_with_lambda(expr, λ = nothing) + if λ === nothing + λ = only(@parameters __homotopy_λ = 1.0) + end + x = unwrap(expr) + return _rewrite_with_lambda(x, λ), λ +end + +function _rewrite_with_lambda(x, λ) + if !iscall(x) + return x + end + op = operation(x) + args = arguments(x) + new_args = map(arg -> _rewrite_with_lambda(arg, λ), args) + if op === homotopy + a, s = new_args[1], new_args[2] + return (1 - λ) * s + λ * a + end + return maketerm(typeof(x), op, new_args, metadata(x)) +end + +""" + rewrite_with_lambda_in_equations(eqs, λ = nothing) + +Return `(new_eqs, λ)`. All `homotopy(a, s)` nodes across all equations are +replaced using a single shared `λ` parameter. Original `eqs` not mutated. +""" +function rewrite_with_lambda_in_equations(eqs, λ = nothing) + if λ === nothing + λ = only(@parameters __homotopy_λ = 1.0) + end + new_eqs = [Equation(_rewrite_with_lambda(unwrap(eq.lhs), λ), + _rewrite_with_lambda(unwrap(eq.rhs), λ)) for eq in eqs] + return new_eqs, λ +end diff --git a/lib/ModelingToolkitBase/test/homotopy_lowering.jl b/lib/ModelingToolkitBase/test/homotopy_lowering.jl new file mode 100644 index 0000000000..ce1944506e --- /dev/null +++ b/lib/ModelingToolkitBase/test/homotopy_lowering.jl @@ -0,0 +1,60 @@ +using Test +using ModelingToolkitBase +using ModelingToolkitBase: rewrite_with_lambda, rewrite_with_lambda_in_equations, + has_homotopy, has_homotopy_in_equations +using Symbolics + +@testset "homotopy operator — L1 lowering rewrite_with_lambda" begin + @testset "L1-Q1 basic single node" begin + @variables x + @parameters p + expr = homotopy(x^2 - p, x - sqrt(p)) + rewritten, λ = rewrite_with_lambda(expr) + @test !has_homotopy(rewritten) + at_one = Symbolics.simplify(Symbolics.substitute(rewritten, Dict(λ => 1.0))) + at_zero = Symbolics.simplify(Symbolics.substitute(rewritten, Dict(λ => 0.0))) + # At λ=1 the lowered expression reduces to `actual`; at λ=0 to `simplified`. + @test isequal(Symbolics.unwrap(at_one), Symbolics.unwrap(x^2 - p)) + @test isequal(Symbolics.unwrap(at_zero), Symbolics.unwrap(x - sqrt(p))) + end + + @testset "L1-Q2 nested homotopy collapses with single λ" begin + @variables x + @parameters p + inner = homotopy(x^2 - p, x - sqrt(p)) + outer = homotopy(inner, x - 1) + rewritten, λ = rewrite_with_lambda(outer) + @test !has_homotopy(rewritten) + at_one = Symbolics.simplify(Symbolics.substitute(rewritten, Dict(λ => 1.0))) + # Nested homotopy at λ=1 collapses fully to the innermost `actual`. + @test isequal(Symbolics.unwrap(at_one), Symbolics.unwrap(x^2 - p)) + end + + @testset "L1-Q3 multiple nodes share the same λ" begin + @variables x y + @parameters p q + e1 = homotopy(x^2 - p, x - 1) + e2 = homotopy(y^2 - q, y - 1) + eqs = [0 ~ e1, 0 ~ e2] + new_eqs, λ = rewrite_with_lambda_in_equations(eqs) + @test length(new_eqs) == 2 + params_in_rhs = Set{Symbol}() + for eq in new_eqs + for var in Symbolics.get_variables(eq.rhs) + push!(params_in_rhs, nameof(Symbolics.unwrap(var))) + end + end + @test :__homotopy_λ in params_in_rhs + @test !has_homotopy_in_equations(new_eqs) + end + + @testset "L1-Q4 no-homotopy equations pass through unchanged" begin + @variables x + @parameters p + eqs = [0 ~ x^2 - p, 0 ~ x + 1] + new_eqs, λ = rewrite_with_lambda_in_equations(eqs) + @test isequal(new_eqs[1].rhs, eqs[1].rhs) + @test isequal(new_eqs[2].rhs, eqs[2].rhs) + @test λ !== nothing + end +end diff --git a/lib/ModelingToolkitBase/test/runtests.jl b/lib/ModelingToolkitBase/test/runtests.jl index 732d61dfe5..9eb123cafb 100644 --- a/lib/ModelingToolkitBase/test/runtests.jl +++ b/lib/ModelingToolkitBase/test/runtests.jl @@ -64,6 +64,7 @@ end @safetestset "Guess Propagation" include("guess_propagation.jl") @safetestset "InitializationSystem Test" include("initializationsystem.jl") @safetestset "Homotopy Operator L0 Trivial" include("homotopy_init.jl") + @safetestset "Homotopy lowering" begin include("homotopy_lowering.jl") end @safetestset "Initial Values Test" include("initial_values.jl") end From 2ef5255fdd24a57b5a68df55deb0b9dd5052a917 Mon Sep 17 00:00:00 2001 From: lf Date: Mon, 25 May 2026 01:53:35 +0800 Subject: [PATCH 13/28] feat: HomotopySweep, TrivialHomotopy, TrivialThenSweep algorithms --- .../src/ModelingToolkitBase.jl | 3 +- .../src/systems/homotopy_sweep.jl | 129 ++++++++++++++++++ .../test/homotopy_sweep.jl | 99 ++++++++++++++ lib/ModelingToolkitBase/test/runtests.jl | 1 + 4 files changed, 231 insertions(+), 1 deletion(-) create mode 100644 lib/ModelingToolkitBase/src/systems/homotopy_sweep.jl create mode 100644 lib/ModelingToolkitBase/test/homotopy_sweep.jl diff --git a/lib/ModelingToolkitBase/src/ModelingToolkitBase.jl b/lib/ModelingToolkitBase/src/ModelingToolkitBase.jl index 3a71cb2124..cd8429a57a 100644 --- a/lib/ModelingToolkitBase/src/ModelingToolkitBase.jl +++ b/lib/ModelingToolkitBase/src/ModelingToolkitBase.jl @@ -195,6 +195,7 @@ include("systems/codegen_utils.jl") include("problems/docs.jl") include("systems/codegen.jl") include("systems/homotopy.jl") +include("systems/homotopy_sweep.jl") include("systems/problem_utils.jl") include("problems/compatibility.jl") @@ -271,7 +272,7 @@ export liouville_transform, change_independent_variable, export respecialize export PDESystem export Differential, expand_derivatives, @derivatives -export homotopy +export homotopy, HomotopySweep, TrivialHomotopy, TrivialThenSweep export Equation export Term export SymScope, LocalScope, ParentScope, GlobalScope diff --git a/lib/ModelingToolkitBase/src/systems/homotopy_sweep.jl b/lib/ModelingToolkitBase/src/systems/homotopy_sweep.jl new file mode 100644 index 0000000000..c2e26eba17 --- /dev/null +++ b/lib/ModelingToolkitBase/src/systems/homotopy_sweep.jl @@ -0,0 +1,129 @@ +using CommonSolve: CommonSolve +using SciMLBase: SciMLBase, NonlinearProblem, ReturnCode, remake, + successful_retcode, state_values, parameter_values +using Setfield: @set + +""" + HomotopySweep(; inner = NewtonRaphson(), schedule = 0.0:0.1:1.0, + set_λ!, maxiters_per_step = nothing) + +Parameter-sweep continuation algorithm for `NonlinearProblem` lowered by +`rewrite_with_lambda`. Each step `λ_k ∈ schedule` is solved by `inner` with +`u0` warm-started from the previous step's solution. λ is written into the +problem's parameter vector via `set_λ!`, normally a `SymbolicIndexingInterface.setp` +handle (or a `(p, v) -> new_p` closure in tests). + +Failure policy: if any step's inner solve returns an unsuccessful retcode, +`HomotopySweep` returns that step's solution unchanged. No auto-fallback. +Auto-fallback to trivial is the job of `TrivialThenSweep`, not this struct. +""" +struct HomotopySweep{Inner, Sched, Setter} <: SciMLBase.AbstractNonlinearAlgorithm + inner::Inner + schedule::Sched + set_λ!::Setter + maxiters_per_step::Union{Int, Nothing} +end + +function HomotopySweep(; inner = nothing, schedule = 0.0:0.1:1.0, + set_λ!, maxiters_per_step = nothing) + return HomotopySweep(inner === nothing ? _default_inner() : inner, + schedule, set_λ!, maxiters_per_step) +end + +# Lazy-load NewtonRaphson from NonlinearSolve, which is a test-extras dependency +# of ModelingToolkitBase (not a runtime [deps]). Callers can avoid this path by +# passing `inner = ...` explicitly; the default is only invoked when the test +# env (or a downstream user env) has already loaded NonlinearSolve. +function _default_inner() + nls = Base.get_extension(@__MODULE__, :NonlinearSolve) + if nls !== nothing + return nls.NewtonRaphson() + end + # Try via Base.require if NonlinearSolve is available in the current env. + try + mod = Base.require(Base.PkgId(Base.UUID("8913a72c-1f9b-4ce2-8d82-65094dcecaec"), + "NonlinearSolve")) + return mod.NewtonRaphson() + catch err + throw(ArgumentError( + "HomotopySweep/TrivialHomotopy default `inner` requires NonlinearSolve " * + "to be loaded. Either `using NonlinearSolve` first, or pass `inner = ...` " * + "explicitly. (Underlying error: $err)")) + end +end + +function CommonSolve.solve(prob::NonlinearProblem, alg::HomotopySweep; kwargs...) + u_curr = copy(state_values(prob)) + last_sol = nothing + for λ in alg.schedule + p_in = copy(parameter_values(prob)) + ret = alg.set_λ!(p_in, λ) + new_p = ret === nothing ? p_in : ret + step_prob = remake(prob; u0 = u_curr, p = new_p) + inner_kwargs = alg.maxiters_per_step === nothing ? + kwargs : (; kwargs..., maxiters = alg.maxiters_per_step) + step_sol = CommonSolve.solve(step_prob, alg.inner; inner_kwargs...) + last_sol = step_sol + if !successful_retcode(step_sol) + return step_sol + end + u_curr = copy(step_sol.u) + end + return last_sol +end + +""" + TrivialHomotopy(; inner = NewtonRaphson()) + +Trivial-form init: leaves `__homotopy_λ` at its default 1.0 (so the lowered +system is `actual`) and solves once with `inner`. Equivalent to OMC's +`-noHomotopyOnFirstTry`. Use when the guess is reliably close to the actual +root and sweep cost is unwarranted. +""" +struct TrivialHomotopy{Inner} <: SciMLBase.AbstractNonlinearAlgorithm + inner::Inner +end + +TrivialHomotopy(; inner = nothing) = + TrivialHomotopy(inner === nothing ? _default_inner() : inner) + +function CommonSolve.solve(prob::NonlinearProblem, alg::TrivialHomotopy; kwargs...) + return CommonSolve.solve(prob, alg.inner; kwargs...) +end + +""" + TrivialThenSweep(; trivial = TrivialHomotopy(), sweep = HomotopySweep(...)) + +Composite algorithm matching OpenModelica's default user experience: attempt +the trivial (single-Newton) solve first, and on unsuccessful retcode fall +back to a full parameter sweep. The returned `sol.original` records which +path succeeded under `:path` (`:trivial` or `:sweep_fallback`). + +This is MTK's default `initializealg.nlsolve` when a system contains +`homotopy(...)` nodes. Users override by passing `OverrideInit(nlsolve = ...)` +explicitly. +""" +struct TrivialThenSweep{T, S} <: SciMLBase.AbstractNonlinearAlgorithm + trivial::T + sweep::S +end + +TrivialThenSweep(; trivial = TrivialHomotopy(), sweep) = + TrivialThenSweep(trivial, sweep) + +function CommonSolve.solve(prob::NonlinearProblem, alg::TrivialThenSweep; kwargs...) + trivial_sol = CommonSolve.solve(prob, alg.trivial; kwargs...) + if successful_retcode(trivial_sol) + return _annotate_path(trivial_sol, :trivial) + end + sweep_sol = CommonSolve.solve(prob, alg.sweep; kwargs...) + return _annotate_path(sweep_sol, :sweep_fallback) +end + +# Attach a path marker to sol via the `original` field. The inner solver's +# previous `original` (if any) is preserved under `:inner` so callers that need +# the underlying object can still reach it. +function _annotate_path(sol, path::Symbol) + inner_original = hasproperty(sol, :original) ? getproperty(sol, :original) : nothing + return @set sol.original = (; path = path, inner = inner_original) +end diff --git a/lib/ModelingToolkitBase/test/homotopy_sweep.jl b/lib/ModelingToolkitBase/test/homotopy_sweep.jl new file mode 100644 index 0000000000..cf1120bb2e --- /dev/null +++ b/lib/ModelingToolkitBase/test/homotopy_sweep.jl @@ -0,0 +1,99 @@ +using Test +using ModelingToolkitBase +using ModelingToolkitBase: HomotopySweep, TrivialHomotopy, TrivialThenSweep +using SciMLBase +using NonlinearSolve: NewtonRaphson + +@testset "Homotopy nlsolve algorithms" begin + # Build the canonical homotopy fixture: f(u, p) = (1-λ)*(u-1) + λ*(u^2-4) + # - At λ=0: u = 1 (simplified). At λ=1: u = ±2 (actual). + # - From u0 = 1, sweep should track to u ≈ 2. + f_canonical! = (du, u, p) -> begin + λ = p[1] + du[1] = (1 - λ) * (u[1] - 1) + λ * (u[1]^2 - 4) + end + set_λ_explicit = (p, v) -> begin + p2 = copy(p) + p2[1] = v + p2 + end + + @testset "S1 HomotopySweep linear schedule converges to actual" begin + prob = NonlinearProblem(f_canonical!, [1.0], [0.0]) + alg = HomotopySweep(; inner = NewtonRaphson(), + schedule = 0.0:0.1:1.0, + set_λ! = set_λ_explicit) + sol = solve(prob, alg) + @test SciMLBase.successful_retcode(sol) + @test abs(sol.u[1] - 2.0) < 1e-6 + end + + @testset "S2 HomotopySweep step failure returns unsuccessful retcode" begin + # Force impossible homotopy: f = u^2 + 1, no real root → Newton diverges. + f_bad! = (du, u, p) -> begin + λ = p[1] + du[1] = (1 - λ) * u[1] + λ * (u[1]^2 + 1) + end + prob = NonlinearProblem(f_bad!, [0.5], [0.0]) + alg = HomotopySweep(; inner = NewtonRaphson(), + schedule = 0.0:0.1:1.0, + set_λ! = set_λ_explicit, + maxiters_per_step = 20) + sol = solve(prob, alg) + @test !SciMLBase.successful_retcode(sol) + end + + @testset "S3 HomotopySweep default constructor" begin + alg = HomotopySweep(; set_λ! = ((p, v) -> p)) + @test alg.schedule == 0.0:0.1:1.0 + # NonlinearSolve 5+ exposes `NewtonRaphson` as a constructor function + # (not a type), so compare against the type of an instance. + @test alg.inner isa typeof(NewtonRaphson()) + end + + @testset "S4 TrivialHomotopy dispatches to inner once" begin + f! = (du, u, p) -> du[1] = u[1]^2 - 4 + prob = NonlinearProblem(f!, [1.5], Float64[]) + alg = TrivialHomotopy() + sol = solve(prob, alg) + @test SciMLBase.successful_retcode(sol) + @test abs(sol.u[1] - 2.0) < 1e-6 + end + + @testset "S5 TrivialThenSweep succeeds via trivial path when guess is good" begin + # u0 = 1.9 is close to actual root u=2 → Newton converges from trivial. + # Sweep should NOT be invoked. We detect by `sol.original.path === :trivial`. + prob = NonlinearProblem(f_canonical!, [1.9], [1.0]) # p[1]=1.0 → already trivial + alg = TrivialThenSweep(; + trivial = TrivialHomotopy(), + sweep = HomotopySweep(; inner = NewtonRaphson(), + schedule = 0.0:0.1:1.0, + set_λ! = set_λ_explicit)) + sol = solve(prob, alg) + @test SciMLBase.successful_retcode(sol) + @test abs(sol.u[1] - 2.0) < 1e-6 + @test sol.original isa NamedTuple && sol.original.path === :trivial + end + + @testset "S6 TrivialThenSweep falls back to sweep when trivial fails" begin + # Fixture: simplified u (linear, root u=0), actual atan(u) (root u=0). + # Newton on atan(u) from u0=10 is the classic diverging case (stalls/diverges + # because |u_{n+1}| > |u_n| outside the basin of attraction). + # Sweep from u0=10 at λ=0 immediately drives u→0 via the simplified branch, + # then atan walks smoothly to root u=0 across the schedule. + f_hard! = (du, u, p) -> begin + λ = p[1] + du[1] = (1 - λ) * u[1] + λ * atan(u[1]) + end + prob = NonlinearProblem(f_hard!, [10.0], [1.0]) # start at λ=1 (trivial) + alg = TrivialThenSweep(; + trivial = TrivialHomotopy(; inner = NewtonRaphson()), + sweep = HomotopySweep(; inner = NewtonRaphson(), + schedule = 0.0:0.1:1.0, + set_λ! = set_λ_explicit)) + sol = solve(prob, alg) + @test SciMLBase.successful_retcode(sol) + @test abs(atan(sol.u[1])) < 1e-4 + @test sol.original.path === :sweep_fallback + end +end diff --git a/lib/ModelingToolkitBase/test/runtests.jl b/lib/ModelingToolkitBase/test/runtests.jl index 9eb123cafb..b9f5b9cbc5 100644 --- a/lib/ModelingToolkitBase/test/runtests.jl +++ b/lib/ModelingToolkitBase/test/runtests.jl @@ -65,6 +65,7 @@ end @safetestset "InitializationSystem Test" include("initializationsystem.jl") @safetestset "Homotopy Operator L0 Trivial" include("homotopy_init.jl") @safetestset "Homotopy lowering" begin include("homotopy_lowering.jl") end + @safetestset "Homotopy nlsolve algorithms" begin include("homotopy_sweep.jl") end @safetestset "Initial Values Test" include("initial_values.jl") end From aa9ce905c4eb17478b70b76769209247ca1155b5 Mon Sep 17 00:00:00 2001 From: lf Date: Mon, 25 May 2026 02:19:46 +0800 Subject: [PATCH 14/28] =?UTF-8?q?refactor:=20polish=20HomotopySweep=20?= =?UTF-8?q?=E2=80=94=20drop=20dead=20get=5Fextension=20branch,=20doc=20set?= =?UTF-8?q?p=20contract?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- lib/ModelingToolkitBase/src/systems/homotopy_sweep.jl | 11 +++++------ lib/ModelingToolkitBase/test/homotopy_sweep.jl | 3 +++ 2 files changed, 8 insertions(+), 6 deletions(-) diff --git a/lib/ModelingToolkitBase/src/systems/homotopy_sweep.jl b/lib/ModelingToolkitBase/src/systems/homotopy_sweep.jl index c2e26eba17..d83efb83de 100644 --- a/lib/ModelingToolkitBase/src/systems/homotopy_sweep.jl +++ b/lib/ModelingToolkitBase/src/systems/homotopy_sweep.jl @@ -1,5 +1,5 @@ using CommonSolve: CommonSolve -using SciMLBase: SciMLBase, NonlinearProblem, ReturnCode, remake, +using SciMLBase: SciMLBase, NonlinearProblem, remake, successful_retcode, state_values, parameter_values using Setfield: @set @@ -35,11 +35,6 @@ end # passing `inner = ...` explicitly; the default is only invoked when the test # env (or a downstream user env) has already loaded NonlinearSolve. function _default_inner() - nls = Base.get_extension(@__MODULE__, :NonlinearSolve) - if nls !== nothing - return nls.NewtonRaphson() - end - # Try via Base.require if NonlinearSolve is available in the current env. try mod = Base.require(Base.PkgId(Base.UUID("8913a72c-1f9b-4ce2-8d82-65094dcecaec"), "NonlinearSolve")) @@ -55,6 +50,10 @@ end function CommonSolve.solve(prob::NonlinearProblem, alg::HomotopySweep; kwargs...) u_curr = copy(state_values(prob)) last_sol = nothing + # `set_λ!` is expected to be a `SymbolicIndexingInterface.setp` handle in + # production (mutates `prob.p` in place via `setindex!`, returns nothing). + # Tests pass an OOP `(p, v) -> new_p` closure; both paths are supported + # by the `ret === nothing ? p_in : ret` normalization below. for λ in alg.schedule p_in = copy(parameter_values(prob)) ret = alg.set_λ!(p_in, λ) diff --git a/lib/ModelingToolkitBase/test/homotopy_sweep.jl b/lib/ModelingToolkitBase/test/homotopy_sweep.jl index cf1120bb2e..4e12a64657 100644 --- a/lib/ModelingToolkitBase/test/homotopy_sweep.jl +++ b/lib/ModelingToolkitBase/test/homotopy_sweep.jl @@ -63,6 +63,9 @@ using NonlinearSolve: NewtonRaphson @testset "S5 TrivialThenSweep succeeds via trivial path when guess is good" begin # u0 = 1.9 is close to actual root u=2 → Newton converges from trivial. # Sweep should NOT be invoked. We detect by `sol.original.path === :trivial`. + # NOTE: p[1] = 1.0 is intentional — TrivialHomotopy does NOT touch λ, + # so the canonical fixture must already be in actual-form. Don't "fix" + # this to 0.0; that would silently flip the trivial path to simplified. prob = NonlinearProblem(f_canonical!, [1.9], [1.0]) # p[1]=1.0 → already trivial alg = TrivialThenSweep(; trivial = TrivialHomotopy(), From c63e47ffa2ba6ec24bcbac06e9c94d0c40aa24a2 Mon Sep 17 00:00:00 2001 From: lf Date: Mon, 25 May 2026 12:08:57 +0800 Subject: [PATCH 15/28] refactor: always lower homotopy(a,s) in init hook --- .../src/problems/initializationproblem.jl | 27 ++++++++++++++----- 1 file changed, 20 insertions(+), 7 deletions(-) diff --git a/lib/ModelingToolkitBase/src/problems/initializationproblem.jl b/lib/ModelingToolkitBase/src/problems/initializationproblem.jl index 7056591350..3f801f0f44 100644 --- a/lib/ModelingToolkitBase/src/problems/initializationproblem.jl +++ b/lib/ModelingToolkitBase/src/problems/initializationproblem.jl @@ -82,15 +82,28 @@ All other keyword arguments are forwarded to the wrapped problem constructor. binds[get_iv(sys)::SymbolicT] = Symbolics.COMMON_ZERO @set! isys.bindings = ROSymmapT(binds) end - # L0 trivial homotopy rewrite (Modelica spec 3.7.4 trivial form). - # Only runs when `homotopy(...)` is actually present, so non-homotopy - # systems pay zero overhead. PR2 (L1 parameter sweep) will replace this - # branch with λ-lowering when an explicit HomotopySweep algorithm is in use. - # Applied before `mtkcompile` so that downstream structural transforms - # operate on the actual equations, not on opaque `homotopy(...)` calls. + # Modelica homotopy operator (spec 3.7.4). Lower every `homotopy(a, s)` + # node to `(1 - λ)*s + λ*a` and inject `__homotopy_λ` with default 1.0. + # At λ=1 the lowered system is numerically `actual` (trivial form); + # `HomotopySweep` walks λ from 0 → 1 to obtain `actual`'s root from a + # `simplified` starting point. Applied before `mtkcompile` so structural + # transforms see real equations, not opaque `homotopy(...)` calls. if has_homotopy_in_equations(equations(isys)) - new_eqs = rewrite_trivial_in_equations(equations(isys)) + new_eqs, homotopy_λ = rewrite_with_lambda_in_equations(equations(isys)) @set! isys.eqs = new_eqs + homotopy_λ_sym = unwrap(homotopy_λ) + current_ps = get_ps(isys) + if !any(p -> isequal(p, homotopy_λ_sym), current_ps) + @set! isys.ps = vcat(current_ps, homotopy_λ_sym) + end + # `@parameters __homotopy_λ = 1.0` carries the default in `Num` metadata, + # but pushing the symbol straight into `isys.ps` here bypasses the System + # constructor that normally copies metadata defaults into `bindings`. So + # we add the binding explicitly to avoid `MTKParameters` complaining + # about a missing value at problem-construction time. + binds = copy(parent(bindings(isys))) + binds[homotopy_λ_sym] = Symbolics.SConst(1.0) + @set! isys.bindings = ROSymmapT(binds) end if simplify_system From 01e0262dc3a1610c9e2e026825b37d975d9fba46 Mon Sep 17 00:00:00 2001 From: lf Date: Mon, 25 May 2026 19:00:34 +0800 Subject: [PATCH 16/28] =?UTF-8?q?feat:=20lower=20homotopy=20in=20complete(?= =?UTF-8?q?)=20so=20parent=20and=20init=20share=20=CE=BB=20identity?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Lift the homotopy lowering pass from the InitializationProblem hook up to complete(), mirroring how MTK already injects Initial(x) parameters. The parent system and init system now both carry the same `__homotopy_λ` parameter, which lets MTKParametersReconstructor resolve λ when it appears in init-system observed expressions (the bug that broke PR1's Buildings PressureDrop fixture under the previous attempt). - add_homotopy_parameter(sys) lowers equations + observed and injects λ; hooked into complete() right after add_initialization_parameters - InitializationProblem hook simplified to (a) plant op[λ]=1.0 for runtime default (avoids mtkcompile substituting λ away via sys.bindings) and (b) defensive fallback for systems built without complete() --- .../src/problems/initializationproblem.jl | 46 ++++++++----- .../src/systems/abstractsystem.jl | 8 +++ .../src/systems/homotopy.jl | 69 +++++++++++++++++++ 3 files changed, 105 insertions(+), 18 deletions(-) diff --git a/lib/ModelingToolkitBase/src/problems/initializationproblem.jl b/lib/ModelingToolkitBase/src/problems/initializationproblem.jl index 3f801f0f44..19a03194b5 100644 --- a/lib/ModelingToolkitBase/src/problems/initializationproblem.jl +++ b/lib/ModelingToolkitBase/src/problems/initializationproblem.jl @@ -82,28 +82,38 @@ All other keyword arguments are forwarded to the wrapped problem constructor. binds[get_iv(sys)::SymbolicT] = Symbolics.COMMON_ZERO @set! isys.bindings = ROSymmapT(binds) end - # Modelica homotopy operator (spec 3.7.4). Lower every `homotopy(a, s)` - # node to `(1 - λ)*s + λ*a` and inject `__homotopy_λ` with default 1.0. - # At λ=1 the lowered system is numerically `actual` (trivial form); - # `HomotopySweep` walks λ from 0 → 1 to obtain `actual`'s root from a - # `simplified` starting point. Applied before `mtkcompile` so structural - # transforms see real equations, not opaque `homotopy(...)` calls. - if has_homotopy_in_equations(equations(isys)) + # Modelica `homotopy(actual, simplified)` (spec 3.7.4) lowering happens + # in the parent system inside `complete()` via `add_homotopy_parameter` + # (`systems/homotopy.jl`). The init system here inherits the lowered + # equations + `__homotopy_λ` parameter from the parent. + # + # `MTKParameters` requires every parameter to have a value in `op` + # (sys.bindings is intentionally NOT used — that would let `mtkcompile` + # substitute λ away and break the sweep path). So we plant the runtime + # default `λ = 1.0` into `op` here, which is op-level metadata and not + # consulted by `mtkcompile`'s substitution. `HomotopySweep` overrides + # this at solve time via the SII setp handle. + homotopy_λ_idx = findfirst(get_ps(isys)) do p + sym = unwrap(p) + SymbolicUtils.issym(sym) && nameof(sym) === :__homotopy_λ + end + if homotopy_λ_idx !== nothing + homotopy_λ_sym = unwrap(get_ps(isys)[homotopy_λ_idx]) + if !haskey(op, homotopy_λ_sym) + op[homotopy_λ_sym] = 1.0 + end + elseif has_homotopy_in_equations(equations(isys)) + # Defensive fallback: caller skipped `complete()` (rare), do the + # lowering inline so `mtkcompile` doesn't see opaque homotopy nodes. new_eqs, homotopy_λ = rewrite_with_lambda_in_equations(equations(isys)) @set! isys.eqs = new_eqs homotopy_λ_sym = unwrap(homotopy_λ) - current_ps = get_ps(isys) - if !any(p -> isequal(p, homotopy_λ_sym), current_ps) - @set! isys.ps = vcat(current_ps, homotopy_λ_sym) + if !any(p -> isequal(p, homotopy_λ_sym), get_ps(isys)) + @set! isys.ps = vcat(get_ps(isys), homotopy_λ_sym) + end + if !haskey(op, homotopy_λ_sym) + op[homotopy_λ_sym] = 1.0 end - # `@parameters __homotopy_λ = 1.0` carries the default in `Num` metadata, - # but pushing the symbol straight into `isys.ps` here bypasses the System - # constructor that normally copies metadata defaults into `bindings`. So - # we add the binding explicitly to avoid `MTKParameters` complaining - # about a missing value at problem-construction time. - binds = copy(parent(bindings(isys))) - binds[homotopy_λ_sym] = Symbolics.SConst(1.0) - @set! isys.bindings = ROSymmapT(binds) end if simplify_system diff --git a/lib/ModelingToolkitBase/src/systems/abstractsystem.jl b/lib/ModelingToolkitBase/src/systems/abstractsystem.jl index aa7e6928f5..1576b8ed69 100644 --- a/lib/ModelingToolkitBase/src/systems/abstractsystem.jl +++ b/lib/ModelingToolkitBase/src/systems/abstractsystem.jl @@ -698,6 +698,14 @@ function complete( _unhack_sys = unhack_system(sys) if add_initial_parameters sys = add_initialization_parameters(sys; split, _unhack_sys)::T + # Modelica `homotopy(actual, simplified)` operator (spec 3.7.4): + # if the system carries any homotopy nodes, lower them now to + # `(1-λ)*s + λ*a` and inject `__homotopy_λ` into the parent so + # `MTKParametersReconstructor` can resolve λ in init-system + # observed expressions. Lives at the same lifecycle point as + # `add_initialization_parameters` so the index_cache rebuild at + # line 723+ picks up λ uniformly. + sys = add_homotopy_parameter(sys)::T end cb_alg_eqs = Equation[alg_equations(_unhack_sys); observed(_unhack_sys)] if has_continuous_events(sys) && is_time_dependent(sys) diff --git a/lib/ModelingToolkitBase/src/systems/homotopy.jl b/lib/ModelingToolkitBase/src/systems/homotopy.jl index 3ce95d5f2d..888cae2ba5 100644 --- a/lib/ModelingToolkitBase/src/systems/homotopy.jl +++ b/lib/ModelingToolkitBase/src/systems/homotopy.jl @@ -132,3 +132,72 @@ function rewrite_with_lambda_in_equations(eqs, λ = nothing) _rewrite_with_lambda(unwrap(eq.rhs), λ)) for eq in eqs] return new_eqs, λ end + +""" + has_any_homotopy(sys) + +Return `true` iff `sys` contains a `homotopy(...)` node anywhere in its +equations OR its observed equations. Used by `add_homotopy_parameter` to +decide whether to inject `__homotopy_λ` into a parent system at `complete` +time. +""" +function has_any_homotopy(sys) + has_homotopy_in_equations(equations(sys)) && return true + obs = observed(sys) + obs === nothing && return false + for eq in obs + (has_homotopy(eq.lhs) || has_homotopy(eq.rhs)) && return true + end + return false +end + +""" + add_homotopy_parameter(sys) + +If `sys` contains `homotopy(a, s)` nodes (in equations or observed), lower +every node to `(1 - λ)*s + λ*a` and inject a shared `__homotopy_λ` parameter +(default value `1.0`) into the system's parameter list. Mirrors the +`add_initialization_parameters` pattern (`abstractsystem.jl:559`) — same +lifecycle point inside `complete()`, so the parent and the init system end +up sharing one identity for `__homotopy_λ`. This lets downstream codepaths +like `MTKParametersReconstructor` resolve λ when it appears in init-system +observed expressions. + +At λ=1 (the default) the lowered expression reduces numerically to `actual`; +`HomotopySweep` walks λ from 0 → 1 to obtain the actual root from a +`simplified` starting point. + +No-op if `sys` doesn't support initialization (Modelica spec 3.7.4 only +applies during init), is itself an initializesystem, or contains no +`homotopy(...)` nodes. +""" +function add_homotopy_parameter(sys::AbstractSystem) + supports_initialization(sys) || return sys + is_initializesystem(sys) && return sys + has_any_homotopy(sys) || return sys + + λ = only(@parameters __homotopy_λ = 1.0) + λ_sym = unwrap(λ) + + # Lower homotopy nodes in equations + new_eqs, _ = rewrite_with_lambda_in_equations(equations(sys), λ_sym) + @set! sys.eqs = new_eqs + + # Lower homotopy nodes in observed equations (PressureDrop-style: m_flow + # is eliminated by `mtkcompile` and its definition lives in observed — + # downstream observed codegen must see lowered form, not opaque homotopy). + obs = observed(sys) + if obs !== nothing && !isempty(obs) + new_obs = [Equation(_rewrite_with_lambda(unwrap(eq.lhs), λ_sym), + _rewrite_with_lambda(unwrap(eq.rhs), λ_sym)) for eq in obs] + @set! sys.observed = new_obs + end + + # Inject λ as a parameter on the parent — idempotent + current_ps = get_ps(sys) + if !any(p -> isequal(p, λ_sym), current_ps) + @set! sys.ps = vcat(current_ps, λ_sym) + end + + return sys +end From 9441577963ad785ae4122520f83d34859ab36d7b Mon Sep 17 00:00:00 2001 From: lf Date: Mon, 25 May 2026 19:25:39 +0800 Subject: [PATCH 17/28] feat: stash setp handle + TrivialThenSweep default in OverrideInitData MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit When the init system contains __homotopy_λ (injected by add_homotopy_parameter in complete()), populate two new InitializationMetadata fields: - homotopy_set_λ! — a SymbolicIndexingInterface.setp wrapper that HomotopySweep uses to mutate λ between continuation steps - homotopy_default_initializealg — a complete OverrideInit(nlsolve= TrivialThenSweep(...)) instance for Task 5 to inject as the default initializealg when no homotopy nodes are present, both fields stay nothing and the metadata behaves exactly as before. Adds an L1-Q5 testset covering the round-trip and default-alg shape. --- .../src/systems/problem_utils.jl | 57 +++++++++++++++++++ .../test/homotopy_lowering.jl | 46 +++++++++++++++ 2 files changed, 103 insertions(+) diff --git a/lib/ModelingToolkitBase/src/systems/problem_utils.jl b/lib/ModelingToolkitBase/src/systems/problem_utils.jl index 6e1036c028..5d7ecebc41 100644 --- a/lib/ModelingToolkitBase/src/systems/problem_utils.jl +++ b/lib/ModelingToolkitBase/src/systems/problem_utils.jl @@ -1249,6 +1249,37 @@ struct InitializationMetadata{R <: ReconstructInitializeprob, GUU, SIU} The value of the `missing_guess_value` keyword indicating how to handle missing guesses. """ missing_guess_value::MissingGuessValue.Type + """ + Writer closure `(p, λ) -> p_new` that sets `__homotopy_λ` in a parameter + container. Built from `SymbolicIndexingInterface.setp` when the init + system carries a `__homotopy_λ` parameter (i.e. the homotopy lowering + in `complete()` fired); `nothing` otherwise. Consumed by `HomotopySweep` + during parameter-sweep continuation. + """ + homotopy_set_λ!::Any + """ + The default `OverrideInit(nlsolve = TrivialThenSweep(...))` instance MTK + can inject when the user does not supply `initializealg` and the system + contains homotopy nodes; `nothing` otherwise. Mirrors OpenModelica's + default user experience (try trivial first, fall back to sweep). + """ + homotopy_default_initializealg::Any +end + +# Backward-compatible 9-positional outer constructor: callers that predate +# the homotopy fields keep working; the new fields default to `nothing`. +function InitializationMetadata( + op::SymmapT, guesses::SymmapT, + additional_initialization_eqs::Vector{Equation}, + use_scc::Bool, time_dependent_init::Bool, + oop_reconstruct_u0_p::R, get_updated_u0::GUU, + set_initial_unknowns!::SIU, missing_guess_value::MissingGuessValue.Type, + ) where {R <: ReconstructInitializeprob, GUU, SIU} + return InitializationMetadata{R, GUU, SIU}( + op, guesses, additional_initialization_eqs, use_scc, time_dependent_init, + oop_reconstruct_u0_p, get_updated_u0, set_initial_unknowns!, missing_guess_value, + nothing, nothing, + ) end """ @@ -1482,6 +1513,32 @@ function maybe_build_initialization_problem( get_initial_unknowns, SetInitialUnknowns(sys), missing_guess_value ) + # If the init system carries `__homotopy_λ` (injected by the homotopy + # lowering pass in `complete()` via `add_homotopy_parameter`), populate + # the homotopy-related metadata so downstream solvers (HomotopySweep) + # can mutate λ during continuation and Task 5's default-injection hook + # can find the OMC-aligned `TrivialThenSweep` default. When the system + # has no homotopy nodes, both fields stay `nothing`. + initsys_for_homotopy = initializeprob.f.sys + homotopy_λ_idx = findfirst(parameters(initsys_for_homotopy)) do p + sym = SymbolicUtils.unwrap(p) + SymbolicUtils.issym(sym) && nameof(sym) === :__homotopy_λ + end + if homotopy_λ_idx !== nothing + raw_setter = SymbolicIndexingInterface.setp(initsys_for_homotopy, :__homotopy_λ) + # Normalize both mutating (returns Nothing) and OOP setp contracts. + set_λ! = (p, v) -> begin + p2 = copy(p) + ret = raw_setter(p2, v) + ret === nothing ? p2 : ret + end + default_sweep = HomotopySweep(; schedule = 0.0:0.1:1.0, set_λ! = set_λ!) + default_alg = SciMLBase.OverrideInit(; + nlsolve = TrivialThenSweep(; sweep = default_sweep)) + @set! meta.homotopy_set_λ! = set_λ! + @set! meta.homotopy_default_initializealg = default_alg + end + if time_dependent_init all_init_syms = Set(all_symbols(initializeprob)) solved_unknowns = filter(var -> var in all_init_syms, unknowns(sys)) diff --git a/lib/ModelingToolkitBase/test/homotopy_lowering.jl b/lib/ModelingToolkitBase/test/homotopy_lowering.jl index ce1944506e..c6008e2e1a 100644 --- a/lib/ModelingToolkitBase/test/homotopy_lowering.jl +++ b/lib/ModelingToolkitBase/test/homotopy_lowering.jl @@ -57,4 +57,50 @@ using Symbolics @test isequal(new_eqs[2].rhs, eqs[2].rhs) @test λ !== nothing end + + @testset "L1-Q5 OverrideInitData metadata exposes setp handle + default alg" begin + using ModelingToolkitBase: System, mtkcompile, HomotopySweep, TrivialThenSweep + using ModelingToolkitBase: t_nounits as t, D_nounits as D + using OrdinaryDiffEqRosenbrock: Rodas5P + using OrdinaryDiffEqNonlinearSolve + using SymbolicIndexingInterface + using SciMLBase + + @variables x(t) y(t) + @parameters p + eqs = [D(x) ~ -x, + 0 ~ homotopy(y^2 - p, y - 1)] + @named sys = System(eqs, t; guesses = [y => 2.5]) + sys = mtkcompile(sys) + prob = ODEProblem(sys, Dict(x => 1.0, p => 9.0), (0.0, 1.0)) + + meta = prob.f.initialization_data.metadata + @test meta.homotopy_set_λ! !== nothing + @test meta.homotopy_default_initializealg !== nothing + + # setp handle round-trip — write λ = 0.3, read it back + initprob = prob.f.initialization_data.initializeprob + set_λ! = meta.homotopy_set_λ! + new_p = set_λ!(parameter_values(initprob), 0.3) + getter = SymbolicIndexingInterface.getp(initprob, :__homotopy_λ) + @test isapprox(getter(new_p), 0.3; atol = 1e-12) + + # Default initializealg shape + default_alg = meta.homotopy_default_initializealg + @test default_alg isa SciMLBase.OverrideInit + @test default_alg.nlsolve isa TrivialThenSweep + @test default_alg.nlsolve.sweep isa HomotopySweep + + # Non-homotopy system: metadata fields stay `nothing`. + @variables a(t) + @parameters q + @named sys2 = System([D(a) ~ -q * a], t; guesses = [a => 1.0]) + sys2 = mtkcompile(sys2) + prob2 = ODEProblem(sys2, Dict(a => 1.0, q => 2.0), (0.0, 1.0)) + if prob2.f.initialization_data !== nothing + meta2 = prob2.f.initialization_data.metadata + @test meta2.homotopy_set_λ! === nothing + @test meta2.homotopy_default_initializealg === nothing + end + end end From 31943c764ad85a8af7df85341d45b30328f19f42 Mon Sep 17 00:00:00 2001 From: lf Date: Tue, 26 May 2026 13:52:10 +0800 Subject: [PATCH 18/28] feat: inject TrivialThenSweep as default initializealg in prob.kwargs MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit When add_homotopy_parameter has lowered a homotopy node, the init metadata exposes an OMC-aligned OverrideInit(TrivialThenSweep(...)) instance. process_kwargs now reads it (via a new initialization_data kwarg) and merges initializealg into prob.kwargs — but only when the caller did not pass initializealg themselves. Explicit user-supplied initializealg overrides via the standard merge(kwargs1, kwargs) order. Wired up for ODEProblem; other problem constructors can opt in by passing initialization_data = f.initialization_data to process_kwargs. Adds L1-Q6 testset covering auto-inject, explicit override, and non-homotopy-no-injection paths. --- .../src/problems/odeproblem.jl | 4 +- .../src/systems/problem_utils.jl | 20 +++++++++- .../test/homotopy_lowering.jl | 38 +++++++++++++++++++ 3 files changed, 60 insertions(+), 2 deletions(-) diff --git a/lib/ModelingToolkitBase/src/problems/odeproblem.jl b/lib/ModelingToolkitBase/src/problems/odeproblem.jl index 42d041208e..7f9e6fc1fe 100644 --- a/lib/ModelingToolkitBase/src/problems/odeproblem.jl +++ b/lib/ModelingToolkitBase/src/problems/odeproblem.jl @@ -116,7 +116,9 @@ Base.@nospecializeinfer @fallback_iip_specialize function SciMLBase.ODEProblem{i ) kwargs = process_kwargs( - sys; expression, callback, eval_expression, eval_module, op, _skip_events, tspan, kwargs... + sys; expression, callback, eval_expression, eval_module, op, _skip_events, tspan, + initialization_data = (hasproperty(f, :initialization_data) ? f.initialization_data : nothing), + kwargs... ) ptype = getmetadata(sys, ProblemTypeCtx, StandardODEProblem()) diff --git a/lib/ModelingToolkitBase/src/systems/problem_utils.jl b/lib/ModelingToolkitBase/src/systems/problem_utils.jl index 5d7ecebc41..4cfae663bd 100644 --- a/lib/ModelingToolkitBase/src/systems/problem_utils.jl +++ b/lib/ModelingToolkitBase/src/systems/problem_utils.jl @@ -2010,7 +2010,8 @@ end function process_kwargs( sys::System; expression = Val{false}, callback = nothing, eval_expression = false, eval_module = @__MODULE__, - _skip_events = false, _skip_tstops = false, tspan = nothing, kwargs... + _skip_events = false, _skip_tstops = false, tspan = nothing, + initialization_data = nothing, kwargs... ) kwargs = filter_kwargs(kwargs) kwargs1 = (;) @@ -2033,6 +2034,23 @@ function process_kwargs( end end + # Default `initializealg` injection for homotopy systems (Modelica spec + # 3.7.4). When `add_homotopy_parameter` lowered a homotopy node during + # `complete()`, the init metadata carries an OMC-aligned + # `OverrideInit(nlsolve = TrivialThenSweep(...))` instance. Inject it as + # the default `initializealg` ONLY if the user did not pass one + # explicitly — `solve(prob, alg; initializealg = X)` always wins because + # `prob.kwargs.initializealg` is overridden by an explicit solve-time + # kwarg per SciMLBase conventions. + if !haskey(kwargs, :initializealg) && + initialization_data isa SciMLBase.OverrideInitData && + initialization_data.metadata isa InitializationMetadata + default_alg = initialization_data.metadata.homotopy_default_initializealg + if default_alg !== nothing + kwargs1 = merge(kwargs1, (; initializealg = default_alg)) + end + end + return merge(kwargs1, kwargs) end diff --git a/lib/ModelingToolkitBase/test/homotopy_lowering.jl b/lib/ModelingToolkitBase/test/homotopy_lowering.jl index c6008e2e1a..967f420013 100644 --- a/lib/ModelingToolkitBase/test/homotopy_lowering.jl +++ b/lib/ModelingToolkitBase/test/homotopy_lowering.jl @@ -103,4 +103,42 @@ using Symbolics @test meta2.homotopy_default_initializealg === nothing end end + + @testset "L1-Q6 default initializealg injected into prob.kwargs for homotopy systems" begin + using ModelingToolkitBase: System, mtkcompile, TrivialThenSweep, TrivialHomotopy + using ModelingToolkitBase: t_nounits as t, D_nounits as D + using OrdinaryDiffEqRosenbrock + using OrdinaryDiffEqNonlinearSolve + using SciMLBase + + @variables x(t) y(t) + @parameters p + eqs = [D(x) ~ -x, + 0 ~ homotopy(y^2 - p, y - 1)] + @named sys = System(eqs, t; guesses = [y => 2.5]) + sys = mtkcompile(sys) + + # Auto-injection path: user passes no `initializealg` + prob_default = ODEProblem(sys, Dict(x => 1.0, p => 9.0), (0.0, 1.0)) + @test haskey(prob_default.kwargs, :initializealg) + injected = prob_default.kwargs[:initializealg] + @test injected isa SciMLBase.OverrideInit + @test injected.nlsolve isa TrivialThenSweep + + # Explicit override path: user passes their own initializealg + explicit_alg = SciMLBase.OverrideInit(; nlsolve = TrivialHomotopy()) + prob_explicit = ODEProblem(sys, Dict(x => 1.0, p => 9.0), (0.0, 1.0); + initializealg = explicit_alg) + @test haskey(prob_explicit.kwargs, :initializealg) + # User's explicit alg must win over MTK's default + @test prob_explicit.kwargs[:initializealg] === explicit_alg + + # Non-homotopy system: no `initializealg` injection + @variables a(t) + @parameters q + @named sys2 = System([D(a) ~ -q * a], t; guesses = [a => 1.0]) + sys2 = mtkcompile(sys2) + prob_nohomotopy = ODEProblem(sys2, Dict(a => 1.0, q => 2.0), (0.0, 1.0)) + @test !haskey(prob_nohomotopy.kwargs, :initializealg) + end end From 648aa6583e6d3479ad2f96cb902fc360012784fa Mon Sep 17 00:00:00 2001 From: lf Date: Tue, 26 May 2026 14:13:36 +0800 Subject: [PATCH 19/28] test: integration tests for default, explicit-sweep, explicit-trivial, default-injection wiring --- .../test/homotopy_integration.jl | 73 +++++++++++++++++++ lib/ModelingToolkitBase/test/runtests.jl | 1 + 2 files changed, 74 insertions(+) create mode 100644 lib/ModelingToolkitBase/test/homotopy_integration.jl diff --git a/lib/ModelingToolkitBase/test/homotopy_integration.jl b/lib/ModelingToolkitBase/test/homotopy_integration.jl new file mode 100644 index 0000000000..742bd27153 --- /dev/null +++ b/lib/ModelingToolkitBase/test/homotopy_integration.jl @@ -0,0 +1,73 @@ +using Test +using ModelingToolkitBase +using ModelingToolkitBase: System, mtkcompile, + HomotopySweep, TrivialHomotopy, TrivialThenSweep +using ModelingToolkitBase: t_nounits as t, D_nounits as D +using OrdinaryDiffEqRosenbrock: Rodas5P +using OrdinaryDiffEqNonlinearSolve +using SciMLBase +using NonlinearSolve: NewtonRaphson + +# Shared fixture: algebraic 0 = homotopy(y^2 - p, y - 1). +# Simplified at λ=0: y = 1. Actual at λ=1: y = ±√p. At p=9, guess y=2.5 → y=3. +function build_homotopy_fixture(; guess = 2.5, p_val = 9.0) + @variables x(t) y(t) + @parameters p + eqs = [D(x) ~ -x, + 0 ~ homotopy(y^2 - p, y - 1)] + @named sys = System(eqs, t; guesses = [y => guess]) + sys = mtkcompile(sys) + prob = ODEProblem(sys, Dict(x => 1.0, p => p_val), (0.0, 1.0)) + return prob, y +end + +@testset "homotopy operator — integration end-to-end" begin + @testset "I1 default (TrivialThenSweep) on good guess → trivial path wins" begin + prob, y = build_homotopy_fixture(; guess = 2.5) + sol = solve(prob, Rodas5P(); abstol = 1e-10, reltol = 1e-10) + @test SciMLBase.successful_retcode(sol) + @test abs(sol[y][1] - 3.0) < 1e-4 + end + + @testset "I2 explicit HomotopySweep" begin + prob, y = build_homotopy_fixture(; guess = 2.5) + meta = prob.f.initialization_data.metadata + sweep = HomotopySweep(; inner = NewtonRaphson(), + schedule = 0.0:0.1:1.0, + set_λ! = meta.homotopy_set_λ!) + initalg = SciMLBase.OverrideInit(; nlsolve = sweep) + sol = solve(prob, Rodas5P(); initializealg = initalg, + abstol = 1e-10, reltol = 1e-10) + @test SciMLBase.successful_retcode(sol) + @test abs(sol[y][1] - 3.0) < 1e-4 + end + + @testset "I3 explicit TrivialHomotopy" begin + prob, y = build_homotopy_fixture(; guess = 2.5) + initalg = SciMLBase.OverrideInit(; nlsolve = TrivialHomotopy()) + sol = solve(prob, Rodas5P(); initializealg = initalg, + abstol = 1e-10, reltol = 1e-10) + @test SciMLBase.successful_retcode(sol) + @test abs(sol[y][1] - 3.0) < 1e-4 + end + + @testset "I4 default-injection wires TrivialThenSweep all the way to solve" begin + # Structural assertion: the OMC-aligned default `TrivialThenSweep` + # really reaches `prob.kwargs.initializealg` so `solve` consumes it + # without explicit user opt-in. Trivial-vs-sweep dispatch within + # `TrivialThenSweep.solve` is tested directly in `homotopy_sweep.jl` + # (S5 trivial-success, S6 sweep-fallback) — no need to re-prove it + # through the ODEProblem pipeline. + prob, y = build_homotopy_fixture(; guess = 2.5) + @test haskey(prob.kwargs, :initializealg) + injected = prob.kwargs[:initializealg] + @test injected isa SciMLBase.OverrideInit + @test injected.nlsolve isa TrivialThenSweep + @test injected.nlsolve.trivial isa TrivialHomotopy + @test injected.nlsolve.sweep isa HomotopySweep + # And the actual end-to-end solve must succeed using only this default + sol = solve(prob, Rodas5P(); abstol = 1e-10, reltol = 1e-10) + @test SciMLBase.successful_retcode(sol) + @test abs(sol[y][1] - 3.0) < 1e-4 + end +end diff --git a/lib/ModelingToolkitBase/test/runtests.jl b/lib/ModelingToolkitBase/test/runtests.jl index b9f5b9cbc5..a4f9e03578 100644 --- a/lib/ModelingToolkitBase/test/runtests.jl +++ b/lib/ModelingToolkitBase/test/runtests.jl @@ -66,6 +66,7 @@ end @safetestset "Homotopy Operator L0 Trivial" include("homotopy_init.jl") @safetestset "Homotopy lowering" begin include("homotopy_lowering.jl") end @safetestset "Homotopy nlsolve algorithms" begin include("homotopy_sweep.jl") end + @safetestset "Homotopy integration" begin include("homotopy_integration.jl") end @safetestset "Initial Values Test" include("initial_values.jl") end From 4ff664a83681f481617cfac38db02d6443c11e03 Mon Sep 17 00:00:00 2001 From: lf Date: Tue, 26 May 2026 14:35:55 +0800 Subject: [PATCH 20/28] test: OMC numerical parity for Buildings PressureDrop fixture MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Frozen OMC reference: m_flow(0) = sqrt(5) ≈ 2.2360679774997896, x(0) = 1.0 (captured 2026-05-26 from OMC v1.27.0-dev via the modelica-OMEdit sandbox). All three nlsolve paths — default TrivialThenSweep, explicit HomotopySweep, explicit TrivialHomotopy — converge to the OMC value within 1e-6. The .mo fixture lives in the OMC sandbox at ~/code/modelica-OMEdit/homotopy-tests/PressureDropParityFixture.mo and is NOT committed to this repo per CLAUDE.md (no .mo files in workspace). --- .../test/homotopy_omc_parity.jl | 60 +++++++++++++++++++ lib/ModelingToolkitBase/test/runtests.jl | 1 + 2 files changed, 61 insertions(+) create mode 100644 lib/ModelingToolkitBase/test/homotopy_omc_parity.jl diff --git a/lib/ModelingToolkitBase/test/homotopy_omc_parity.jl b/lib/ModelingToolkitBase/test/homotopy_omc_parity.jl new file mode 100644 index 0000000000..1cfedf7924 --- /dev/null +++ b/lib/ModelingToolkitBase/test/homotopy_omc_parity.jl @@ -0,0 +1,60 @@ +using Test +using ModelingToolkitBase +using ModelingToolkitBase: System, mtkcompile, HomotopySweep, TrivialHomotopy +using ModelingToolkitBase: t_nounits as t, D_nounits as D +using OrdinaryDiffEqRosenbrock: Rodas5P +using OrdinaryDiffEqNonlinearSolve +using SciMLBase +using NonlinearSolve: NewtonRaphson + +@testset "homotopy OMC parity — Buildings PressureDrop fixture" begin + # Reference values captured 2026-05-26 from OMC v1.27.0-dev running + # `~/code/modelica-OMEdit/homotopy-tests/PressureDropParityFixture.mo` + # via the sandbox's `run_mo.sh`. m_flow(0) is √5 (the actual root); + # x(0) is the user-supplied start value, untouched by init. + OMC_M_FLOW_AT_T0 = 2.2360679774997896 # = sqrt(5) + OMC_X_AT_T0 = 1.0 + PARITY_TOL = 1e-6 + + @variables x(t) m_flow(t) + @parameters dp k m_flow_nominal dp_nominal x_decay + + basic_flow(dp_val, k_val) = sign(dp_val) * sqrt(abs(dp_val)) * k_val + + eqs = [D(x) ~ -x_decay * x, + 0 ~ m_flow - homotopy( + basic_flow(dp, k), + m_flow_nominal * dp / dp_nominal)] + @named sys = System(eqs, t; guesses = [m_flow => 2.0]) + sys = mtkcompile(sys) + op = Dict(x => 1.0, dp => 5.0, k => 1.0, + m_flow_nominal => 1.0, dp_nominal => 5.0, x_decay => 1.0) + prob = ODEProblem(sys, op, (0.0, 1.0)) + + @testset "default (TrivialThenSweep) matches OMC" begin + sol = solve(prob, Rodas5P(); abstol = 1e-12, reltol = 1e-12) + @test SciMLBase.successful_retcode(sol) + @test abs(sol[m_flow][1] - OMC_M_FLOW_AT_T0) < PARITY_TOL + @test abs(sol[x][1] - OMC_X_AT_T0) < PARITY_TOL + end + + @testset "explicit HomotopySweep matches OMC" begin + meta = prob.f.initialization_data.metadata + sweep = HomotopySweep(; inner = NewtonRaphson(), + schedule = 0.0:0.1:1.0, + set_λ! = meta.homotopy_set_λ!) + initalg = SciMLBase.OverrideInit(; nlsolve = sweep) + sol = solve(prob, Rodas5P(); initializealg = initalg, + abstol = 1e-12, reltol = 1e-12) + @test SciMLBase.successful_retcode(sol) + @test abs(sol[m_flow][1] - OMC_M_FLOW_AT_T0) < PARITY_TOL + end + + @testset "explicit TrivialHomotopy matches OMC" begin + initalg = SciMLBase.OverrideInit(; nlsolve = TrivialHomotopy()) + sol = solve(prob, Rodas5P(); initializealg = initalg, + abstol = 1e-12, reltol = 1e-12) + @test SciMLBase.successful_retcode(sol) + @test abs(sol[m_flow][1] - OMC_M_FLOW_AT_T0) < PARITY_TOL + end +end diff --git a/lib/ModelingToolkitBase/test/runtests.jl b/lib/ModelingToolkitBase/test/runtests.jl index a4f9e03578..0d01f30dc9 100644 --- a/lib/ModelingToolkitBase/test/runtests.jl +++ b/lib/ModelingToolkitBase/test/runtests.jl @@ -67,6 +67,7 @@ end @safetestset "Homotopy lowering" begin include("homotopy_lowering.jl") end @safetestset "Homotopy nlsolve algorithms" begin include("homotopy_sweep.jl") end @safetestset "Homotopy integration" begin include("homotopy_integration.jl") end + @safetestset "Homotopy OMC parity" begin include("homotopy_omc_parity.jl") end @safetestset "Initial Values Test" include("initial_values.jl") end From f98b3be84938abdb9e4cda7d9c546d0a823e9b48 Mon Sep 17 00:00:00 2001 From: lf Date: Mon, 1 Jun 2026 10:15:24 +0800 Subject: [PATCH 21/28] test: end-to-end sweep-rescue from out-of-basin guess (TrivialThenSweep continuation) --- .../test/homotopy_integration.jl | 24 +++++++++++++++++++ 1 file changed, 24 insertions(+) diff --git a/lib/ModelingToolkitBase/test/homotopy_integration.jl b/lib/ModelingToolkitBase/test/homotopy_integration.jl index 742bd27153..390b45f4cb 100644 --- a/lib/ModelingToolkitBase/test/homotopy_integration.jl +++ b/lib/ModelingToolkitBase/test/homotopy_integration.jl @@ -70,4 +70,28 @@ end @test SciMLBase.successful_retcode(sol) @test abs(sol[y][1] - 3.0) < 1e-4 end + + @testset "I5 out-of-basin guess: continuation rescues init the trivial solve can't" begin + # The sweep's reason-for-being, exercised THROUGH the real ODEProblem + # pipeline (homotopy_sweep.jl S6 only proves it for a hand-built + # NonlinearProblem). `actual = atan(y - 3)` has its single root at y = 3 + # but Newton diverges from any guess outside |y-3| ≲ 1.39 (the classic + # atan basin-escape). `simplified = y` has root y = 0. The two roots + # differ, so the landing value distinguishes "continuation tracked + # λ:0→1 onto the actual root" (y = 3) from "stuck at simplified" (y = 0). + @variables x(t) y(t) + eqs = [D(x) ~ -x, + 0 ~ homotopy(atan(y - 3.0), y)] + # guess y = 12 is far outside the actual root's Newton basin: a single + # Newton at λ=1 diverges, so reaching y = 3 is only possible via the + # default TrivialThenSweep continuation walking λ from 0 to 1. + @named sys = System(eqs, t; guesses = [y => 12.0]) + sys = mtkcompile(sys) + + prob = ODEProblem(sys, Dict(x => 1.0), (0.0, 1.0)) + sol = solve(prob, Rodas5P(); abstol = 1e-10, reltol = 1e-10) + @test SciMLBase.successful_retcode(sol) + @test abs(sol[y][1] - 3.0) < 1e-6 # actual root via the sweep + @test abs(sol[y][1]) > 0.5 # definitively NOT the simplified root y = 0 + end end From 6c5181a59bd20fd3b1e7b3de588e62bb5ae0316f Mon Sep 17 00:00:00 2001 From: lf Date: Mon, 1 Jun 2026 10:40:40 +0800 Subject: [PATCH 22/28] fix: homotopy inner solves must not re-trigger DAE init MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit An explicit `initializealg = OverrideInit(...)` on the parent problem is baked into the init problem's own `prob.kwargs`; `solve_up` merges it back into the inner nonlinear solve even if dropped from the forwarded call kwargs. The inner solver then runs NonlinearSolveBase's `OverrideInit` init path on a problem with `initialization_data === nothing` (that path is unguarded, unlike the default init path) and throws a `FieldError` on the divergence/reinit path — instead of surfacing a clean InitialFailure retcode. Force `initializealg = NoInit()` on the inner `TrivialHomotopy`/`HomotopySweep` solves so the call kwarg overrides any baked-in OverrideInit. Add regression test I6 (explicit TrivialHomotopy on an unsolvable init fails cleanly). --- .../src/systems/homotopy_sweep.jl | 17 ++++++++++++++-- .../test/homotopy_integration.jl | 20 +++++++++++++++++++ 2 files changed, 35 insertions(+), 2 deletions(-) diff --git a/lib/ModelingToolkitBase/src/systems/homotopy_sweep.jl b/lib/ModelingToolkitBase/src/systems/homotopy_sweep.jl index d83efb83de..b0e9e5dbfb 100644 --- a/lib/ModelingToolkitBase/src/systems/homotopy_sweep.jl +++ b/lib/ModelingToolkitBase/src/systems/homotopy_sweep.jl @@ -3,6 +3,18 @@ using SciMLBase: SciMLBase, NonlinearProblem, remake, successful_retcode, state_values, parameter_values using Setfield: @set +# An inner residual/continuation solve must never re-trigger DAE initialization. +# When the parent is built with an explicit `initializealg`, that `OverrideInit` +# is baked into the init problem's own `prob.kwargs` and `solve_up` merges it back +# even if we drop it from the forwarded call kwargs. NonlinearSolveBase then runs +# its `OverrideInit` init path on the inner problem — whose `f.initialization_data` +# is `nothing` — and throws a `FieldError` (that path is unguarded, unlike its +# default init path). Force `initializealg = NoInit()` on the inner solve so the +# call kwarg overrides any baked-in `OverrideInit` and no inner re-init runs. +_inner_solve_kwargs(kwargs) = (; + Base.structdiff((; kwargs...), NamedTuple{(:initializealg,)})..., + initializealg = SciMLBase.NoInit()) + """ HomotopySweep(; inner = NewtonRaphson(), schedule = 0.0:0.1:1.0, set_λ!, maxiters_per_step = nothing) @@ -59,8 +71,9 @@ function CommonSolve.solve(prob::NonlinearProblem, alg::HomotopySweep; kwargs... ret = alg.set_λ!(p_in, λ) new_p = ret === nothing ? p_in : ret step_prob = remake(prob; u0 = u_curr, p = new_p) + base_kwargs = _inner_solve_kwargs(kwargs) inner_kwargs = alg.maxiters_per_step === nothing ? - kwargs : (; kwargs..., maxiters = alg.maxiters_per_step) + base_kwargs : (; base_kwargs..., maxiters = alg.maxiters_per_step) step_sol = CommonSolve.solve(step_prob, alg.inner; inner_kwargs...) last_sol = step_sol if !successful_retcode(step_sol) @@ -87,7 +100,7 @@ TrivialHomotopy(; inner = nothing) = TrivialHomotopy(inner === nothing ? _default_inner() : inner) function CommonSolve.solve(prob::NonlinearProblem, alg::TrivialHomotopy; kwargs...) - return CommonSolve.solve(prob, alg.inner; kwargs...) + return CommonSolve.solve(prob, alg.inner; _inner_solve_kwargs(kwargs)...) end """ diff --git a/lib/ModelingToolkitBase/test/homotopy_integration.jl b/lib/ModelingToolkitBase/test/homotopy_integration.jl index 390b45f4cb..6b9400bf6f 100644 --- a/lib/ModelingToolkitBase/test/homotopy_integration.jl +++ b/lib/ModelingToolkitBase/test/homotopy_integration.jl @@ -94,4 +94,24 @@ end @test abs(sol[y][1] - 3.0) < 1e-6 # actual root via the sweep @test abs(sol[y][1]) > 0.5 # definitively NOT the simplified root y = 0 end + + @testset "I6 explicit TrivialHomotopy on an unsolvable init fails cleanly" begin + # Regression: an init the single-Newton trivial path cannot solve must + # surface as an unsuccessful retcode, NOT a cryptic internal error. The + # DAE-init driver forwards its `initializealg` (an OverrideInit) into the + # nlsolve `solve`; if that leaks into the inner solver it re-runs the + # OverrideInit on a problem with no initialization_data and throws. Same + # out-of-basin fixture as I5, but pinning the trivial path (no sweep + # fallback) so init genuinely cannot converge. + @variables x(t) y(t) + eqs = [D(x) ~ -x, + 0 ~ homotopy(atan(y - 3.0), y)] + @named sys = System(eqs, t; guesses = [y => 12.0]) + sys = mtkcompile(sys) + + prob = ODEProblem(sys, Dict(x => 1.0), (0.0, 1.0); + initializealg = SciMLBase.OverrideInit(nlsolve = TrivialHomotopy())) + sol = solve(prob, Rodas5P()) # must NOT throw an internal error + @test !SciMLBase.successful_retcode(sol) # clean init failure instead + end end From 4111bce71de448eb5983eabc97a4193e887ab2ad Mon Sep 17 00:00:00 2001 From: lf Date: Mon, 1 Jun 2026 10:47:52 +0800 Subject: [PATCH 23/28] fix: default homotopy inner solver falls back to SimpleNonlinearSolve NonlinearSolve is a test-only dependency of ModelingToolkitBase, not a runtime [deps]. `_default_inner` previously `Base.require`d it and threw if it was not loaded, so merely *constructing* a homotopy `ODEProblem` (which eagerly builds the default `TrivialThenSweep`) failed for any downstream user who had not run `using NonlinearSolve`. Fall back to `SimpleNonlinearSolve.SimpleNewtonRaphson` (a runtime [deps]) so homotopy systems build and solve out of the box; loading NonlinearSolve still restores its richer `NewtonRaphson` as the default. --- .../src/systems/homotopy_sweep.jl | 15 ++++++++++----- 1 file changed, 10 insertions(+), 5 deletions(-) diff --git a/lib/ModelingToolkitBase/src/systems/homotopy_sweep.jl b/lib/ModelingToolkitBase/src/systems/homotopy_sweep.jl index b0e9e5dbfb..cbd73c69c2 100644 --- a/lib/ModelingToolkitBase/src/systems/homotopy_sweep.jl +++ b/lib/ModelingToolkitBase/src/systems/homotopy_sweep.jl @@ -2,6 +2,7 @@ using CommonSolve: CommonSolve using SciMLBase: SciMLBase, NonlinearProblem, remake, successful_retcode, state_values, parameter_values using Setfield: @set +using SimpleNonlinearSolve: SimpleNewtonRaphson # An inner residual/continuation solve must never re-trigger DAE initialization. # When the parent is built with an explicit `initializealg`, that `OverrideInit` @@ -51,11 +52,15 @@ function _default_inner() mod = Base.require(Base.PkgId(Base.UUID("8913a72c-1f9b-4ce2-8d82-65094dcecaec"), "NonlinearSolve")) return mod.NewtonRaphson() - catch err - throw(ArgumentError( - "HomotopySweep/TrivialHomotopy default `inner` requires NonlinearSolve " * - "to be loaded. Either `using NonlinearSolve` first, or pass `inner = ...` " * - "explicitly. (Underlying error: $err)")) + catch + # NonlinearSolve is a test-only dependency of ModelingToolkitBase, not a + # runtime [deps]. When it isn't available, fall back to SimpleNonlinearSolve + # (a runtime [deps]) so a system carrying `homotopy(...)` can be built and + # solved out of the box — without this, merely constructing a homotopy + # `ODEProblem` would throw at problem-construction time. Load NonlinearSolve + # (`using NonlinearSolve`) to restore its richer `NewtonRaphson` as the + # default, or pass `inner = ...` explicitly. + return SimpleNewtonRaphson() end end From 3d2d5ed4a24e036e66506dbb595a9158ff3f70b6 Mon Sep 17 00:00:00 2001 From: lf Date: Mon, 1 Jun 2026 10:47:52 +0800 Subject: [PATCH 24/28] docs: document the Modelica homotopy operator (NEWS + API page) --- NEWS.md | 13 +++++++++++ docs/pages.jl | 1 + docs/src/API/homotopy.md | 48 ++++++++++++++++++++++++++++++++++++++++ 3 files changed, 62 insertions(+) create mode 100644 docs/src/API/homotopy.md diff --git a/NEWS.md b/NEWS.md index c1f09df1a8..89862a1de2 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,5 +1,18 @@ # ModelingToolkit v11 Release Notes +## Modelica `homotopy(actual, simplified)` operator + +ModelingToolkit now supports Modelica's `homotopy(actual, simplified)` operator +(Modelica specification §3.7.4) as an initialization aid. Annotate a hard-to-solve +relation with an easy `simplified` form, and initialization can start from the +`simplified` problem and continuously deform it into the `actual` one. Each node is +lowered to `(1 - λ)*simplified + λ*actual` with a shared `__homotopy_λ` parameter; +for systems containing `homotopy` nodes the default initialization becomes +`OverrideInit(nlsolve = TrivialThenSweep(...))`, which tries the trivial solve at +`λ = 1` first and falls back to a `λ`-sweep continuation (mirroring OpenModelica). +`HomotopySweep`, `TrivialHomotopy` and `TrivialThenSweep` are exported for explicit +use via `initializealg`. See the "Homotopy continuation for initialization" API page. + ## Symbolics@7 and SymbolicUtils@4 compatibility SymbolicUtils version 4 involved a major overhaul of the core symbolic infrastructure, which diff --git a/docs/pages.jl b/docs/pages.jl index 4d01d67324..8e8a68ccbe 100644 --- a/docs/pages.jl +++ b/docs/pages.jl @@ -40,6 +40,7 @@ pages = [ "API/variables.md", "API/model_building.md", "API/problems.md", + "API/homotopy.md", "API/dynamic_opt.md", "API/codegen.md", "API/PDESystem.md", diff --git a/docs/src/API/homotopy.md b/docs/src/API/homotopy.md new file mode 100644 index 0000000000..7f8f35114e --- /dev/null +++ b/docs/src/API/homotopy.md @@ -0,0 +1,48 @@ +```@meta +CollapsedDocStrings = true +``` + +# Homotopy continuation for initialization (Modelica `homotopy`) + +ModelingToolkit implements Modelica's `homotopy(actual, simplified)` operator +(Modelica specification §3.7.4). It is an initialization aid: where the `actual` +equations are hard to solve from the available guesses, the operator lets +initialization start from an easy `simplified` problem and continuously deform it +into the `actual` one. + +During `complete`/`mtkcompile`, every `homotopy(actual, simplified)` node is lowered +to + +```math +(1 - \lambda) \cdot \mathrm{simplified} + \lambda \cdot \mathrm{actual} +``` + +with a single shared parameter `__homotopy_λ` (default `1.0`). At `λ = 1` the system +reduces to `actual`; at `λ = 0` to `simplified`. For systems that contain `homotopy` +nodes, the default initialization algorithm becomes +`OverrideInit(nlsolve = TrivialThenSweep(...))`: it first attempts the trivial +single solve at `λ = 1` and, if that fails, runs a parameter-sweep continuation that +walks `λ` from 0 to 1 (mirroring OpenModelica's default user experience). Pass an +explicit `initializealg` to the problem constructor or to `solve` to override this. + +At runtime, outside initialization, `homotopy(actual, simplified)` evaluates to +`actual`, as required by the specification. + +!!! note + This operator is unrelated to the polynomial `HomotopyContinuationProblem`, which + solves polynomial systems via homotopy continuation. They only share the word + "homotopy". + +## Operator + +```@docs +ModelingToolkit.homotopy +``` + +## Initialization algorithms + +```@docs +ModelingToolkit.TrivialThenSweep +ModelingToolkit.HomotopySweep +ModelingToolkit.TrivialHomotopy +``` From 1bf8b4a8054f70e16df93b6322a13ca5598b6838 Mon Sep 17 00:00:00 2001 From: lf Date: Mon, 1 Jun 2026 11:26:21 +0800 Subject: [PATCH 25/28] fix: auto-inject homotopy default initializealg for non-ODE problems PIPE-1: the OMC-aligned default `OverrideInit(nlsolve = TrivialThenSweep)` was only threaded into `process_kwargs` by odeproblem.jl, so a homotopy System lowered to a NonlinearProblem / SteadyStateProblem / DAEProblem / SDEProblem never received the default initializealg and continuation did not engage without explicit user opt-in. Thread `initialization_data = f.initialization_data` into `process_kwargs` for those four constructors, mirroring odeproblem.jl, so the homotopy default reaches `prob.kwargs[:initializealg]` uniformly. Test I7 asserts the injection for a NonlinearProblem and a SteadyStateProblem. --- .../src/problems/daeproblem.jl | 5 +++- .../src/problems/nonlinearproblem.jl | 10 ++++++- .../src/problems/odeproblem.jl | 7 ++++- .../src/problems/sdeproblem.jl | 5 +++- .../test/homotopy_integration.jl | 30 +++++++++++++++++++ 5 files changed, 53 insertions(+), 4 deletions(-) diff --git a/lib/ModelingToolkitBase/src/problems/daeproblem.jl b/lib/ModelingToolkitBase/src/problems/daeproblem.jl index 872ba19efc..0869be8139 100644 --- a/lib/ModelingToolkitBase/src/problems/daeproblem.jl +++ b/lib/ModelingToolkitBase/src/problems/daeproblem.jl @@ -86,7 +86,10 @@ end kwargs = process_kwargs( sys; expression, callback, eval_expression, eval_module, - op, tspan, kwargs... + op, tspan, + initialization_data = (hasproperty(f, :initialization_data) ? f.initialization_data : + nothing), + kwargs... ) diffvars = collect_differential_variables(sys) diff --git a/lib/ModelingToolkitBase/src/problems/nonlinearproblem.jl b/lib/ModelingToolkitBase/src/problems/nonlinearproblem.jl index 7cae109703..acfb82cb52 100644 --- a/lib/ModelingToolkitBase/src/problems/nonlinearproblem.jl +++ b/lib/ModelingToolkitBase/src/problems/nonlinearproblem.jl @@ -78,7 +78,15 @@ end check_length, check_compatibility, expression, kwargs... ) - kwargs = process_kwargs(sys; kwargs...) + # Thread `initialization_data` so the homotopy default `initializealg` + # (OMC-aligned `OverrideInit(nlsolve = TrivialThenSweep)`) is auto-injected + # for homotopy systems lowered to a `NonlinearProblem` (PIPE-1). + kwargs = process_kwargs( + sys; + initialization_data = (hasproperty(f, :initialization_data) ? f.initialization_data : + nothing), + kwargs... + ) ptype = getmetadata(sys, ProblemTypeCtx, StandardNonlinearProblem()) args = (; f, u0, p, ptype) diff --git a/lib/ModelingToolkitBase/src/problems/odeproblem.jl b/lib/ModelingToolkitBase/src/problems/odeproblem.jl index 7f9e6fc1fe..fa97c78216 100644 --- a/lib/ModelingToolkitBase/src/problems/odeproblem.jl +++ b/lib/ModelingToolkitBase/src/problems/odeproblem.jl @@ -141,7 +141,12 @@ end is_steadystateprob = true, kwargs... ) - kwargs = process_kwargs(sys; expression, tspan = (0, Inf), kwargs...) + kwargs = process_kwargs( + sys; expression, tspan = (0, Inf), + initialization_data = (hasproperty(f, :initialization_data) ? f.initialization_data : + nothing), + kwargs... + ) args = (; f, u0, p) maybe_codegen_scimlproblem(expression, SteadyStateProblem{_iip}, args; kwargs...) diff --git a/lib/ModelingToolkitBase/src/problems/sdeproblem.jl b/lib/ModelingToolkitBase/src/problems/sdeproblem.jl index 6325def1bf..c10f38f502 100644 --- a/lib/ModelingToolkitBase/src/problems/sdeproblem.jl +++ b/lib/ModelingToolkitBase/src/problems/sdeproblem.jl @@ -108,7 +108,10 @@ end kwargs = process_kwargs( sys; expression, callback, eval_expression, eval_module, - op, _skip_events, tspan, kwargs... + op, _skip_events, tspan, + initialization_data = (hasproperty(f, :initialization_data) ? f.initialization_data : + nothing), + kwargs... ) args = (; f, u0, tspan, p) diff --git a/lib/ModelingToolkitBase/test/homotopy_integration.jl b/lib/ModelingToolkitBase/test/homotopy_integration.jl index 6b9400bf6f..5499548787 100644 --- a/lib/ModelingToolkitBase/test/homotopy_integration.jl +++ b/lib/ModelingToolkitBase/test/homotopy_integration.jl @@ -95,6 +95,36 @@ end @test abs(sol[y][1]) > 0.5 # definitively NOT the simplified root y = 0 end + @testset "I7 default injection reaches NonlinearProblem and SteadyStateProblem" begin + # PIPE-1: the OMC-aligned default `OverrideInit(nlsolve = TrivialThenSweep)` + # was only threaded into `process_kwargs` by `odeproblem.jl`. A homotopy + # System lowered to a `NonlinearProblem` / `SteadyStateProblem` must get the + # same default in `prob.kwargs[:initializealg]` so continuation engages + # without explicit user opt-in. + + # NonlinearProblem: time-independent algebraic homotopy system. + @variables y + @parameters p + @named nlsys = System([0 ~ homotopy(y^2 - p, y - 1)]; guesses = [y => 2.5]) + nlsys = mtkcompile(nlsys) + nlprob = NonlinearProblem(nlsys, Dict(y => 2.5, p => 9.0)) + @test haskey(nlprob.kwargs, :initializealg) + @test nlprob.kwargs[:initializealg] isa SciMLBase.OverrideInit + @test nlprob.kwargs[:initializealg].nlsolve isa TrivialThenSweep + + # SteadyStateProblem: time-dependent ODE homotopy system. + @variables x(t) yy(t) + @parameters q + eqs = [D(x) ~ -x, + 0 ~ homotopy(yy^2 - q, yy - 1)] + @named ssys = System(eqs, t; guesses = [yy => 2.5]) + ssys = mtkcompile(ssys) + ssprob = SteadyStateProblem(ssys, Dict(x => 1.0, q => 9.0)) + @test haskey(ssprob.kwargs, :initializealg) + @test ssprob.kwargs[:initializealg] isa SciMLBase.OverrideInit + @test ssprob.kwargs[:initializealg].nlsolve isa TrivialThenSweep + end + @testset "I6 explicit TrivialHomotopy on an unsolvable init fails cleanly" begin # Regression: an init the single-Newton trivial path cannot solve must # surface as an unsuccessful retcode, NOT a cryptic internal error. The From f9838fb0005848fea78d88bce8e8686d1b23be41 Mon Sep 17 00:00:00 2001 From: lf Date: Mon, 1 Jun 2026 11:31:35 +0800 Subject: [PATCH 26/28] refactor: provide homotopy default inner via NonlinearSolve extension MERGE-3: replace the construction-time `Base.require()` in `_default_inner` with a package extension. A weak dependency on NonlinearSolve plus the new `MTKNonlinearSolveExt` populates a factory Ref at load time with `NonlinearSolve.NewtonRaphson`; when NonlinearSolve is absent `_default_inner` falls back to `SimpleNonlinearSolve.SimpleNewtonRaphson` (a runtime dependency), so a `homotopy(...)` system stays buildable out of the box. Removes the eager `Base.require` and its hardcoded UUID. Test S3b asserts the factory Ref is populated and resolves to NewtonRaphson when NonlinearSolve is loaded. --- lib/ModelingToolkitBase/Project.toml | 3 +- .../ext/MTKNonlinearSolveExt.jl | 18 +++++++++++ .../src/systems/homotopy_sweep.jl | 31 ++++++++----------- .../test/homotopy_sweep.jl | 11 +++++++ 4 files changed, 44 insertions(+), 19 deletions(-) create mode 100644 lib/ModelingToolkitBase/ext/MTKNonlinearSolveExt.jl diff --git a/lib/ModelingToolkitBase/Project.toml b/lib/ModelingToolkitBase/Project.toml index a33f98b39e..3986fad0e5 100644 --- a/lib/ModelingToolkitBase/Project.toml +++ b/lib/ModelingToolkitBase/Project.toml @@ -72,6 +72,7 @@ JuliaFormatter = "98e50ef6-434e-11e9-1051-2b60c6c9e899" LabelledArrays = "2ee39098-c373-598a-b85f-a56591580800" Latexify = "23fbe1c1-3f47-55db-b15f-69d7ec21a316" Mooncake = "da2b9cff-9c12-43a0-ae48-6db2b0edb7d6" +NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec" Pyomo = "0e8e1daf-01b5-4eba-a626-3897743a3816" Tracker = "9f7883ad-71c0-57eb-9f7f-b5c9e6d3789c" @@ -86,6 +87,7 @@ MTKJuliaFormatterExt = "JuliaFormatter" MTKLabelledArraysExt = "LabelledArrays" MTKLatexifyExt = "Latexify" MTKMooncakeExt = "Mooncake" +MTKNonlinearSolveExt = "NonlinearSolve" MTKPyomoDynamicOptExt = "Pyomo" MTKTrackerExt = "Tracker" @@ -193,7 +195,6 @@ Latexify = "23fbe1c1-3f47-55db-b15f-69d7ec21a316" LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae" Logging = "56ddb016-857b-54e1-b83d-db4d58db5568" ModelingToolkitStandardLibrary = "16a59e39-deab-5bd0-87e4-056b12336739" -NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec" Optimization = "7f7a1694-90dd-40f0-9382-eb1efda571ba" OptimizationBase = "bca83a33-5cc9-4baa-983d-23429ab6bcbb" OptimizationIpopt = "43fad042-7963-4b32-ab19-e2a4f9a67124" diff --git a/lib/ModelingToolkitBase/ext/MTKNonlinearSolveExt.jl b/lib/ModelingToolkitBase/ext/MTKNonlinearSolveExt.jl new file mode 100644 index 0000000000..4579f4fc4f --- /dev/null +++ b/lib/ModelingToolkitBase/ext/MTKNonlinearSolveExt.jl @@ -0,0 +1,18 @@ +module MTKNonlinearSolveExt + +using ModelingToolkitBase: ModelingToolkitBase +using NonlinearSolve: NewtonRaphson + +# Make `NonlinearSolve.NewtonRaphson` the default inner solver for the homotopy +# continuation algorithms (`HomotopySweep`/`TrivialHomotopy`/`TrivialThenSweep`) +# whenever NonlinearSolve is loaded. NonlinearSolve is not a runtime dependency +# of ModelingToolkitBase, so `_default_inner` cannot reference it directly; this +# extension populates the factory Ref at load time. Without NonlinearSolve loaded +# ModelingToolkitBase falls back to `SimpleNonlinearSolve.SimpleNewtonRaphson` +# (a hard `[deps]`), so a `homotopy(...)` system is still buildable out of the box. +function __init__() + ModelingToolkitBase._DEFAULT_INNER_FACTORY[] = () -> NewtonRaphson() + return nothing +end + +end diff --git a/lib/ModelingToolkitBase/src/systems/homotopy_sweep.jl b/lib/ModelingToolkitBase/src/systems/homotopy_sweep.jl index cbd73c69c2..71a3acf346 100644 --- a/lib/ModelingToolkitBase/src/systems/homotopy_sweep.jl +++ b/lib/ModelingToolkitBase/src/systems/homotopy_sweep.jl @@ -43,25 +43,20 @@ function HomotopySweep(; inner = nothing, schedule = 0.0:0.1:1.0, schedule, set_λ!, maxiters_per_step) end -# Lazy-load NewtonRaphson from NonlinearSolve, which is a test-extras dependency -# of ModelingToolkitBase (not a runtime [deps]). Callers can avoid this path by -# passing `inner = ...` explicitly; the default is only invoked when the test -# env (or a downstream user env) has already loaded NonlinearSolve. +# Default inner solver for the homotopy continuation algorithms. NonlinearSolve +# is NOT a runtime `[deps]` of ModelingToolkitBase, so its richer `NewtonRaphson` +# cannot be referenced directly here. The `MTKNonlinearSolveExt` package +# extension (a `[weakdeps]` extension) populates this Ref with a thunk returning +# `NonlinearSolve.NewtonRaphson()` whenever NonlinearSolve is loaded. Until then +# `_default_inner` falls back to `SimpleNonlinearSolve.SimpleNewtonRaphson` (a +# runtime `[deps]`) so a system carrying `homotopy(...)` can be built and solved +# out of the box. Load NonlinearSolve (`using NonlinearSolve`) to restore its +# `NewtonRaphson` as the default, or pass `inner = ...` explicitly. +const _DEFAULT_INNER_FACTORY = Ref{Union{Nothing, Function}}(nothing) + function _default_inner() - try - mod = Base.require(Base.PkgId(Base.UUID("8913a72c-1f9b-4ce2-8d82-65094dcecaec"), - "NonlinearSolve")) - return mod.NewtonRaphson() - catch - # NonlinearSolve is a test-only dependency of ModelingToolkitBase, not a - # runtime [deps]. When it isn't available, fall back to SimpleNonlinearSolve - # (a runtime [deps]) so a system carrying `homotopy(...)` can be built and - # solved out of the box — without this, merely constructing a homotopy - # `ODEProblem` would throw at problem-construction time. Load NonlinearSolve - # (`using NonlinearSolve`) to restore its richer `NewtonRaphson` as the - # default, or pass `inner = ...` explicitly. - return SimpleNewtonRaphson() - end + factory = _DEFAULT_INNER_FACTORY[] + return factory === nothing ? SimpleNewtonRaphson() : factory() end function CommonSolve.solve(prob::NonlinearProblem, alg::HomotopySweep; kwargs...) diff --git a/lib/ModelingToolkitBase/test/homotopy_sweep.jl b/lib/ModelingToolkitBase/test/homotopy_sweep.jl index 4e12a64657..493878c5eb 100644 --- a/lib/ModelingToolkitBase/test/homotopy_sweep.jl +++ b/lib/ModelingToolkitBase/test/homotopy_sweep.jl @@ -51,6 +51,17 @@ using NonlinearSolve: NewtonRaphson @test alg.inner isa typeof(NewtonRaphson()) end + @testset "S3b default inner supplied by NonlinearSolve extension" begin + # MERGE-3: the default inner solver is provided by the + # `MTKNonlinearSolveExt` package extension (populated when NonlinearSolve + # is loaded), not by an eager construction-time `Base.require`. With + # NonlinearSolve loaded the factory Ref is set and resolves to + # `NonlinearSolve.NewtonRaphson`; without it, `_default_inner` falls back + # to `SimpleNonlinearSolve.SimpleNewtonRaphson` (a hard dependency). + @test ModelingToolkitBase._DEFAULT_INNER_FACTORY[] !== nothing + @test ModelingToolkitBase._default_inner() isa typeof(NewtonRaphson()) + end + @testset "S4 TrivialHomotopy dispatches to inner once" begin f! = (du, u, p) -> du[1] = u[1]^2 - 4 prob = NonlinearProblem(f!, [1.5], Float64[]) From a599a48dc4260c14710b7717436ebcfbe0432e42 Mon Sep 17 00:00:00 2001 From: lf Date: Mon, 1 Jun 2026 11:36:57 +0800 Subject: [PATCH 27/28] =?UTF-8?q?fix:=20reset=20homotopy=20=CE=BB=20to=20a?= =?UTF-8?q?ctual=20form=20on=20failed=20sweep=20step?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit SWEEP-2: when a `HomotopySweep` step returns an unsuccessful retcode the algorithm returns early, but the returned solution's parameter vector still held the intermediate λ where continuation stalled, exposing a half-homotopied system to any consumer reading those parameters. Reset λ to 1.0 (actual form) on the early-return path before surfacing the failed solution; the unknowns are left untouched so they still carry the diagnostic failed iterate. Test S2b drives a fixture whose real root vanishes around λ ≈ 0.5 and asserts the failed solution reports λ = 1.0. Test L1-Q7 is a regression guard that the observed equations of a lowered homotopy System are homotopy-free (an eliminated variable's definition living in observed must be lowered, not left as an opaque homotopy node). --- .../src/systems/homotopy_sweep.jl | 10 +++++- .../test/homotopy_lowering.jl | 36 +++++++++++++++++++ .../test/homotopy_sweep.jl | 22 ++++++++++++ 3 files changed, 67 insertions(+), 1 deletion(-) diff --git a/lib/ModelingToolkitBase/src/systems/homotopy_sweep.jl b/lib/ModelingToolkitBase/src/systems/homotopy_sweep.jl index 71a3acf346..afc98d12c8 100644 --- a/lib/ModelingToolkitBase/src/systems/homotopy_sweep.jl +++ b/lib/ModelingToolkitBase/src/systems/homotopy_sweep.jl @@ -77,7 +77,15 @@ function CommonSolve.solve(prob::NonlinearProblem, alg::HomotopySweep; kwargs... step_sol = CommonSolve.solve(step_prob, alg.inner; inner_kwargs...) last_sol = step_sol if !successful_retcode(step_sol) - return step_sol + # SWEEP-2: reset λ to actual form (1.0) before surfacing the failed + # solution. The sweep stalled at this step's intermediate λ; leaving + # λ there would expose a half-homotopied system to any consumer that + # reads the returned solution's parameter vector. The unknowns are + # left untouched (they carry the diagnostic failed iterate). + p_reset = copy(parameter_values(step_prob)) + ret = alg.set_λ!(p_reset, one(λ)) + new_p_reset = ret === nothing ? p_reset : ret + return @set step_sol.prob = remake(step_prob; p = new_p_reset) end u_curr = copy(step_sol.u) end diff --git a/lib/ModelingToolkitBase/test/homotopy_lowering.jl b/lib/ModelingToolkitBase/test/homotopy_lowering.jl index 967f420013..383797f9c4 100644 --- a/lib/ModelingToolkitBase/test/homotopy_lowering.jl +++ b/lib/ModelingToolkitBase/test/homotopy_lowering.jl @@ -141,4 +141,40 @@ using Symbolics prob_nohomotopy = ODEProblem(sys2, Dict(a => 1.0, q => 2.0), (0.0, 1.0)) @test !haskey(prob_nohomotopy.kwargs, :initializealg) end + + @testset "L1-Q7 observed equations are homotopy-free after lowering" begin + # Regression guard: `add_homotopy_parameter` lowers homotopy nodes in + # observed equations too (PressureDrop-style: an eliminated variable's + # definition lives in observed). Downstream observed codegen must never + # see an opaque `homotopy(...)` node — assert observed is homotopy-free + # after lowering a System whose homotopy-defined variable is eliminated. + using ModelingToolkitBase: System, mtkcompile, observed, has_homotopy + using ModelingToolkitBase: t_nounits as t, D_nounits as D + using OrdinaryDiffEqRosenbrock + using OrdinaryDiffEqNonlinearSolve + + # `w ~ homotopy(...)` defines w explicitly from a state, so mtkcompile + # eliminates w into observed — the codepath where an eliminated + # variable's definition (PressureDrop's m_flow) lands in observed. + @variables x(t) w(t) + @parameters p + eqs = [D(x) ~ -x, + w ~ homotopy(x^2 - p, x - 1)] + @named sys = System(eqs, t; guesses = [w => 2.5]) + sys = mtkcompile(sys) + + obs = observed(sys) + @test obs !== nothing && !isempty(obs) + for eq in obs + @test !has_homotopy(eq.lhs) + @test !has_homotopy(eq.rhs) + end + + # The lowered ODEProblem's observed must likewise be homotopy-free. + prob = ODEProblem(sys, Dict(x => 1.0, p => 9.0), (0.0, 1.0)) + for eq in observed(prob.f.sys) + @test !has_homotopy(eq.lhs) + @test !has_homotopy(eq.rhs) + end + end end diff --git a/lib/ModelingToolkitBase/test/homotopy_sweep.jl b/lib/ModelingToolkitBase/test/homotopy_sweep.jl index 493878c5eb..44de8a36a3 100644 --- a/lib/ModelingToolkitBase/test/homotopy_sweep.jl +++ b/lib/ModelingToolkitBase/test/homotopy_sweep.jl @@ -43,6 +43,28 @@ using NonlinearSolve: NewtonRaphson @test !SciMLBase.successful_retcode(sol) end + @testset "S2b failed step resets λ to actual form (1.0) before returning" begin + # SWEEP-2: when a sweep step fails mid-schedule, the returned solution's + # parameter vector must hold λ = 1.0 (actual form), not the intermediate + # λ where continuation stalled. Otherwise a downstream consumer reading + # the failed solution's parameters sees a half-homotopied system. + # Fixture: f = (1-λ)*u + λ*(u^2+1). A real root exists only for small λ; + # around λ ≈ 0.5 the quadratic loses its real root and Newton diverges, + # so the failing step's λ is strictly between 0 and 1. + f_bad! = (du, u, p) -> begin + λ = p[1] + du[1] = (1 - λ) * u[1] + λ * (u[1]^2 + 1) + end + prob = NonlinearProblem(f_bad!, [0.5], [0.0]) + alg = HomotopySweep(; inner = NewtonRaphson(), + schedule = 0.0:0.1:1.0, + set_λ! = set_λ_explicit, + maxiters_per_step = 20) + sol = solve(prob, alg) + @test !SciMLBase.successful_retcode(sol) + @test sol.prob.p[1] == 1.0 + end + @testset "S3 HomotopySweep default constructor" begin alg = HomotopySweep(; set_λ! = ((p, v) -> p)) @test alg.schedule == 0.0:0.1:1.0 From e1c85d23ad5ba7fc3a6969d50b6bfde14c6a899f Mon Sep 17 00:00:00 2001 From: lf Date: Tue, 2 Jun 2026 09:35:30 +0800 Subject: [PATCH 28/28] docs: update homotopy operator docstring to reflect shipped sweep MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The `homotopy` operator docstring still said the parameter-sweep continuation was deferred to "future PRs" — but this branch ships it. Since the docstring renders on the documented API page via `@docs ModelingToolkit.homotopy`, the stale wording would advertise the feature as not-yet-implemented. Rewrite it to describe the actual behavior: runtime evaluates to `actual`; during init the node is lowered to `(1-λ)*simplified + λ*actual` and solved trivial-then-sweep. --- lib/ModelingToolkitBase/src/systems/homotopy.jl | 16 +++++++++++++--- 1 file changed, 13 insertions(+), 3 deletions(-) diff --git a/lib/ModelingToolkitBase/src/systems/homotopy.jl b/lib/ModelingToolkitBase/src/systems/homotopy.jl index 888cae2ba5..d3688ca16a 100644 --- a/lib/ModelingToolkitBase/src/systems/homotopy.jl +++ b/lib/ModelingToolkitBase/src/systems/homotopy.jl @@ -3,9 +3,19 @@ """ homotopy(actual, simplified) -Modelica homotopy operator (spec 3.7.4). At runtime and at L0 trivial rewrite -time, equivalent to `actual`. Future PRs may use `simplified` as a starting -point for parameter-sweep continuation during initialization. +Modelica homotopy operator (spec 3.7.4), an initialization aid. At runtime, +outside initialization, it evaluates to `actual`. + +During `complete`/`mtkcompile` every `homotopy(actual, simplified)` node is +lowered to `(1 - λ)*simplified + λ*actual` with a single shared `__homotopy_λ` +parameter (default `1.0`, so `λ=1` reduces to `actual`). A system containing +`homotopy(...)` nodes is then initialized either by solving `actual` directly +(the spec's *trivial form*) or, when that fails, by continuation that sweeps `λ` +from 0 (the easy-to-solve `simplified` equation) to 1 (`actual`) — the spec's +*transformation form*. The default initialization algorithm is +`OverrideInit(nlsolve = TrivialThenSweep(...))`: trivial first, sweep on failure. + +See [`TrivialThenSweep`](@ref), [`HomotopySweep`](@ref) and [`TrivialHomotopy`](@ref). """ homotopy(actual::Real, simplified::Real) = actual