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 +``` 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/ModelingToolkitBase.jl b/lib/ModelingToolkitBase/src/ModelingToolkitBase.jl index da845565bf..cd8429a57a 100644 --- a/lib/ModelingToolkitBase/src/ModelingToolkitBase.jl +++ b/lib/ModelingToolkitBase/src/ModelingToolkitBase.jl @@ -194,6 +194,8 @@ include("systems/ir_info.jl") 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") @@ -270,6 +272,7 @@ export liouville_transform, change_independent_variable, export respecialize export PDESystem export Differential, expand_derivatives, @derivatives +export homotopy, HomotopySweep, TrivialHomotopy, TrivialThenSweep export Equation export Term export SymScope, LocalScope, ParentScope, GlobalScope 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/initializationproblem.jl b/lib/ModelingToolkitBase/src/problems/initializationproblem.jl index e24346cfea..19a03194b5 100644 --- a/lib/ModelingToolkitBase/src/problems/initializationproblem.jl +++ b/lib/ModelingToolkitBase/src/problems/initializationproblem.jl @@ -82,6 +82,40 @@ 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(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_λ) + 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 + end + if simplify_system isys = mtkcompile(isys; fully_determined, split = is_split(sys), initsys_mtkcompile_kwargs...) end 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 42d041208e..fa97c78216 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()) @@ -139,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/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 new file mode 100644 index 0000000000..d3688ca16a --- /dev/null +++ b/lib/ModelingToolkitBase/src/systems/homotopy.jl @@ -0,0 +1,213 @@ +@register_symbolic homotopy(actual, simplified) + +""" + homotopy(actual, simplified) + +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 + +""" + rewrite_trivial(expr) + +Recursively replace every `homotopy(a, s)` node in `expr` with `a`. Preserves +`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 = unwrap(expr) + return _rewrite_trivial(x) +end + +function _rewrite_trivial(x) + if !iscall(x) + return x + end + op = operation(x) + args = arguments(x) + new_args = map(_rewrite_trivial, args) + if op === homotopy + return new_args[1] + end + return maketerm(typeof(x), op, new_args, metadata(x)) +end + +""" + has_homotopy(expr) + +Return `true` iff `expr` contains at least one `homotopy(...)` node. +""" +function has_homotopy(expr) + x = unwrap(expr) + return _has_homotopy(x) +end + +function _has_homotopy(x) + if !iscall(x) + return false + end + if operation(x) === homotopy + return true + 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 + +""" + 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 + +""" + 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 diff --git a/lib/ModelingToolkitBase/src/systems/homotopy_sweep.jl b/lib/ModelingToolkitBase/src/systems/homotopy_sweep.jl new file mode 100644 index 0000000000..afc98d12c8 --- /dev/null +++ b/lib/ModelingToolkitBase/src/systems/homotopy_sweep.jl @@ -0,0 +1,149 @@ +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` +# 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) + +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 + +# 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() + factory = _DEFAULT_INNER_FACTORY[] + return factory === nothing ? SimpleNewtonRaphson() : factory() +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, λ) + 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 ? + 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) + # 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 + 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; _inner_solve_kwargs(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/src/systems/problem_utils.jl b/lib/ModelingToolkitBase/src/systems/problem_utils.jl index 6e1036c028..4cfae663bd 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)) @@ -1953,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 = (;) @@ -1976,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_init.jl b/lib/ModelingToolkitBase/test/homotopy_init.jl new file mode 100644 index 0000000000..3f1bb048cb --- /dev/null +++ b/lib/ModelingToolkitBase/test/homotopy_init.jl @@ -0,0 +1,159 @@ +using Test +using ModelingToolkitBase +using ModelingToolkitBase: rewrite_trivial, has_homotopy +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 !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 + + @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 + + @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 + + @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(t) y(t) + @parameters 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, 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 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 diff --git a/lib/ModelingToolkitBase/test/homotopy_integration.jl b/lib/ModelingToolkitBase/test/homotopy_integration.jl new file mode 100644 index 0000000000..5499548787 --- /dev/null +++ b/lib/ModelingToolkitBase/test/homotopy_integration.jl @@ -0,0 +1,147 @@ +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 + + @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 + + @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 + # 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 diff --git a/lib/ModelingToolkitBase/test/homotopy_lowering.jl b/lib/ModelingToolkitBase/test/homotopy_lowering.jl new file mode 100644 index 0000000000..383797f9c4 --- /dev/null +++ b/lib/ModelingToolkitBase/test/homotopy_lowering.jl @@ -0,0 +1,180 @@ +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 + + @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 + + @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 + + @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_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/homotopy_sweep.jl b/lib/ModelingToolkitBase/test/homotopy_sweep.jl new file mode 100644 index 0000000000..44de8a36a3 --- /dev/null +++ b/lib/ModelingToolkitBase/test/homotopy_sweep.jl @@ -0,0 +1,135 @@ +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 "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 + # 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 "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[]) + 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`. + # 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(), + 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 e024cc0200..0d01f30dc9 100644 --- a/lib/ModelingToolkitBase/test/runtests.jl +++ b/lib/ModelingToolkitBase/test/runtests.jl @@ -63,6 +63,11 @@ 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 "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