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

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

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",

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
+```

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"

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

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

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)

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

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)
 

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...)

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)

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)