|
| 1 | +# [Homotopy](@id homotopy) |
| 2 | + |
| 3 | +ModelingToolkit implements the Modelica `homotopy(actual, simplified)` operator |
| 4 | +([Modelica Specification 3.7.4.2](https://specification.modelica.org/master/operators-and-expressions.html#homotopy)) |
| 5 | +as a way to robustly solve nonlinear systems that are hard to solve from a cold |
| 6 | +start. It is most commonly reached for during initialization, but it is a general |
| 7 | +nonlinear-solving construct: any nonlinear system whose equations carry a |
| 8 | +`homotopy` annotation can be solved by continuation. |
| 9 | + |
| 10 | +## What Is Homotopy? |
| 11 | + |
| 12 | +The `simplified` equations are a set of equations for which the nonlinear system |
| 13 | +is easier to get a convergent (Newton) iteration for, and the `actual` equations |
| 14 | +are the more complex equations you actually want to solve. A homotopy solver |
| 15 | +starts by solving the nonlinear system with the `simplified` equations, and uses |
| 16 | +that solution as the starting point for solving the `actual` problem, deforming |
| 17 | +continuously from one to the other. |
| 18 | + |
| 19 | +There are many ways this can be used. For example, if you have equations with |
| 20 | +multiple solutions — like a quadratic equation with a positive and a negative |
| 21 | +root — you can simplify it down to an approximating linear problem that has a |
| 22 | +single (say, positive) solution, and then deform it to the actual equation to |
| 23 | +stabilise the process of converging to that positive solution. |
| 24 | + |
| 25 | +The operator encodes both expressions in a single annotation: |
| 26 | + |
| 27 | +```julia |
| 28 | +homotopy(actual, simplified) |
| 29 | +``` |
| 30 | + |
| 31 | +Concretely, the continuation introduces a scalar parameter ``\lambda`` and solves |
| 32 | + |
| 33 | +```math |
| 34 | +(1 - \lambda)\,\text{simplified} + \lambda\,\text{actual} |
| 35 | +``` |
| 36 | + |
| 37 | +sweeping ``\lambda`` from `0` (the easy `simplified` system) to `1` (the |
| 38 | +`actual` system), warm-starting each step from the previous solution. |
| 39 | + |
| 40 | +## Runtime Semantics |
| 41 | + |
| 42 | +The operator stays an opaque symbolic function through `System` construction, |
| 43 | +`mtkcompile`, and runtime code generation — no continuation parameter is added to |
| 44 | +the system. Wherever the operator is evaluated numerically outside a continuation |
| 45 | +solve, the generated code calls the numeric fallback |
| 46 | + |
| 47 | +```julia |
| 48 | +homotopy(actual::Real, simplified::Real) = actual |
| 49 | +``` |
| 50 | + |
| 51 | +so the operator evaluates to `actual`, as the Modelica specification prescribes. |
| 52 | +Note the honest cost: the `simplified` argument expression is still evaluated and |
| 53 | +its value discarded (arguments are evaluated before the call). This small |
| 54 | +overhead is borne only by systems that use `homotopy` — systems without the |
| 55 | +operator go through a byte-identical pipeline and are completely unaffected. |
| 56 | + |
| 57 | +Symbolic differentiation works through the operator: nodewise derivative rules |
| 58 | +keep symbolic jacobians, `tgrad`, and index reduction consistent. At runtime, |
| 59 | +differentiated equations reproduce `actual`'s derivative; along the continuation |
| 60 | +they follow the derivative of the blended expression above. |
| 61 | + |
| 62 | +## Building a `HomotopyProblem` |
| 63 | + |
| 64 | +A system whose equations contain `homotopy` nodes is built into a |
| 65 | +[`SciMLBase.HomotopyProblem`](@ref) whose residual is the blended expression |
| 66 | +above, compiled as `f(u, p, λ)`. `λ` is an explicit trailing argument — it is |
| 67 | +never added to the system's parameters, and your parameter object `p` passes |
| 68 | +through untouched. All `homotopy` calls in a system share the single `λ`, per the |
| 69 | +Modelica spec's recommendation of (conceptually) one homotopy iteration over the |
| 70 | +whole model; this includes nested `homotopy` calls. |
| 71 | + |
| 72 | +There are two ways to construct it: |
| 73 | + |
| 74 | +```julia |
| 75 | +# Explicit: always returns a HomotopyProblem (errors if `sys` has no `homotopy`). |
| 76 | +prob = HomotopyProblem(sys, op) |
| 77 | + |
| 78 | +# Automatic: returns a HomotopyProblem when `sys` contains `homotopy` nodes, and |
| 79 | +# a plain NonlinearProblem otherwise. |
| 80 | +prob = AbstractNonlinearProblem(sys, op) |
| 81 | +``` |
| 82 | + |
| 83 | +The `HomotopyProblem`'s `λspan` defaults to `(0.0, 1.0)`. It can be solved with |
| 84 | +any algorithm that supports the problem type; `solve(prob)` with no algorithm |
| 85 | +picks a suitable default. The current default is natural-parameter continuation |
| 86 | +(`NonlinearSolveBase.HomotopySweep`), which sweeps `λ` from `0` to `1`, solving a |
| 87 | +standard nonlinear problem at each step and warm-starting from the previous |
| 88 | +step's solution. |
| 89 | + |
| 90 | +## Example: Out-of-Basin Rescue |
| 91 | + |
| 92 | +The equation `0 = atan(y - 3)` has a root at `y = 3`, but a Newton solver |
| 93 | +starting from `y = 12` diverges because `atan` saturates. Using `homotopy` with |
| 94 | +`simplified = y` (whose root is `y = 0`) lets the continuation walk from the easy |
| 95 | +root to the true one: |
| 96 | + |
| 97 | +```julia |
| 98 | +using ModelingToolkit, NonlinearSolve |
| 99 | + |
| 100 | +@variables y |
| 101 | +@mtkcompile sys = System([0 ~ homotopy(atan(y - 3), y)]) |
| 102 | +prob = HomotopyProblem(sys, [y => 12.0]) |
| 103 | +sol = solve(prob) |
| 104 | +sol[y] # ≈ 3.0 — the continuation rescued the out-of-basin guess |
| 105 | +``` |
| 106 | + |
| 107 | +The operating point (`[y => 12.0]`) provides the starting point of the |
| 108 | +continuation; the sweep deforms the equations so the solver reaches `y ≈ 3` at |
| 109 | +`λ = 1`. |
| 110 | + |
| 111 | +## Broadcasting Over Arrays |
| 112 | + |
| 113 | +`homotopy` is a scalar operator. For array equations, broadcast it elementwise: |
| 114 | + |
| 115 | +```julia |
| 116 | +eqs = 0 .~ homotopy.(actual_array, simplified_array) |
| 117 | +``` |
| 118 | + |
| 119 | +This creates one `homotopy` node per element; the continuation lowering rewrites |
| 120 | +each node independently, and all of them share the single continuation parameter |
| 121 | +`λ`. |
| 122 | + |
| 123 | +## Customizing the Continuation Solver |
| 124 | + |
| 125 | +To tune the sweep, pass your own continuation algorithm to `solve`: |
| 126 | + |
| 127 | +```julia |
| 128 | +sol = solve(prob, HomotopySweep(nsteps = 30)) |
| 129 | +``` |
| 130 | + |
| 131 | +`HomotopySweep` accepts the keyword arguments `inner` (the nonlinear algorithm |
| 132 | +used at each step), `nsteps`, `adaptive`, `initial_step_factor`, and `min_dλ`; |
| 133 | +see the `NonlinearSolveBase.HomotopySweep` docstring for their meanings and |
| 134 | +defaults. |
| 135 | + |
| 136 | +## Limitations |
| 137 | + |
| 138 | + - **`expression = Val{true}` is not yet supported** for the homotopy |
| 139 | + constructor; build the problem directly (the default |
| 140 | + `expression = Val{false}`). This can be added in a future PR. |
| 141 | + - **The jacobian/sparsity of the standard build are dropped.** They encode the |
| 142 | + `λ = 1` (opaque-`actual`) system and would be wrong mid-sweep; continuation |
| 143 | + steps solve with a freshly differentiated residual. Per-problem analytic |
| 144 | + jacobians for the swept residual are future work. |
| 145 | + - **Scalar `Real` expressions only** for a single `homotopy` node, matching |
| 146 | + Modelica's restriction; use broadcasting (above) for arrays. |
| 147 | + - **Only equations and observed equations are rewritten.** A `homotopy` call |
| 148 | + inside a parameter binding or default value is left as-is and evaluates as |
| 149 | + `actual` at all `λ`. |
| 150 | + |
| 151 | +## API Reference |
| 152 | + |
| 153 | +```@docs |
| 154 | +ModelingToolkit.homotopy |
| 155 | +``` |
| 156 | + |
| 157 | +## See Also |
| 158 | + |
| 159 | + - [Modelica Specification 3.7.4.2](https://specification.modelica.org/master/operators-and-expressions.html#homotopy) |
| 160 | + — the upstream specification this operator implements. |
0 commit comments