Homotopy

ModelingToolkit implements the Modelica homotopy(actual, simplified) operator (Modelica Specification 3.7.4.2) 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:

homotopy(actual, simplified)

Concretely, the continuation introduces a scalar parameter $\lambda$ and solves

\[(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

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.

Opting Out of the Lowering

Some targets cannot lower to a continuation solver — for example a GPU ensemble kernel with a fixed set of solvers. For those, pass homotopy = false to mtkcompile:

sys = mtkcompile(sys; homotopy = false)

Every homotopy(actual, simplified) node is then replaced by actual before compilation. The compiled system contains no homotopy nodes, so the generated code evaluates actual directly (the simplified expression is never emitted), AbstractNonlinearProblem(sys, op) and the initialization problem are plain nonlinear problems rather than a HomotopyProblem, and the initialization and event affect systems derived from the compiled system are compiled the same way. Only the simplified starting heuristic is discarded; the solved system is unchanged. ModelingToolkitBase.homotopy_enabled(sys) reports the setting a compiled system was built with.

Building a `HomotopyProblem`

A system whose equations contain homotopy nodes is built into a SciMLBase.HomotopyProblem 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:

# 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 can be solved with any algorithm that supports the problem type; solve(prob) with no algorithm picks a suitable default. The current default is natural-parameter continuation (NonlinearSolveBase.HomotopySweep), which sweeps λ from 0 to 1, solving a standard nonlinear problem at each step and warm-starting from the previous step's solution.

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:

using ModelingToolkit, NonlinearSolve

@variables y
@mtkcompile sys = System([0 ~ homotopy(atan(y - 3), y)])
prob = HomotopyProblem(sys, [y => 12.0])
sol = solve(prob)
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:

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:

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

ModelingToolkitBase.homotopy — Function

homotopy(actual, simplified)

The Modelica homotopy operator (Modelica Specification 3.7.4.2). Annotating an expression as homotopy(actual, simplified) declares that simplified is an easy-to-solve approximation of actual: a continuation (homotopy) solver can start from the simplified equations and continuously deform them into the actual ones. This stabilises the solution of nonlinear systems that are hard to solve from a cold start — pressure-driven flow networks, power-flow equations, chemical equilibria, equations with multiple roots — by walking from the easy solution to the true one. Wherever the operator is evaluated numerically (outside a continuation solve) it is simply actual, as the Modelica spec prescribes; initialization is one common consumer, not the only one.

Arguments

actual: the real expression. Used wherever the operator is evaluated numerically, and reached at the end of the continuation (λ = 1).
simplified: an approximation of actual whose solution is easy to reach from the available guess. Only influences the continuation path.

Both arguments are scalar Real expressions, matching Modelica's restriction of homotopy to scalar Reals; the runtime method is homotopy(actual::Real, simplified::Real) = actual. For array equations, broadcast the operator elementwise — homotopy.(actual, simplified) — which creates one homotopy node per element; the continuation lowering rewrites each node independently, and all of them share the single continuation parameter λ.

Behavior

homotopy(actual, simplified) stays an opaque symbolic operator through System construction, mtkcompile, and runtime code generation — no continuation parameter is injected into the system, and systems that do not use the operator go through a byte-identical pipeline. Targets that cannot lower to a continuation solver can opt out entirely with mtkcompile(sys; homotopy = false), which replaces every node by its actual branch before compilation (see strip_homotopy). 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, and along the continuation they follow the blended expression below.

Building a SciMLBase.HomotopyProblem from a system that contains homotopy nodes — directly with HomotopyProblem(sys, op), or automatically through AbstractNonlinearProblem(sys, op) — regenerates the residual with every homotopy(actual, simplified) replaced by the convex blend

(1 - λ) * simplified + λ * actual

compiled as f(u, p, λ) with λ an explicit trailing argument — λ is never added to the system's parameters, and the user's parameter object p passes through untouched. All homotopy calls in the lowered system share that single λ. Note that initialization applies this per strongly connected component: an init system that tears into SCC blocks builds a HomotopyProblem only for those blocks whose equations carry homotopy nodes, so each such block sweeps its own λ (in dependency order, reaching λ = 1 before the next block begins) while the remaining blocks keep their plain Newton solves. An init system with no SCC decomposition to hang that off is swept as a whole with one λ. The resulting HomotopyProblem (with λspan defaulting to (0.0, 1.0)) can be solved with any algorithm that supports it; solve(prob) with no algorithm picks a default that sweeps λ from 0 (simplified) to 1 (actual).

Example

The equation 0 = atan(y - 3) has its root at y = 3, but Newton from y = 12 diverges because atan saturates. Annotating with simplified = y (root at y = 0) lets the continuation walk to the true root:

using ModelingToolkit, NonlinearSolve

@variables y
@mtkcompile sys = System([0 ~ homotopy(atan(y - 3), y)])
prob = HomotopyProblem(sys, [y => 12.0])
sol = solve(prob)
sol[y] # ≈ 3.0 — the continuation rescued the out-of-basin guess

Notes

Runtime cost: outside a continuation solve the generated code calls the numeric fallback, so the simplified argument expression is evaluated and its value discarded (arguments are evaluated before the call). This small overhead is borne only by systems that use homotopy; all other systems are unaffected.
See the Homotopy documentation page for construction details and the continuation solver's options.

Reference: Modelica Specification 3.7.4.2 https://specification.modelica.org/master/operators-and-expressions.html#homotopy

source

ModelingToolkitBase.strip_homotopy — Function

strip_homotopy(sys::System)

Return a copy of sys in which every Modelica homotopy(actual, simplified) node has been replaced by its actual branch, recursively through subsystems. This is what mtkcompile(sys; homotopy = false) applies before compilation: the returned system contains no homotopy nodes, so generated code evaluates actual directly (the simplified expression is neither emitted nor evaluated) and problem construction never selects a SciMLBase.HomotopyProblem.

The rewrite covers the equations, observed equations, initialization equations, noise equations, costs, constraints and the values of bindings, initial conditions and guesses. Symbolic events and jumps are left as-is; a homotopy node inside them evaluates as actual through the numeric fallback.

source

ModelingToolkitBase.homotopy_enabled — Function

homotopy_enabled(sys::AbstractSystem)

Whether homotopy(actual, simplified) nodes in sys are lowered to a continuation solve. Returns false iff sys was compiled with mtkcompile(sys; homotopy = false) (see HomotopyCtx).

source

ModelingToolkitBase.HomotopyCtx — Type

HomotopyCtx

System metadata key recording the homotopy keyword argument that mtkcompile was called with. mtkcompile(sys; homotopy = false) strips every homotopy(actual, simplified) node down to actual and stores false under this key so that systems derived from the compiled one (the initialization system, event affect systems) are compiled the same way. Query it with homotopy_enabled.

source

ModelingToolkitBase.NonPolynomialReason — Module

Possible reasons why a term is not polynomial

source