Modelica homotopy(actual, simplified) operator for initialization (spec 3.7.4)#4601
Modelica homotopy(actual, simplified) operator for initialization (spec 3.7.4)#4601BatchZero wants to merge 34 commits into
homotopy(actual, simplified) operator for initialization (spec 3.7.4)#4601Conversation
|
Status note: per the review on SciML/SciMLBase.jl#1375 and SciML/NonlinearSolve.jl#962, the homotopy parameter λ is being separated from the parameter vector (the residual becomes |
Ports the Modelica homotopy(actual, simplified) operator registration
and the rewrite_with_lambda symbolic-lowering core from the phase1a
reference branch into current master (feature/homotopy-lowering).
- @register_symbolic homotopy(actual, simplified) + Real runtime fallback
- rewrite_with_lambda / _rewrite_with_lambda: replaces homotopy nodes with
(1-λ)*simplified + λ*actual using a shared __homotopy_λ parameter
- rewrite_with_lambda_in_equations: batch version for equation vectors
- has_homotopy / has_homotopy_in_equations / has_any_homotopy: predicate
helpers needed by the upcoming add_homotopy_parameter (Task 2)
- include("systems/homotopy.jl") wired into ModelingToolkitBase.jl
- test/homotopy_lowering.jl: 3 assertions covering λ=1 (actual), λ=0
(simplified), and __homotopy_λ name; all pass green
…t param buffer resolves it
…naffected + easy-guess init
Add four @safetestset entries to lib/ModelingToolkitBase/test/runtests.jl under the `Initialization` group so CI picks them up once the SciMLBase + NonlinearSolve deps (HomotopyProblem / HomotopySweep) land in the registry.
…ext on HomotopySweep, NEWS entry - Rewrite homotopy(actual,simplified) docstring: remove nonexistent TrivialThenSweep / TrivialHomotopy mentions and their broken @refs; describe the actual behavior (HomotopySweep continuation, spec §3.7.4). - MTKNonlinearSolveExt: guard factory assignment with isdefined(NonlinearSolve, :HomotopySweep) so the ext loads cleanly against pre-release NonlinearSolve versions without HomotopySweep. Import NonlinearSolve as a module instead of importing the symbol directly. - Add Unreleased section to NEWS.md with homotopy operator entry.
…dle Val{true} codegen path
A4: forward `initprob.kwargs...` when wrapping NonlinearProblem into
HomotopyProblem so any kwargs the init problem carried (callbacks, etc.)
are not silently dropped.
A3: drop three dead-code idioms —
- `observed(sys)` always returns `Vector{Equation}`, never `nothing`;
remove the `obs === nothing` guard in `has_any_homotopy` and the
`obs !== nothing &&` guard in `add_homotopy_parameter`.
- Belt-and-suspenders `unique!` on the λ parameter-injection in
`add_homotopy_parameter` to match the `add_initialization_parameters`
idiom (guard already prevents the common double-inject).
- Replace the `getfield_safe` indirection in
`maybe_inject_homotopy_initializealg` with a direct `f.initialization_data`
access; `initialization_data` is a concrete field on every SciML*Function
type, so the `hasproperty` dance was always-true defense. Remove the
now-dead `getfield_safe` helper entirely.
A1: document (not fix) the known limitation that the homotopy
`initializealg` auto-injection is only wired for the `Val{false}`
(default) problem-construction path. The `Val{true}` path builds an
`Expr` AST via `build_scimlproblem_expr`; `args.f.initialization_data`
cannot be inspected at expression-building time, so injecting the
`initializealg` there would require generating conditional AST nodes
referencing runtime state — non-trivial and risky. Users of
`expression = Val{true}` with a homotopy system must pass
`initializealg` explicitly.
…ise derivatives Per the merged SciMLBase#1375 / NonlinearSolve#962 f(u,p,λ) convention, no __homotopy_λ parameter is injected anywhere; user parameters pass through untouched. Differentiation keeps the two-branch structure for OMC-parity under index reduction.
…rings Review fixes on top of the opaque-homotopy rework: - add a testset differentiating through the SECOND argument, pinning both the runtime endpoint (compiled derivative of actual w.r.t. b is 0) and the blend side via _rewrite_with_lambda ((1-λ)*2b = 2.8 at λ=0.3, b=2), so the ∂₂ = homotopy(0,1) registration is no longer untested - correct the rationale comment: SConst wrapping constructs the node from symbolic constants explicitly; it is not mandated for term() arguments - reword docstrings (homotopy, has_any_homotopy, _rewrite_with_lambda) to describe current behavior only, dropping references to init-pipeline machinery that is not wired at this commit - drop the isdefined(:add_homotopy_parameter) assertion and the unused SymbolicUtils import from the test file Note for bisectability: the sibling homotopy test files (homotopy_integration.jl, homotopy_init.jl, homotopy_omc_parity.jl) are knowingly stale mid-rework at this commit and are rewritten in subsequent commits.
… coverage, nested lowering via public API; intent-level docstring
…opy_for_init; test tidy-up
…et; warn on unsweepable
…etection; test/comment accuracy fixes
… drop extension HomotopySweep lives in NonlinearSolveBase 2.31, which is already in MTKBase's unconditional load closure via SimpleNonlinearSolve; gating the default behind a NonlinearSolve weakdep added load-order state for nothing — a hard dep adds zero packages to the resolved set, deletes the type-unstable factory Ref, and the 2.31 compat floor replaces the isdefined guard.
… floor 2 SimpleNonlinearSolve >= 2 is what makes the load-closure claim (MTKBase -> SimpleNonlinearSolve -> NonlinearSolveBase) readable from the Project.toml alone — 0.1/1.x don't depend on NonlinearSolveBase and were only excluded indirectly via the SciMLBase floor. The new testset-2 assertion passes initializealg at ODEProblem construction and checks nsteps == 30 survives (the injected default has nsteps === nothing), pinning the haskey gate in maybe_inject_homotopy_initializealg rather than just DiffEqBase's solve-kwarg precedence.
…eep, remake, aliasing
…eep failure contract
Runic v1.7.0 (matches FormatCheck's runic-action version spec "1") applied to every .jl file in the branch diff. Formatting-only; no semantic changes. Hunks outside the branch's own regions in codegen.jl, problem_utils.jl and both runtests.jl are pre-existing format debt verified Runic-dirty on upstream master as well.
NonlinearSolveBase 2.31 co-resolves only with SimpleNonlinearSolve 2.9.0+ (within the 2.x series), and SciMLBase 3.19 — the floor this branch already declares — only admits SimpleNonlinearSolve 2.11.1+. With the floor at "2", julia-downgrade-compat pins 2.0.0 and the DowngradeSublibraries resolve fails. 2.11.1 is the verified minimal co-resolvable floor (checked together with SciMLBase 3.19.0, NonlinearSolveBase 2.31.0, NonlinearSolve 4.19.0).
BatchZero
force-pushed
the
feature/homotopy-lowering
branch
from
6f2690c to
4706e06
Compare
|
Both upstream dependencies are merged and registered: SciML/SciMLBase.jl#1375 → SciMLBase v3.19.0 ( |
| @@ -0,0 +1,167 @@ | |||
| # [Homotopy Initialization](@id homotopy) | |||
There was a problem hiding this comment.
Why even speak to initialization here? It's just a nonlinear solver. The fact that it shows up commonly in initialization is a detail, but there are many other uses for it.
| Problems built through the default path then receive | ||
|
|
||
| ```julia | ||
| initializealg = SciMLBase.OverrideInit(nlsolve = NonlinearSolveBase.HomotopySweep()) |
There was a problem hiding this comment.
Any nonlinear solve of equations which have a homotopy operator should default to defining a homotopy problem.
There was a problem hiding this comment.
Or rather, there should be an AbstractNonlinearProblem constructor that is automatic here.
| - If the initialization would otherwise target an `SCCNonlinearProblem`, a | ||
| homotopy-using system is forced to the monolithic `NonlinearProblem` target, | ||
| since the continuation deforms the coupled system along one `λ`. |
There was a problem hiding this comment.
It should probably split that to SCCs.
| falls back to least squares), a one-time warning is emitted and initialization | ||
| solves the `actual` (`λ = 1`) form instead. Least-squares and per-SCC sweeps | ||
| are future work. | ||
| - **Trivial initialization is silently correct.** If the initialization system |
There was a problem hiding this comment.
Obviously. Why is this even worth mentioning?
| - **`expression = Val{true}` does not auto-inject the `initializealg`.** When | ||
| generating problem-construction expressions, pass | ||
| `initializealg = SciMLBase.OverrideInit(nlsolve = NonlinearSolveBase.HomotopySweep())` | ||
| explicitly. |
| - **Scalar `Real` expressions only**, matching Modelica's restriction; the | ||
| numeric fallback method is `homotopy(::Real, ::Real)`. |
| - **Scalar `Real` expressions only**, matching Modelica's restriction; the | ||
| numeric fallback method is `homotopy(::Real, ::Real)`. | ||
| - **Failure surfaces as a failed initialization.** If the continuation cannot |
There was a problem hiding this comment.
Obviously. Also, not relevant to the MTK side.
| homotopic = has_any_homotopy(isys) | ||
| if homotopic && TProb === SCCNonlinearProblem | ||
| # A continuation sweep deforms the COUPLED system along one λ; per-SCC | ||
| # sweeping is a separate future improvement. Root ModelingToolkit's | ||
| # `get_initialization_problem_type` prefers `SCCNonlinearProblem` for split | ||
| # fully-determined systems — force the monolithic target here. | ||
| TProb = NonlinearProblem | ||
| end | ||
| initprob = TProb{_iip}(isys, op; kwargs..., build_initializeprob = false, is_initializeprob = true) | ||
| if homotopic | ||
| initprob = _maybe_sweep_homotopy_initprob(initprob, isys, _iip, kwargs) | ||
| end | ||
| return initprob |
There was a problem hiding this comment.
Split to the AbstractNonlinearProblem constructor.
| rhss = [_iszero(eq.lhs) ? eq.rhs : eq.rhs - eq.lhs for eq in eqs] | ||
| end | ||
| else | ||
| # NOTE: generate_homotopy_residual in homotopy.jl mirrors the time-independent path here; keep in sync. |
| ) | ||
| ) | ||
| end | ||
| λname = gensym(:__homotopy_λ) |
There was a problem hiding this comment.
gensym will break precompilation. Instead use a hard to hit sentinal with unicodes.
| @@ -1,3 +1,22 @@ | |||
| # Unreleased | |||
There was a problem hiding this comment.
? it'll release right away? Why write that here?
|
Do one PR first with just the HomotopyProblem constructor and direct interface, and the AbstractNonlinearProblem stuff. Then do a follow up with initialization. |
| homotopy(actual, simplified) | ||
| ``` | ||
|
|
||
| During initialization the solver may start from the easy `simplified` form and | ||
| continuously deform it into `actual`; outside initialization the operator | ||
| evaluates to `actual`, as the Modelica specification prescribes. |
There was a problem hiding this comment.
It should be more clear. The simplified equations are a set of equations for which the set of nonlinear equations is easier to get a convergent (Newton) iteration for, and the actual equations are the more complex equations you want to actually solve for. Homotopy solvers start by solving the nonlinear system with the simplified equations, and uses that solution as a starting point to solve the actual problem. There are many ways this can be used. For example if you have equations that have multiple solutions, like a quadratic equation with positive and negative solutions, you can simplify it down to an approximating linear problem which has one problem (the positive solution), and then deform it to the actual in order to stabilize the process of getting that positive solution.
There was a problem hiding this comment.
This also highlights that your implementation of homotopy has a bug since it doesn't start at lambda=0 😅
2 participants
Summary
Implements Modelica's
homotopy(actual, simplified)operator (Modelica Specification §3.7.4.2) for initialization: outside initialization the operator is semanticallyactual; during initialization the system may be solved by a homotopy deformingsimplified(λ = 0) intoactual(λ = 1). This is the standard Modelica mechanism for robustly initializing pressure-driven flow networks and other hard nonlinear init systems (the Modelica Buildings library relies on it heavily).MTK's part is an independent lowering pass targeting an MTK-independent API, per the review verdict on #4576. The solver stack merged and registered upstream first:
SciMLBase.HomotopyProblem(AddHomotopyProblem: a homotopy/continuation problem type SciMLBase.jl#1375, registered in SciMLBase v3.19.0) — residual conventionf(u, p, λ)out-of-place /f(du, u, p, λ)in-place; λ is a trailing scalar argument, never an entry ofp;NonlinearSolveBase.HomotopySweep(AddHomotopySweep: natural-parameter continuation solver forHomotopyProblemNonlinearSolve.jl#962, registered in NonlinearSolveBase v2.31.0, re-exported by NonlinearSolve) — the natural-parameter continuation solver dispatching onHomotopyProblem.This PR contains lowering + targeting only. A grep over the full
srcdiff finds zero sweep/λ-stepping logic (9 doc/comment mentions plus the singleHomotopySweep()construction at the injection site).What changed since the original push
As committed in the status note on this PR, the lowering was reworked once #1375/#962 settled the
f(u, p, λ)convention. The branch is rebased onto current master (d95960f), 34 commits. Concretely:__homotopy_λparameter atcomplete()and shipped ahomotopy_parameterlocator intop; per the fixcollect_difference_variables#1375/Fix length ofis_linear_equations#962 reviews, λ is now a trailing scalar argument of the generated residual. Userp/MTKParameterspass through untouched; the locator design is gone.complete()to the targeting boundary. The parentSystemis never mutated; a shadow rewrite regenerates only the initialization residual atInitializationProblembuild time (post-mtkcompile).isdefinedguard and type-unstable factoryRef), replaced by a hardNonlinearSolveBasedependency — see Dependencies & compat for why this is free.How it works
homotopy(actual, simplified)is registered and opaque throughSystemconstruction,mtkcompile, and runtime codegen. The numeric fallback isactual—simplifiedis evaluated and its value discarded, an honest cost borne only by annotated systems (eliding the dead branch is listed as an optional follow-up optimization). Systems that don't use the operator are byte-unaffected.homotopy(1, 0), ∂₂ =homotopy(0, 1). At runtime these fold to the derivative ofactual; under the swept rewrite they blend to λ and 1 − λ respectively — so equations that index reduction has differentiated keep the same blend structure as the undifferentiated ones, matching OMC's treatment.InitializationProblembuild, a detection scan runs over equations ∪ observed of the compiled init system. For the squareNonlinearProblem-with-unknowns case, the residual is regenerated from a shadow rewrite(1 − λ)·simplified + λ·actualwith one sharedgensymλ across all operator occurrences (a bare symbolic, never a parameter, so it cannot clash with or capture user-chosen names), compiled asf(u, p, λ)with λ in the t-slot of the generated function; an explicitp_endkeeps λ out of theMTKParametersdestructuring. The generator mirrorsgenerate_rhswith exactly two deltas and a keep-in-sync back-pointer comment incodegen.jl. The result is wrapped as aSciMLBase.HomotopyProblem;λspanis inherited from the upstreamHomotopyProblemdefault — this PR introduces no numeric defaults of its own.SCCNonlinearProblemfor split, fully-determined init systems. Homotopy systems are forced to the monolithicNonlinearProblemtarget (the sweep is global over λ; a per-SCC sweep is a follow-up). Validated by a root-MTK end-to-end test against an SCC-routed twin system, with a negative control confirming the non-homotopy SCC fast path is preserved.Activation rules
Pinned by tests, except the three marked otherwise:
HomotopyProblem+ defaultHomotopySweep).initializealgalways wins — gated at both the constructor andsolveprecedence.homotopyfolded away by simplification ⇒ standard init, no sweep machinery (no dedicated test; covered by the detection gate — after simplification there is no operator occurrence left to detect).homotopycalls ⇒ one shared λ, one sweep.actual, which is correct.@warn+ ordinary solve ofactual(λ = 1); an NLLS sweep is a follow-up.expression = Val{true}⇒ no auto-injection — documented limitation, no dedicated test (passinitializealgexplicitly; follow-up).homotopyin parameter bindings/defaults is not rewritten — it evaluates asactualat all λ (the operator targets the equation system being initialized); documented limitation, no dedicated test.Why the default injection exists
When (and only when)
initializeprob isa HomotopyProblemand the user passed noinitializealg, the standard (Val{false}) build path injectsOverrideInit(nlsolve = NonlinearSolveBase.HomotopySweep())into the problem kwargs. The reason injection is needed at all: OrdinaryDiffEqNonlinearSolve'sdefault_nlsolve(::Nothing, …, ::AbstractNonlinearProblem)substitutes a poly-algorithm that has noHomotopyProblemmethod, whereas an explicitOverrideInitalg passes throughdefault_nlsolveas the identity; tolerances are inherited from the integrator. The clean long-term hook is upstream and declared below as a follow-up PR —default_nlsolve(::Nothing, …, ::HomotopyProblem) = HomotopySweep()in OrdinaryDiffEqNonlinearSolve — after which this injection branch shrinks away.@code_warntypeon the injected-alg construction is concrete.Dependencies & compat
NonlinearSolveBasebecomes a hard dependency of ModelingToolkitBase, floor2.31(the registered release carryingHomotopySweep). It is already in MTKBase's unconditional load closure via the existing hard depSimpleNonlinearSolve, so the hard dep adds zero packages to the resolved set (verified withPkg.resolve). Precedent: the existing hardimport SCCNonlinearSolve. This is what allows deleting the weakdep extension, itsisdefinedguard, and the type-unstable factoryRef.SciMLBase = "3.19"(carriesHomotopyProblem),NonlinearSolveBase = "2.31", andSimpleNonlinearSolve = "2.11.1"— all three hard deps in[deps]; test depNonlinearSolve = "4.19". TheSimpleNonlinearSolvefloor moves from"0.1.0, 1, 2"to"2.11.1"because that is the lowest version co-resolvable with the declared SciMLBase ≥ 3.19 / NonlinearSolveBase ≥ 2.31 floors (DowngradeSublibraries pins[deps]floors; 2.11.1 verified co-resolvable, while older 2.x releases are not admitted by SciMLBase 3.19). The root NonlinearSolve is deliberately not floored at 4.20: 4.20.0 is unregistered, while the registered 4.19.1 re-exportsHomotopySweepand co-resolves with NonlinearSolveBase 2.31 — co-resolution of all four floors verified.Tests
All green locally:
Initializationgroup, full run — green; the only non-green entries are 12 pre-existingBrokenatinitializationsystem.jl:959(SimpleNewtonRaphson matrix entries), plus 2 version-gated ones at:1735-1736off Julia 1.12 — all pre-existing on upstream master. Zero new failures.lib/ModelingToolkitBase/test/homotopy_lowering.jl(31): operator opacity throughSystem/mtkcompile, ∂₁/∂₂ rules incl. mutation-resistant pins, the rewrite pass, gensym non-capture, nested operators.lib/ModelingToolkitBase/test/homotopy_integration.jl(32): generated-residual endpoints and blend (oop + iip), misuse guards, wrapper assertions, folded-away, trivial-init silence, non-square warn, end-to-end default rescue, structured parameters (incl. aninteg.ps[c]round-trip check), multi-homotopy single sweep,remake, u0 aliasing.lib/ModelingToolkitBase/test/homotopy_init.jl(7): injection chain, constructor-gate discriminator, solve precedence, absence of all of it on non-homotopy systems.lib/ModelingToolkitBase/test/homotopy_omc_parity.jl(16): OpenModelica v1.27 numerical parity on the Buildings PressureDrop pattern —m_flow(0) = √5to 1e-6 withm_floweliminated into observed (i.e. the swept-observed path), through both the default-injection and explicit-alg routes.test/initialization_homotopy.jl(6): SCC forcing under full MTK with an SCC-routed twin + negative control; non-homotopy SCC fast path preserved.Touched files are formatted with Runic 1.7.0 (the CI version), which picks up a few pre-existing format-debt hunks in those files that are Runic-dirty on master as well.
Documentation
New docs page
docs/src/basics/Homotopy.md(registered inpages.jl), a NEWS.md entry, and SciMLStyle docstrings on the exportedhomotopyoperator and the new internals.Maintenance surface
One feature-local file —
lib/ModelingToolkitBase/src/systems/homotopy.jl(operator + derivative rules + rewrite pass + residual codegen) — plus one named helper and a forcing line ininitializationproblem.jl, one gated injection branch inproblem_utils.jl, and a one-line keep-in-sync comment incodegen.jl. All continuation/solver maintenance lives in NonlinearSolve(Base).Limitations & follow-ups
Each of these can be a follow-up PR:
actualsolve).default_nlsolve(::Nothing, …, ::HomotopyProblem) = HomotopySweep()in OrdinaryDiffEqNonlinearSolve (retires the injection branch here).expression = Val{true}(currently a documented limitation).simplifiedbranch (optional optimization; the current cost is borne only by annotated systems).Relationship to #4576
Supersedes #4576, per the review there:
The solver now lives in NonlinearSolve(Base) and the problem type in SciMLBase; this PR is the targeting side only. #4576 will be closed in favor of this PR, as promised there.
Progresses #1385 (Modelica's homotopy operator), #3988 (homotopy initialization), and #4039 (cyclic guesses).