diff --git a/docs/pages.jl b/docs/pages.jl index 4d01d67324..13ae9dfb0e 100644 --- a/docs/pages.jl +++ b/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", diff --git a/docs/src/API/problems.md b/docs/src/API/problems.md index d8cb7e940c..74f169ce73 100644 --- a/docs/src/API/problems.md +++ b/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 diff --git a/docs/src/basics/Homotopy.md b/docs/src/basics/Homotopy.md new file mode 100644 index 0000000000..4d04621cfb --- /dev/null +++ b/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. diff --git a/lib/ModelingToolkitBase/Project.toml b/lib/ModelingToolkitBase/Project.toml index 53b9915178..c4aeb1e124 100644 --- a/lib/ModelingToolkitBase/Project.toml +++ b/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"] diff --git a/lib/ModelingToolkitBase/src/ModelingToolkitBase.jl b/lib/ModelingToolkitBase/src/ModelingToolkitBase.jl index d4d8995d2c..f6b62fe4e6 100644 --- a/lib/ModelingToolkitBase/src/ModelingToolkitBase.jl +++ b/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 diff --git a/lib/ModelingToolkitBase/src/problems/homotopyproblem.jl b/lib/ModelingToolkitBase/src/problems/homotopyproblem.jl new file mode 100644 index 0000000000..f1f576af00 --- /dev/null +++ b/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 diff --git a/lib/ModelingToolkitBase/src/problems/nonlinearproblem.jl b/lib/ModelingToolkitBase/src/problems/nonlinearproblem.jl index 2eb3d0bf51..bec8f66dc1 100644 --- a/lib/ModelingToolkitBase/src/problems/nonlinearproblem.jl +++ b/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... diff --git a/lib/ModelingToolkitBase/src/systems/codegen.jl b/lib/ModelingToolkitBase/src/systems/codegen.jl index 260ab7d98e..25c98bfff2 100644 --- a/lib/ModelingToolkitBase/src/systems/codegen.jl +++ b/lib/ModelingToolkitBase/src/systems/codegen.jl @@ -30,6 +30,10 @@ $GENERATE_X_KWARGS `is_discrete_system(sys)`. - `scalar`: Whether to generate a single-out-of-place function that returns a scalar for the only equation in the system. +- `extra_args`: Extra trailing symbolic arguments appended after all standard arguments and + kept out of the `MTKParameters` collapse, so they stay live scalar arguments of the + generated function (e.g. the homotopy continuation `λ`, threaded as the "t slot"). Empty + by default, which leaves the standard codegen path byte-identical. All other keyword arguments are forwarded to [`build_function_wrapper`](@ref). """ @@ -37,7 +41,8 @@ function generate_rhs( sys::System; implicit_dae = false, scalar = false, expression = Val{true}, wrap_gfw = Val{false}, eval_expression = false, eval_module = @__MODULE__, override_discrete = false, - cachesyms = nothing, compiler_options::CompilerOptions = CompilerOptions(), + cachesyms = nothing, extra_args::Tuple = (), + compiler_options::CompilerOptions = CompilerOptions(), kwargs... ) dvs = unknowns(sys) @@ -107,10 +112,18 @@ function generate_rhs( args = (ddvs, args...) p_start += 1 end + # Extra trailing symbolic arguments (e.g. the homotopy continuation `λ`), + # appended after all standard arguments and kept out of the `MTKParameters` + # collapse via `p_end`. Empty by default, so every standard caller is + # byte-identical: `p_end` is only forwarded when `extra_args` is non-empty, + # and `nargs` already counts the extra arguments through `length(args)`. + args = (args..., extra_args...) + p_end_kw = isempty(extra_args) ? (;) : + (; p_end = (t === nothing ? length(args) : length(args) - 1) - length(extra_args)) u_arg = scalar ? -1 : (implicit_dae ? 2 : 1) res = build_function_wrapper( - sys, rhss, args...; p_start, extra_assignments, u_arg, + sys, rhss, args...; p_start, extra_assignments, u_arg, p_end_kw..., expression = Val{true}, expression_module = eval_module, kwargs... ) nargs = length(args) - length(p) + 1 diff --git a/lib/ModelingToolkitBase/src/systems/homotopy_operator.jl b/lib/ModelingToolkitBase/src/systems/homotopy_operator.jl new file mode 100644 index 0000000000..3cab9e3196 --- /dev/null +++ b/lib/ModelingToolkitBase/src/systems/homotopy_operator.jl @@ -0,0 +1,307 @@ +# Modelica's `homotopy(actual, simplified)` operator (spec 3.7.4.2) and its +# lowering to a continuation residual. NOT to be confused with +# `systems/nonlinear/homotopy_continuation.jl`, which is the unrelated polynomial +# HomotopyContinuation.jl rooting feature. +@register_symbolic homotopy(actual, simplified) + +""" + 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 `Real`s; 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. 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`](@ref) 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 + +```julia +(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 a system share the single `λ`, per the +Modelica spec's recommendation of (conceptually) one homotopy iteration over the +whole model. The resulting `HomotopyProblem` (with `λspan` defaulting to +`(0.0, 1.0)`) is solved by natural-parameter continuation +(`NonlinearSolveBase.HomotopySweep`), which 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: + +```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 +``` + +# 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](@ref homotopy) documentation page for construction details + and the continuation solver's options. + +**Reference:** Modelica Specification 3.7.4.2 + +""" +homotopy(actual::Real, simplified::Real) = actual + +# Nodewise symbolic derivatives (Modelica-faithful). +# +# ∂homotopy(a, s)/∂a = homotopy(1, 0) and ∂homotopy(a, s)/∂s = homotopy(0, 1), +# as OPAQUE nodes — deliberately NOT folded to the constants 1/0. The chain +# rule then yields `homotopy(1, 0)*a′ + homotopy(0, 1)*s′`, so that: +# * at runtime the numeric fallback folds the factors to 1/0, reproducing +# exactly the derivative of `actual` (Modelica spec: homotopy ≡ actual +# outside a continuation solve); +# * the continuation rewrite `homotopy(a, s) → (1-λ)s + λa` blends the same +# factors to λ/(1-λ), i.e. the true derivative of the blended expression — +# this keeps symbolic jacobians and index reduction (dummy derivatives) +# consistent with OMC's homotopy lowering. +# `@register_derivative` requires the rule to return a symbolic result; +# `SymbolicUtils.term` builds the opaque call node without invoking the numeric +# fallback on the constant arguments, and the constants are wrapped in +# `Symbolics.SConst` to construct the node explicitly from symbolic constants. +Symbolics.@register_derivative homotopy(actual, simplified) 1 SymbolicUtils.term( + homotopy, Symbolics.SConst(1), Symbolics.SConst(0); type = Real +) +Symbolics.@register_derivative homotopy(actual, simplified) 2 SymbolicUtils.term( + homotopy, Symbolics.SConst(0), Symbolics.SConst(1); type = Real +) + +""" + 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 + +""" + has_any_homotopy(sys) + +Return `true` iff `sys` contains a `homotopy(...)` node anywhere in its +equations OR its observed equations. +""" +function has_any_homotopy(sys) + has_homotopy_in_equations(equations(sys)) && return true + # `observed(sys)` always returns a `Vector{Equation}`, never `nothing`. + return any(eq -> has_homotopy(eq.lhs) || has_homotopy(eq.rhs), observed(sys)) +end + +""" + _rewrite_with_lambda(x, λ) + +Recursively replace every `homotopy(a, s)` node in the unwrapped expression `x` +with `(1 - λ)*s + λ*a`, where `λ` is supplied by the caller. At λ=1 the lowered +expression reduces numerically to `actual` (trivial form); at λ=0 it reduces to +`simplified`. + +Hand-written TermInterface walk rather than a SymbolicUtils `@rule` — keeps the +recursion explicit and avoids the rule-rewriter overhead for a single-node rewrite. +""" +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 + +# The continuation parameter is a fixed sentinel symbol, NOT `gensym`. A +# `gensym` name embeds a process-global counter, so the same system would lower +# to a different name (and thus a different generated `Expr`) in the precompile +# process than in the user session — defeating the `RuntimeGeneratedFunctions` +# Expr-hash cache and breaking precompilation. A fixed name makes codegen a pure +# function of the system. The `ₘₜₖ` unicode suffix (matching MTK's other +# internal sentinels, e.g. `__log_assertions_ₘₜₖ`) makes a collision with a +# user-chosen symbol practically impossible; `lower_homotopy` additionally +# guards against it, since `Sym` equality is name-keyed and a collision would +# otherwise silently rewrite `λ` into a state/parameter access and drop the +# `λ`-dependence of the residual. +const HOMOTOPY_LAMBDA = unwrap(only(@variables __homotopy_λₘₜₖ)) + +""" + lower_homotopy(sys::AbstractSystem) -> (shadow, λ) + +Independent lowering pass that feeds homotopy-swept residual code generation. +Returns a shadow copy of `sys` in which every `homotopy(actual, simplified)` +node in the equations and observed equations is replaced by the convex blend +`(1 - λ)*simplified + λ*actual`, plus the single shared bare symbolic `λ` +([`HOMOTOPY_LAMBDA`](@ref)). + +`sys` must be a `complete`d (flattened) system: the pass reads equations and +observed through accessors but reconstructs the toplevel fields, which is only +coherent once subsystems have been flattened away. + +`λ` is a fixed reserved sentinel and is **not** added to the system's parameters +— it is intended to become an explicit trailing argument of the swept residual +rather than an entry of `p`, per `SciMLBase.HomotopyProblem`'s `f(u, p, λ)` +contract. The shadow never passes through `mtkcompile` (`TearingState` would +reject the free λ); it exists only to feed the swept-residual code generation. +""" +function lower_homotopy(sys::AbstractSystem) + if !iscomplete(sys) + throw( + ArgumentError( + "`lower_homotopy` expects a `complete`d (flattened) system." + ) + ) + end + λ = HOMOTOPY_LAMBDA + λname = getname(λ) + # Guard against a user symbol named exactly like the sentinel: a collision in + # the unknowns/parameters would silently rewrite λ into a state/parameter + # access and drop the residual's λ-dependence. Observed lhs are included for + # completeness (a variable eliminated into observed by `mtkcompile`). + for v in Iterators.flatten( + ( + unknowns(sys), parameters(sys; initial_parameters = true), + (eq.lhs for eq in observed(sys)), + ) + ) + if hasname(v) && getname(v) === λname + throw( + ArgumentError( + "the homotopy continuation sentinel `$(λname)` collides with a " * + "symbol of the system; this name is reserved for internal use." + ) + ) + end + end + rew(x) = _rewrite_with_lambda(unwrap(x), λ) + new_eqs = [Equation(rew(eq.lhs), rew(eq.rhs)) for eq in equations(sys)] + new_obs = [Equation(rew(eq.lhs), rew(eq.rhs)) for eq in observed(sys)] + # All other fields (incl. `tearing_state`/`schedule`) are carried over by + # reference and remain valid: the rewrite is structure-preserving — same + # equation count/order/incidence, modulo the free λ. + return ConstructionBase.setproperties(sys; eqs = new_eqs, observed = new_obs), λ +end + +""" + generate_homotopy_residual(shadow, λ; kwargs...) + +Compile the swept residual `f(u, p, λ)` (out-of-place) / `f(du, u, p, λ)` +(in-place) from the homotopy-lowered shadow system. `λ` is the bare symbol +returned by [`lower_homotopy`](@ref); it is appended as an explicit trailing +scalar argument (the "t slot" of the 4-argument convention used by +`SciMLBase.HomotopyProblem`), and the user's parameter object is passed through +untouched (λ is excluded from the `MTKParameters` destructuring). + +This is exactly [`generate_rhs`](@ref) for the time-independent case with `λ` +threaded through its `extra_args`, behind two loud guards: both misuse modes +(a non-lowered system with raw `homotopy(...)` nodes, or a time-dependent +system) were probed to fail SILENTLY with wrong numerics rather than erroring. + +# Keyword Arguments + +$GENERATE_X_KWARGS + +All other keyword arguments are forwarded to [`generate_rhs`](@ref). +""" +function generate_homotopy_residual( + shadow::AbstractSystem, λ; expression = Val{false}, wrap_gfw = Val{true}, + kwargs... + ) + if has_any_homotopy(shadow) + throw( + ArgumentError( + "`generate_homotopy_residual` expects a homotopy-lowered shadow " * + "system as returned by `lower_homotopy`; this system still " * + "contains raw `homotopy(...)` nodes, which codegen would fold " * + "to `actual` via the numeric fallback, producing a " * + "λ-independent residual." + ) + ) + end + if is_time_dependent(shadow) + throw( + ArgumentError( + "`generate_homotopy_residual` expects a time-independent system; " * + "got a time-dependent one." + ) + ) + end + return generate_rhs( + shadow; extra_args = (λ,), expression, wrap_gfw, kwargs... + ) +end diff --git a/lib/ModelingToolkitBase/test/homotopy_lowering.jl b/lib/ModelingToolkitBase/test/homotopy_lowering.jl new file mode 100644 index 0000000000..338b2fb52e --- /dev/null +++ b/lib/ModelingToolkitBase/test/homotopy_lowering.jl @@ -0,0 +1,207 @@ +using ModelingToolkitBase +using ModelingToolkitBase: homotopy, has_homotopy, has_any_homotopy, + _rewrite_with_lambda, lower_homotopy, + t_nounits as t, D_nounits as D +using Symbolics +using Test + +@testset "homotopy stays opaque through complete(): no λ parameter is injected" begin + @variables x(t) y(t) + @named sys = System([D(x) ~ homotopy(-x, -x + 1), y ~ homotopy(x^2, x)], t) + csys = complete(sys) + @test !any(p -> occursin("homotopy", string(p)), parameters(csys)) + @test has_any_homotopy(csys) # nodes survive untouched +end + +@testset "runtime numeric fallback is the actual branch" begin + @test ModelingToolkitBase.homotopy(2.0, 5.0) == 2.0 +end + +@testset "registered derivative is nodewise (Modelica-faithful)" begin + @variables a b + # d/da homotopy(a^2, b) = homotopy(1,0) * 2a — the homotopy(1,0) factor is an + # opaque node (NOT folded to 1) so the continuation rewrite can blend it to λ. + d = Symbolics.derivative(homotopy(a^2, b), a) + @test has_homotopy(d) + # numeric endpoint: compiled evaluation folds homotopy(1,0) via the runtime + # fallback, so the derivative at a=3 must equal actual's derivative 6. + # (`substitute` does not re-fold a node whose arguments are already constant, + # so the genuine runtime path — codegen + numeric fallback — is exercised.) + fd = Symbolics.build_function(d, a, b; expression = Val(false)) + @test fd(3.0, 7.0) ≈ 6.0 +end + +@testset "registered derivative w.r.t. the simplified argument (∂₂)" begin + @variables a b λ + # d/db homotopy(a, b^2) = homotopy(0,1) * 2b — differentiation path through + # the SECOND argument exercises the ∂₂ = homotopy(0,1) registration. + d2 = Symbolics.derivative(homotopy(a, b^2), b) + @test has_homotopy(d2) + # runtime: actual = a does not depend on b, so the compiled derivative is 0. + fd2 = Symbolics.build_function(d2, a, b; expression = Val(false)) + @test fd2(3.0, 2.0) == 0.0 + # blend side: lowering the factor homotopy(0,1) → (1-λ)*1 + λ*0 = 1-λ yields + # the true derivative of the blended expression through the simplified + # branch: d/db [(1-λ)b^2 + λa] = (1-λ)*2b = (1-0.3)*2*2 = 2.8. + blended = _rewrite_with_lambda(Symbolics.unwrap(d2), Symbolics.unwrap(λ)) + fblend = Symbolics.build_function(blended, λ, b; expression = Val(false)) + @test fblend(0.3, 2.0) ≈ 2.8 +end + +@testset "symbolic jacobian of a homotopy system builds and matches actual" begin + @variables x(t) + @named sys = System([D(x) ~ homotopy(-x^2, -x)], t) + csys = complete(sys) + J = ModelingToolkitBase.calculate_jacobian(csys) + # at runtime homotopy ≡ actual: J = d(-x^2)/dx = -2x; check numerically at x=2 + fJ = Symbolics.build_function(only(J), unknowns(csys)[1]; expression = Val(false)) + @test fJ(2.0) ≈ -4.0 +end + +@testset "lower_homotopy: one shared λ across equations and observed" begin + @variables x(t) y(t) z(t) + @named sys = System( + [ + D(x) ~ homotopy(-x, -x + 1), + 0 ~ homotopy(atan(z - 3), z), + ], t; + observed = [y ~ homotopy(x^2, x)] + ) + # un-`complete`d systems are rejected (the pass reconstructs toplevel fields) + @test_throws ArgumentError lower_homotopy(sys) + csys = complete(sys) + shadow, λ = lower_homotopy(csys) + λn = Symbolics.wrap(λ) + # no homotopy nodes remain anywhere + @test !has_any_homotopy(shadow) + # single λ identity: substituting THE returned λ resolves every blend, at + # BOTH endpoints, in EVERY equation and in observed — a per-equation λ + # would leave at least one blend unresolved here. + eq1 = equations(shadow)[1] + @test isequal(Symbolics.substitute(eq1.rhs, Dict(λn => 1.0)), -x) + @test isequal(Symbolics.substitute(eq1.rhs, Dict(λn => 0.0)), -x + 1) + eq2 = equations(shadow)[2] + @test isequal(Symbolics.substitute(eq2.rhs, Dict(λn => 1.0)), atan(z - 3)) + @test isequal(Symbolics.substitute(eq2.rhs, Dict(λn => 0.0)), z) + obs = only(observed(shadow)) + @test isequal(Symbolics.substitute(obs.rhs, Dict(λn => 1.0)), x^2) + @test isequal(Symbolics.substitute(obs.rhs, Dict(λn => 0.0)), x) + # λ is NOT a parameter of the shadow + @test !any(p -> isequal(Symbolics.unwrap(p), λ), Symbolics.unwrap.(parameters(shadow))) +end + +@testset "fixed sentinel λ does not capture a near-name user symbol" begin + @variables x(t) + @parameters __homotopy_λ # close to, but NOT, the `__homotopy_λₘₜₖ` sentinel + @named sys = System([D(x) ~ homotopy(-x * __homotopy_λ, -x + 1)], t) + csys = complete(sys) + shadow, λ = lower_homotopy(csys) + λn = Symbolics.wrap(λ) + userλ = Symbolics.unwrap(__homotopy_λ) + # the sentinel λ is a distinct symbol (Sym equality is name-keyed) + @test !isequal(λ, userλ) + @test isequal(λ, ModelingToolkitBase.HOMOTOPY_LAMBDA) + # the user's parameter survives lowering; the sentinel λ is not a parameter + @test any( + p -> isequal(Symbolics.unwrap(p), userλ), + Symbolics.unwrap.(parameters(shadow)) + ) + @test !any(p -> isequal(Symbolics.unwrap(p), λ), Symbolics.unwrap.(parameters(shadow))) + rhs = only(equations(shadow)).rhs + # substituting the USER symbol does NOT resolve the blend: the result still + # depends on the sentinel λ (its endpoints differ) + after_user = Symbolics.substitute(rhs, Dict(__homotopy_λ => 1.0)) + @test !isequal( + Symbolics.substitute(after_user, Dict(λn => 1.0)), + Symbolics.substitute(after_user, Dict(λn => 0.0)) + ) + # only the sentinel λ resolves it: λ=1 → actual, λ=0 → simplified + @test isequal(Symbolics.substitute(rhs, Dict(λn => 1.0)), -x * __homotopy_λ) + @test isequal(Symbolics.substitute(rhs, Dict(λn => 0.0)), -x + 1) +end + +@testset "guard: a user symbol colliding with the sentinel name errors" begin + # A user symbol spelled EXACTLY like the reserved sentinel would otherwise be + # silently rewritten into a state/parameter access, dropping the λ-dependence + # of the residual. The loud guard turns that into an error. + @variables x(t) + @parameters __homotopy_λₘₜₖ + @named sys = System([D(x) ~ homotopy(-x * __homotopy_λₘₜₖ, -x + 1)], t) + csys = complete(sys) + @test_throws ArgumentError lower_homotopy(csys) +end + +@testset "broadcast: homotopy.(actual, simplified) lowers elementwise, one shared λ" begin + @variables x(t)[1:2] + actual = [atan(x[1] - 3), atan(x[2] - 5)] + simplified = [x[1], x[2]] + eqs = collect(0 .~ homotopy.(actual, simplified)) + @named sys = System(eqs, t) + csys = complete(sys) + @test has_any_homotopy(csys) # one homotopy node survives per element + shadow, λ = lower_homotopy(csys) + λn = Symbolics.wrap(λ) + @test !has_any_homotopy(shadow) # every element lowered + # the SAME λ resolves both elementwise blends, at both endpoints + for (i, eq) in enumerate(equations(shadow)) + @test isequal(Symbolics.substitute(eq.rhs, Dict(λn => 1.0)), actual[i]) + @test isequal(Symbolics.substitute(eq.rhs, Dict(λn => 0.0)), simplified[i]) + end +end + +@testset "nested homotopy lowers through lower_homotopy with one λ" begin + @variables w(t) + @named sys = System([0 ~ homotopy(homotopy(w^3, w^2), w)], t) + csys = complete(sys) + shadow, λ = lower_homotopy(csys) + λn = Symbolics.wrap(λ) + @test !has_any_homotopy(shadow) + rhs = only(equations(shadow)).rhs + # the SAME returned λ resolves both nesting levels + @test isequal(Symbolics.substitute(rhs, Dict(λn => 1.0)), w^3) + @test isequal(Symbolics.substitute(rhs, Dict(λn => 0.0)), w) + @test Symbolics.value(Symbolics.substitute(rhs, Dict(λn => 1.0, w => 2.0))) ≈ 8.0 + @test Symbolics.value(Symbolics.substitute(rhs, Dict(λn => 0.0, w => 2.0))) ≈ 2.0 +end + +@testset "nested homotopy rewrites recursively with the same λ" begin + # worker-level isolation test for `_rewrite_with_lambda`; the system-level + # nesting path through the public API is the testset above. + @variables w λtest + shadowex = ModelingToolkitBase._rewrite_with_lambda( + Symbolics.unwrap(homotopy(homotopy(w^3, w^2), w)), Symbolics.unwrap(λtest) + ) + @test Symbolics.value( + Symbolics.substitute( + Symbolics.wrap(shadowex), Dict(λtest => 1.0, w => 2.0) + ) + ) ≈ 8.0 + @test Symbolics.value( + Symbolics.substitute( + Symbolics.wrap(shadowex), Dict(λtest => 0.0, w => 2.0) + ) + ) ≈ 2.0 +end + +@testset "generate_homotopy_residual guards: raw nodes and time-dependence error" begin + # Both misuse modes were probed to fail SILENTLY with wrong numerics rather than + # erroring; the guards turn them loud. These lock that contract. + # + # Guard 1: a system still carrying raw homotopy(...) nodes (i.e. NOT run through + # `lower_homotopy`) would codegen a λ-independent residual via the numeric + # fallback — caught before any code is generated. + @variables x + @named raw = System([0 ~ homotopy(x^2 - 2, x)]) + craw = complete(raw) + @test has_any_homotopy(craw) + @test_throws ArgumentError ModelingToolkitBase.generate_homotopy_residual( + craw, ModelingToolkitBase.HOMOTOPY_LAMBDA + ) + # Guard 2: a correctly-lowered but time-dependent shadow — codegen would emit + # `f(u, p, t)` and silently misuse the λ slot as `t`. + @variables z(t) + @named tdep = System([D(z) ~ homotopy(-z, -z + 1)], t) + shadow, λ = lower_homotopy(complete(tdep)) + @test !has_any_homotopy(shadow) # lowered: the raw-node guard would pass + @test_throws ArgumentError ModelingToolkitBase.generate_homotopy_residual(shadow, λ) +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..256aba75b8 --- /dev/null +++ b/lib/ModelingToolkitBase/test/homotopy_omc_parity.jl @@ -0,0 +1,87 @@ +using ModelingToolkitBase +using SciMLBase +# narrow import: a blanket `using SymbolicIndexingInterface` would clash with +# ModelingToolkitBase's exported `observed` +using SymbolicIndexingInterface: state_values, parameter_values +using NonlinearSolve # HomotopySweep (re-exported) + inner nonlinear solvers +using Test + +# OMC parity fixture: the Modelica Buildings PressureDrop pattern +# +# m_flow ~ homotopy(sign(dp)*sqrt(abs(dp))*k, m_flow_nominal*dp/dp_nominal) +# +# Written in `var ~ homotopy(...)` form ON PURPOSE: with the variable itself as +# the LHS, `mtkcompile` eliminates `m_flow` into `observed` (a `0 ~ m_flow - ...` +# form would keep it as an unknown). That makes this the suite's ONLY coverage +# where the homotopy node reaches the swept residual through an OBSERVED equation +# — exercising the observed-rewriting in `lower_homotopy` and the +# observed-inlining in `generate_homotopy_residual` on a real Modelica-parity +# model. The single unknown `w` is pinned THROUGH the homotopy-observed `m_flow` +# by `w*|w| ~ m_flow^2/5`: `u*|u|` is monotone, so the root is unique, and at the +# actual (λ = 1) endpoint `m_flow = √5` gives `w = 1`, the exact OMC solution +# point (`w = 1`, `m_flow = √5`, `dp = 5`). +# +# Frozen reference (MTK CI does not run OMC): OMC v1.27, captured 2026-06-06 from +# ~/code/modelica-OMEdit/homotopy-tests/PressureDropParityFixture.mo — +# m_flow = sqrt(5) = 2.2360679774997896, at dp = 5, k = 1, m_flow_nominal = 1, +# dp_nominal = 5. +@testset "homotopy OMC parity — Buildings PressureDrop fixture (m_flow eliminated into observed)" begin + OMC_M_FLOW = 2.2360679774997896 # = √5, OMC v1.27 ref captured 2026-06-06 + OMC_W = 1.0 + PARITY_TOL = 1.0e-6 + + @variables m_flow w + @parameters dp k m_flow_nominal dp_nominal + basic_flow(dpv, kv) = sign(dpv) * sqrt(abs(dpv)) * kv + + @mtkcompile sys = System( + [ + m_flow ~ homotopy(basic_flow(dp, k), m_flow_nominal * dp / dp_nominal), + 0 ~ w * abs(w) - m_flow^2 / 5, + ]; + guesses = [w => 0.5, m_flow => 2.0] + ) + + # The key structural fact: m_flow is eliminated into observed, NOT an unknown + # — the elimination this fixture exists to guard. + @test !any(v -> isequal(v, m_flow), unknowns(sys)) + @test any(eq -> isequal(eq.lhs, m_flow), observed(sys)) + # ...and the homotopy node reaches the system through OBSERVED, not equations. + @test ModelingToolkitBase.has_any_homotopy(sys) + @test !ModelingToolkitBase.has_homotopy_in_equations(equations(sys)) + @test any(eq -> ModelingToolkitBase.has_homotopy(eq.rhs), observed(sys)) + + op = [w => 0.5, dp => 5.0, k => 1.0, m_flow_nominal => 1.0, dp_nominal => 5.0] + prob = HomotopyProblem(sys, op) + @test prob isa SciMLBase.HomotopyProblem + + # λ-dependence of the wired residual: homotopy sits ONLY in observed here, so + # differing endpoints come solely from observed-rewriting (`lower_homotopy`) + # + observed-inlining (`generate_homotopy_residual`) — the codegen path under + # guard. + u = copy(state_values(prob)) + p = parameter_values(prob) + if SciMLBase.isinplace(prob) + r0 = similar(u) + r1 = similar(u) + prob.f(r0, u, p, 0.0) + prob.f(r1, u, p, 1.0) + @test !(r0 ≈ r1) + else + @test !(prob.f(u, p, 0.0) ≈ prob.f(u, p, 1.0)) + end + + @testset "default sweep matches OMC" begin + sol = solve(prob; abstol = 1.0e-12, reltol = 1.0e-12) + @test SciMLBase.successful_retcode(sol) + @test abs(sol[m_flow] - OMC_M_FLOW) < PARITY_TOL + @test abs(sol[w] - OMC_W) < PARITY_TOL + end + + @testset "explicit HomotopySweep matches OMC" begin + sol = solve(prob, HomotopySweep(); abstol = 1.0e-12, reltol = 1.0e-12) + @test SciMLBase.successful_retcode(sol) + @test abs(sol[m_flow] - OMC_M_FLOW) < PARITY_TOL + @test abs(sol[w] - OMC_W) < PARITY_TOL + end +end diff --git a/lib/ModelingToolkitBase/test/homotopy_problem.jl b/lib/ModelingToolkitBase/test/homotopy_problem.jl new file mode 100644 index 0000000000..1f8bcbd7a6 --- /dev/null +++ b/lib/ModelingToolkitBase/test/homotopy_problem.jl @@ -0,0 +1,105 @@ +using ModelingToolkitBase +using ModelingToolkitBase: homotopy, has_any_homotopy, get_nonlinear_problem_type, + t_nounits as t, D_nounits as D +using SciMLBase +using NonlinearSolve # HomotopySweep (re-exported) + inner nonlinear solvers +using Test + +@testset "direct HomotopyProblem(sys, op): construction + residual endpoints" begin + @variables y + @mtkcompile sys = System([0 ~ homotopy(atan(y - 3), y)]) + prob = HomotopyProblem(sys, [y => 12.0]) + @test prob isa SciMLBase.HomotopyProblem + @test prob.λspan == (0.0, 1.0) # inherited SciMLBase default + @test !hasfield(typeof(prob), :homotopy_parameter) # merged API: λ is an arg, not a field + # the swept residual genuinely depends on λ (endpoints differ through the wrap) + @test !(prob.f(copy(prob.u0), prob.p, 0.0) ≈ prob.f(copy(prob.u0), prob.p, 1.0)) + # a custom λspan threads through to the constructed problem + @test HomotopyProblem(sys, [y => 12.0]; λspan = (0.0, 0.5)).λspan == (0.0, 0.5) + # expression = Val{true} is an explicit (documented) error, not silent misgen + @test_throws ArgumentError HomotopyProblem(sys, [y => 12.0]; expression = Val{true}) + + u = copy(prob.u0) + p = prob.p + @test only(prob.f(u, p, 0.0)) ≈ 12.0 atol = 1.0e-10 # simplified endpoint + @test only(prob.f(u, p, 1.0)) ≈ atan(9.0) atol = 1.0e-10 # actual endpoint + @test only(prob.f(u, p, 0.5)) ≈ 0.5 * 12.0 + 0.5 * atan(9.0) atol = 1.0e-10 # convex blend + # in-place residual variant f(du, u, p, λ): same wrapper, arity 4 + if SciMLBase.isinplace(prob) + du = similar(u) + prob.f(du, u, p, 0.5) + @test only(du) ≈ only(prob.f(u, p, 0.5)) atol = 1.0e-10 + end +end + +@testset "AbstractNonlinearProblem(sys, op) selects HomotopyProblem ⟺ homotopy present" begin + @variables y + @mtkcompile hsys = System([0 ~ homotopy(atan(y - 3), y)]) + @test get_nonlinear_problem_type(hsys) === SciMLBase.HomotopyProblem + @test AbstractNonlinearProblem(hsys, [y => 12.0]) isa SciMLBase.HomotopyProblem + + @variables z + @mtkcompile nsys = System([0 ~ atan(z - 3)]) + @test get_nonlinear_problem_type(nsys) === SciMLBase.NonlinearProblem + np = AbstractNonlinearProblem(nsys, [z => 12.0]) + @test np isa SciMLBase.NonlinearProblem + @test !(np isa SciMLBase.HomotopyProblem) + + # HomotopyProblem on a system without `homotopy` is a loud error, not a + # silent λ-independent problem. + @test_throws ArgumentError HomotopyProblem(nsys, [z => 12.0]) +end + +@testset "constructor converts a time-dependent system (ODE -> NonlinearSystem)" begin + # A time-dependent system exercises the is_time_dependent -> NonlinearSystem + # conversion in the constructor (D(x) substituted to 0), which must still leave + # the homotopy node intact and produce a swept HomotopyProblem. + @variables x(t) y(t) + @mtkcompile sys = System( + [D(x) ~ -x, 0 ~ homotopy(atan(y - 3), y)], t; guesses = [y => 12.0] + ) + prob = HomotopyProblem(sys, [x => 1.0, y => 12.0]) + @test prob isa SciMLBase.HomotopyProblem + # residual is λ-dependent through the conversion + lowering + u = copy(prob.u0) + @test !(prob.f(u, prob.p, 0.0) ≈ prob.f(u, prob.p, 1.0)) + sol = solve(prob) + @test SciMLBase.successful_retcode(sol) + @test sol[y] ≈ 3.0 atol = 1.0e-6 +end + +@testset "e2e rescue: continuation walks an out-of-basin guess to the true root" begin + @variables y + @mtkcompile sys = System([0 ~ homotopy(atan(y - 3), y)]) + prob = HomotopyProblem(sys, [y => 12.0]) + # a `HomotopyProblem` with no algorithm defaults to the continuation sweep + sol = solve(prob) + @test SciMLBase.successful_retcode(sol) + @test sol[y] ≈ 3.0 atol = 1.0e-6 + # an explicit `HomotopySweep` gives the same answer + sol2 = solve(prob, HomotopySweep()) + @test sol2[y] ≈ 3.0 atol = 1.0e-6 +end + +@testset "structured p passes through the sweep untouched" begin + @variables y + @parameters c = 3.0 + @mtkcompile sys = System([0 ~ homotopy(atan(y - c), y)]) + prob = HomotopyProblem(sys, [y => 12.0, c => 3.0]) + @test prob.p isa ModelingToolkitBase.MTKParameters + sol = solve(prob) + @test sol[y] ≈ 3.0 atol = 1.0e-6 + @test sol.ps[c] == 3.0 # user parameter unscathed by the sweep +end + +@testset "two homotopy() calls share the single sweep" begin + @variables y z + @mtkcompile sys = System( + [0 ~ homotopy(atan(y - 3), y), 0 ~ homotopy(atan(z + 2), z)] + ) + prob = HomotopyProblem(sys, [y => 12.0, z => -12.0]) + @test prob isa SciMLBase.HomotopyProblem + sol = solve(prob) + @test sol[y] ≈ 3.0 atol = 1.0e-6 + @test sol[z] ≈ -2.0 atol = 1.0e-6 +end diff --git a/lib/ModelingToolkitBase/test/runtests.jl b/lib/ModelingToolkitBase/test/runtests.jl index ab61d90389..6898d461d6 100644 --- a/lib/ModelingToolkitBase/test/runtests.jl +++ b/lib/ModelingToolkitBase/test/runtests.jl @@ -93,6 +93,9 @@ end @safetestset "SDESystem Test" include("sdesystem.jl") @safetestset "DDESystem Test" include("dde.jl") @safetestset "NonlinearSystem Test" include("nonlinearsystem.jl") + @safetestset "Homotopy lowering" include("homotopy_lowering.jl") + @safetestset "Homotopy problem construction & sweep" include("homotopy_problem.jl") + @safetestset "Homotopy OMC parity" include("homotopy_omc_parity.jl") @safetestset "PDE Construction Test" include("pdesystem.jl") @safetestset "JumpSystem Test" include("jumpsystem.jl") @safetestset "Poissonians Test" include("poissonians.jl")