Add the homotopy operator and HomotopyProblem construction

docs/pages.jl

@@ -47,6 +47,7 @@ pages = [
     "Basics" => Any[
         "basics/Composition.md",
         "basics/Events.md",
+        "basics/Homotopy.md",
         "basics/Linearization.md",
         "basics/InputOutput.md",
         "basics/MTKLanguage.md",

docs/src/API/problems.md

@@ -36,6 +36,8 @@ SciMLBase.ImplicitDiscreteProblem
 ```@docs
 SciMLBase.NonlinearFunction
 SciMLBase.NonlinearProblem
+SciMLBase.AbstractNonlinearProblem(::System, ::Any)
+SciMLBase.HomotopyProblem
 SciMLBase.SCCNonlinearProblem
 SciMLBase.NonlinearLeastSquaresProblem
 SciMLBase.SteadyStateProblem

docs/src/basics/Homotopy.md

@@ -0,0 +1,159 @@
+# [Homotopy](@id homotopy)
+
+ModelingToolkit implements the Modelica `homotopy(actual, simplified)` operator
+([Modelica Specification 3.7.4.2](https://specification.modelica.org/master/operators-and-expressions.html#homotopy))
+as a way to robustly solve nonlinear systems that are hard to solve from a cold
+start. It is most commonly reached for during initialization, but it is a general
+nonlinear-solving construct: any nonlinear system whose equations carry a
+`homotopy` annotation can be solved by continuation.
+
+## What Is Homotopy?
+
+The `simplified` equations are a set of equations for which the nonlinear system
+is easier to get a convergent (Newton) iteration for, and the `actual` equations
+are the more complex equations you actually want to solve. A homotopy solver
+starts by solving the nonlinear system with the `simplified` equations, and uses
+that solution as the starting point for solving the `actual` problem, deforming
+continuously from one to the other.
+
+There are many ways this can be used. For example, if you have equations with
+multiple solutions — like a quadratic equation with a positive and a negative
+root — you can simplify it down to an approximating linear problem that has a
+single (say, positive) solution, and then deform it to the actual equation to
+stabilise the process of converging to that positive solution.
+
+The operator encodes both expressions in a single annotation:
+
+```julia
+homotopy(actual, simplified)
+```
+
+Concretely, the continuation introduces a scalar parameter ``\lambda`` and solves
+
+```math
+(1 - \lambda)\,\text{simplified} + \lambda\,\text{actual}
+```
+
+sweeping ``\lambda`` from `0` (the easy `simplified` system) to `1` (the
+`actual` system), warm-starting each step from the previous solution.
+
+## Runtime Semantics
+
+The operator stays an opaque symbolic function through `System` construction,
+`mtkcompile`, and runtime code generation — no continuation parameter is added to
+the system. Wherever the operator is evaluated numerically outside a continuation
+solve, the generated code calls the numeric fallback
+
+```julia
+homotopy(actual::Real, simplified::Real) = actual
+```
+
+so the operator evaluates to `actual`, as the Modelica specification prescribes.
+Note the honest cost: the `simplified` argument expression is still evaluated and
+its value discarded (arguments are evaluated before the call). This small
+overhead is borne only by systems that use `homotopy` — systems without the
+operator go through a byte-identical pipeline and are completely unaffected.
+
+Symbolic differentiation works through the operator: nodewise derivative rules
+keep symbolic jacobians, `tgrad`, and index reduction consistent. At runtime,
+differentiated equations reproduce `actual`'s derivative; along the continuation
+they follow the derivative of the blended expression above.
+
+## Building a `HomotopyProblem`
+
+A system whose equations contain `homotopy` nodes is built into a
+[`SciMLBase.HomotopyProblem`](@ref) whose residual is the blended expression
+above, compiled as `f(u, p, λ)`. `λ` is an explicit trailing argument — it is
+never added to the system's parameters, and your parameter object `p` passes
+through untouched. All `homotopy` calls in a system share the single `λ`, per the
+Modelica spec's recommendation of (conceptually) one homotopy iteration over the
+whole model; this includes nested `homotopy` calls.
+
+There are two ways to construct it:
+
+```julia
+# Explicit: always returns a HomotopyProblem (errors if `sys` has no `homotopy`).
+prob = HomotopyProblem(sys, op)
+
+# Automatic: returns a HomotopyProblem when `sys` contains `homotopy` nodes, and
+# a plain NonlinearProblem otherwise.
+prob = AbstractNonlinearProblem(sys, op)
+```
+
+The `HomotopyProblem`'s `λspan` defaults to `(0.0, 1.0)`. It is solved by
+`NonlinearSolveBase.HomotopySweep`, a natural-parameter continuation solver that
+sweeps `λ` from `0` to `1`, solving a standard nonlinear problem at each step and
+warm-starting from the previous step's solution. A `HomotopyProblem` with no
+algorithm (`solve(prob)`) defaults to this solver.
+
+## Example: Out-of-Basin Rescue
+
+The equation `0 = atan(y - 3)` has a root at `y = 3`, but a Newton solver
+starting from `y = 12` diverges because `atan` saturates. Using `homotopy` with
+`simplified = y` (whose root is `y = 0`) lets the continuation walk from the easy
+root to the true one:
+
+```julia
+using ModelingToolkit, NonlinearSolve
+
+@variables y
+@mtkcompile sys = System([0 ~ homotopy(atan(y - 3), y)])
+prob = HomotopyProblem(sys, [y => 12.0])
+sol = solve(prob, HomotopySweep())
+sol[y] # ≈ 3.0 — the continuation rescued the out-of-basin guess
+```
+
+The operating point (`[y => 12.0]`) provides the starting point of the
+continuation; the sweep deforms the equations so the solver reaches `y ≈ 3` at
+`λ = 1`.
+
+## Broadcasting Over Arrays
+
+`homotopy` is a scalar operator. For array equations, broadcast it elementwise:
+
+```julia
+eqs = 0 .~ homotopy.(actual_array, simplified_array)
+```
+
+This creates one `homotopy` node per element; the continuation lowering rewrites
+each node independently, and all of them share the single continuation parameter
+`λ`.
+
+## Customizing the Continuation Solver
+
+To tune the sweep, pass your own continuation algorithm to `solve`:
+
+```julia
+sol = solve(prob, HomotopySweep(nsteps = 30))
+```
+
+`HomotopySweep` accepts the keyword arguments `inner` (the nonlinear algorithm
+used at each step), `nsteps`, `adaptive`, `initial_step_factor`, and `min_dλ`;
+see the `NonlinearSolveBase.HomotopySweep` docstring for their meanings and
+defaults.
+
+## Limitations
+
+  - **`expression = Val{true}` is not yet supported** for the homotopy
+    constructor; build the problem directly (the default
+    `expression = Val{false}`). This can be added in a future PR.
+  - **The jacobian/sparsity of the standard build are dropped.** They encode the
+    `λ = 1` (opaque-`actual`) system and would be wrong mid-sweep; continuation
+    steps solve with a freshly differentiated residual. Per-problem analytic
+    jacobians for the swept residual are future work.
+  - **Scalar `Real` expressions only** for a single `homotopy` node, matching
+    Modelica's restriction; use broadcasting (above) for arrays.
+  - **Only equations and observed equations are rewritten.** A `homotopy` call
+    inside a parameter binding or default value is left as-is and evaluates as
+    `actual` at all `λ`.
+
+## API Reference
+
+```@docs
+ModelingToolkit.homotopy
+```
+
+## See Also
+
+  - [Modelica Specification 3.7.4.2](https://specification.modelica.org/master/operators-and-expressions.html#homotopy)
+    — the upstream specification this operator implements.

lib/ModelingToolkitBase/Project.toml

@@ -143,7 +143,11 @@ ModelingToolkitStandardLibrary = "2.20"
 Mooncake = "0.5"
 Moshi = "0.3.6"
 NaNMath = "0.3, 1"
-NonlinearSolve = "4.3"
+NonlinearSolve = "4.19"
+# Test-only floor: `HomotopySweep` (which the homotopy tests solve with, via
+# NonlinearSolve's re-export) lives in NonlinearSolveBase >= 2.31; NonlinearSolve
+# 4.19's own bound (>= 2.20) does not force it, so pin it for downgrade CI.
+NonlinearSolveBase = "2.31"
 OffsetArrays = "1"
 OptimizationBase = "5"
 OptimizationIpopt = "1"
@@ -169,12 +173,12 @@ ReferenceTests = "0.10"
 RuntimeGeneratedFunctions = "0.5.12"
 SCCNonlinearSolve = "1.8.1"
 SafeTestsets = "0.1"
-SciMLBase = "3.18"
+SciMLBase = "3.19"
 SciMLPublic = "1.0.0"
 SciMLStructures = "1.7"
 Serialization = "1"
 Setfield = "0.7, 0.8, 1"
-SimpleNonlinearSolve = "0.1.0, 1, 2"
+SimpleNonlinearSolve = "2.11.1"
 SparseArrays = "1"
 SpecialFunctions = "1, 2"
 StableRNGs = "1"
@@ -208,6 +212,7 @@ LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
 Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
 ModelingToolkitStandardLibrary = "16a59e39-deab-5bd0-87e4-056b12336739"
 NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
+NonlinearSolveBase = "be0214bd-f91f-a760-ac4e-3421ce2b2da0"
 Optimization = "7f7a1694-90dd-40f0-9382-eb1efda571ba"
 OptimizationBase = "bca83a33-5cc9-4baa-983d-23429ab6bcbb"
 OptimizationIpopt = "43fad042-7963-4b32-ab19-e2a4f9a67124"
@@ -232,4 +237,4 @@ Sundials = "c3572dad-4567-51f8-b174-8c6c989267f4"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["BenchmarkTools", "BoundaryValueDiffEqMIRK", "BoundaryValueDiffEqAscher", "ControlSystemsBase", "DataInterpolations", "DelayDiffEq", "NonlinearSolve", "ForwardDiff", "ModelingToolkitStandardLibrary", "Optimization", "OptimizationIpopt", "OptimizationOptimJL", "OrdinaryDiffEq", "OrdinaryDiffEqDefault", "REPL", "Random", "ReferenceTests", "SafeTestsets", "StableRNGs", "Statistics", "SteadyStateDiffEq", "Test", "StochasticDiffEq", "Sundials", "Pkg", "OrdinaryDiffEqNonlinearSolve", "Logging", "OptimizationBase", "LinearSolve", "Latexify", "Distributed", "DiffEqNoiseProcess", "DynamicQuantities", "OrdinaryDiffEqSDIRK", "OrdinaryDiffEqRosenbrock", "OrdinaryDiffEqBDF", "OrdinaryDiffEqFunctionMap"]
+test = ["BenchmarkTools", "BoundaryValueDiffEqMIRK", "BoundaryValueDiffEqAscher", "ControlSystemsBase", "DataInterpolations", "DelayDiffEq", "NonlinearSolve", "ForwardDiff", "ModelingToolkitStandardLibrary", "NonlinearSolveBase", "Optimization", "OptimizationIpopt", "OptimizationOptimJL", "OrdinaryDiffEq", "OrdinaryDiffEqDefault", "REPL", "Random", "ReferenceTests", "SafeTestsets", "StableRNGs", "Statistics", "SteadyStateDiffEq", "Test", "StochasticDiffEq", "Sundials", "Pkg", "OrdinaryDiffEqNonlinearSolve", "Logging", "OptimizationBase", "LinearSolve", "Latexify", "Distributed", "DiffEqNoiseProcess", "DynamicQuantities", "OrdinaryDiffEqSDIRK", "OrdinaryDiffEqRosenbrock", "OrdinaryDiffEqBDF", "OrdinaryDiffEqFunctionMap"]

lib/ModelingToolkitBase/src/ModelingToolkitBase.jl

@@ -53,7 +53,7 @@ using Compat
 using AbstractTrees
 using SciMLBase: StandardODEProblem, StandardNonlinearProblem, handle_varmap, TimeDomain,
     PeriodicClock, Clock, SolverStepClock, ContinuousClock, OverrideInit,
-    NoInit
+    NoInit, AbstractNonlinearProblem
 import Moshi
 using Moshi.Data: @data
 using Reexport
@@ -191,6 +191,9 @@ include("systems/codegen_utils.jl")
 include("problems/docs.jl")
 include("systems/codegen.jl")
 include("systems/problem_utils.jl")
+# Operator + lowering layer; must load before the problem constructors that
+# consume it (problems/nonlinearproblem.jl selector + problems/homotopyproblem.jl).
+include("systems/homotopy_operator.jl")
 
 include("problems/compatibility.jl")
 include("problems/odeproblem.jl")
@@ -199,6 +202,7 @@ include("problems/daeproblem.jl")
 include("problems/sdeproblem.jl")
 include("problems/sddeproblem.jl")
 include("problems/nonlinearproblem.jl")
+include("problems/homotopyproblem.jl")
 include("problems/intervalnonlinearproblem.jl")
 include("problems/implicitdiscreteproblem.jl")
 include("problems/discreteproblem.jl")
@@ -246,6 +250,10 @@ export ImplicitDiscreteProblem, ImplicitDiscreteFunction
 export ODEProblem, SDEProblem
 export NonlinearFunction
 export NonlinearProblem
+# `AbstractNonlinearProblem(sys, op)` is the automatic constructor that selects a
+# `HomotopyProblem` when `sys` contains `homotopy(...)` nodes (see
+# `problems/homotopyproblem.jl`); re-exported so it is reachable unqualified.
+export AbstractNonlinearProblem
 export IntervalNonlinearFunction
 export IntervalNonlinearProblem
 export OptimizationProblem, constraints
@@ -300,6 +308,8 @@ export generate_initializesystem, Initial, isinitial, InitializationProblem
 export alg_equations, diff_equations, has_alg_equations, has_diff_equations
 export get_alg_eqs, get_diff_eqs, has_alg_eqs, has_diff_eqs
 
+export homotopy
+
 export @variables, @parameters, @independent_variables, @constants, @brownians, @brownian,
     @poissonians, @discretes
 export @named, @nonamespace, @namespace, extend, compose, complete, toggle_namespacing

lib/ModelingToolkitBase/src/problems/homotopyproblem.jl

@@ -0,0 +1,83 @@
+"""
+    SciMLBase.HomotopyProblem(sys::System, op; λspan = (0.0, 1.0), kwargs...)
+
+Build a [`SciMLBase.HomotopyProblem`](@ref) from a `System` whose equations
+contain Modelica `homotopy(actual, simplified)` nodes (Modelica spec 3.7.4.2).
+
+The standard nonlinear residual is built from `sys` (and used only to obtain
+`u0`, `p`, the observed function, and the residual prototype), then the residual
+is regenerated with every `homotopy(actual, simplified)` replaced by the convex
+blend `(1 - λ)*simplified + λ*actual` and compiled as `f(u, p, λ)`. `λ` is an
+explicit trailing argument — it is never added to the system's parameters, and
+`p` passes through untouched — so the result solves by natural-parameter
+continuation (`NonlinearSolveBase.HomotopySweep`), sweeping `λ` across `λspan`
+(default `(0.0, 1.0)`, i.e. from `simplified` to `actual`).
+
+`sys` must be `complete`d and must contain at least one `homotopy` node (use
+[`SciMLBase.NonlinearProblem`](@ref) otherwise, or `AbstractNonlinearProblem` to
+select automatically). The jacobian/sparsity of the standard build are
+deliberately dropped: they encode the `λ = 1` (opaque-`actual`) system and would
+be wrong mid-sweep.
+
+!!! note
+
+    `expression = Val{true}` (codegen-to-`Expr`) is not yet supported for the
+    homotopy constructor; this can be added in a future PR.
+"""
+function SciMLBase.HomotopyProblem(
+        sys::System, op;
+        expression = Val{false}, λspan = (0.0, 1.0),
+        check_length = true, check_compatibility = true,
+        eval_expression = false, eval_module = @__MODULE__,
+        checkbounds = false, cse = true, kwargs...
+    )
+    if expression !== Val{false}
+        throw(
+            ArgumentError(
+                "`HomotopyProblem(sys, op)` does not yet support " *
+                    "`expression = Val{true}`; build the problem directly " *
+                    "(the default `expression = Val{false}`)."
+            )
+        )
+    end
+    check_complete(sys, SciMLBase.HomotopyProblem)
+    if is_time_dependent(sys)
+        sys = NonlinearSystem(sys)
+    end
+    check_compatibility && check_compatible_system(SciMLBase.NonlinearProblem, sys)
+    if !has_any_homotopy(sys)
+        throw(
+            ArgumentError(
+                "`HomotopyProblem` requires a system containing " *
+                    "`homotopy(actual, simplified)` nodes; none were found. Use " *
+                    "`NonlinearProblem` for systems without `homotopy`."
+            )
+        )
+    end
+
+    _iip = resolve_iip(Both, op)
+    f, u0,
+        p = process_SciMLProblem(
+        SciMLBase.NonlinearFunction{_iip}, sys, op;
+        check_length, check_compatibility, expression,
+        eval_expression, eval_module, checkbounds, cse, kwargs...
+    )
+
+    # Swap the opaque-`actual` residual for the homotopy-swept `f(u, p, λ)`. The
+    # observed function and residual prototype carry over; `initialization_data`,
+    # jacobian, and sparsity are deliberately not carried (the latter two encode
+    # the `λ = 1` system and would be wrong mid-sweep).
+    shadow, λ = lower_homotopy(sys)
+    hf = generate_homotopy_residual(
+        shadow, λ; eval_expression, eval_module, checkbounds, cse
+    )
+    swept_f = SciMLBase.NonlinearFunction{_iip}(
+        hf; sys = f.sys, observed = f.observed, resid_prototype = f.resid_prototype
+    )
+
+    kwargs = process_kwargs(sys; kwargs...)
+    args = (; f = swept_f, u0, p)
+    return maybe_codegen_scimlproblem(
+        expression, SciMLBase.HomotopyProblem{_iip}, args; λspan, kwargs...
+    )
+end

lib/ModelingToolkitBase/src/problems/nonlinearproblem.jl

@@ -118,6 +118,31 @@ end
     )
 end
 
+"""
+    get_nonlinear_problem_type(sys::System)
+
+The concrete `AbstractNonlinearProblem` type that `AbstractNonlinearProblem(sys, op)`
+builds for `sys`: [`SciMLBase.HomotopyProblem`](@ref) when `sys` contains Modelica
+`homotopy(actual, simplified)` nodes, otherwise [`SciMLBase.NonlinearProblem`](@ref).
+"""
+function get_nonlinear_problem_type(sys::System)
+    return has_any_homotopy(sys) ? SciMLBase.HomotopyProblem : SciMLBase.NonlinearProblem
+end
+
+"""
+    SciMLBase.AbstractNonlinearProblem(sys::System, op; kwargs...)
+
+Build a nonlinear problem from `sys`, automatically selecting the concrete type via
+`get_nonlinear_problem_type`: a [`SciMLBase.HomotopyProblem`](@ref) when `sys`
+contains Modelica `homotopy(actual, simplified)` nodes (so the equations are solved by
+continuation from the `simplified` form), otherwise a plain
+[`SciMLBase.NonlinearProblem`](@ref). Keyword arguments are forwarded to the selected
+constructor; `λspan` only applies to the `HomotopyProblem` branch.
+"""
+function SciMLBase.AbstractNonlinearProblem(sys::System, op; kwargs...)
+    return get_nonlinear_problem_type(sys)(sys, op; kwargs...)
+end
+
 @fallback_iip_specialize function SciMLBase.NonlinearLeastSquaresProblem{iip, spec}(
         sys::System, op; check_length = false, lb = nothing, ub = nothing,
         check_compatibility = true, expression = Val{false}, kwargs...