Use the initial point to avoid singular dummy-derivative selections and solve directions

Project.toml

@@ -1,7 +1,7 @@
 name = "StateSelection"
 uuid = "64909d44-ed92-46a8-bbd9-f047dfbdc84b"
 authors = ["JuliaHub", "Inc. and other contributors"]
-version = "1.11.1"
+version = "1.12.0"
 
 [deps]
 BipartiteGraphs = "caf10ac8-0290-4205-88aa-f15908547e8d"

lib/ModelingToolkitTearing/Project.toml

@@ -1,6 +1,6 @@
 name = "ModelingToolkitTearing"
 uuid = "6bb917b9-1269-42b9-9f7c-b0dca72083ab"
-version = "1.20.7"
+version = "1.21.0"
 authors = ["Aayush Sabharwal <aayush.sabharwal@juliahub.com>"]
 
 [deps]
@@ -30,7 +30,7 @@ ForwardDiff = "1.3"
 Graphs = "1"
 LinearAlgebra = "1"
 ModelingToolkit = "11"
-ModelingToolkitBase = "1.57.1"
+ModelingToolkitBase = "1.69.0"
 Moshi = "0.3.7"
 OffsetArrays = "1"
 OrderedCollections = "1.8.1, 2"

lib/ModelingToolkitTearing/src/stateselection_interface.jl

@@ -315,6 +315,22 @@ function _check_allow_symbolic_parameter(
 end
 
 
+"""
+    $TYPEDSIGNATURES
+
+Check whether the coefficient `coeff` (a symbolic or array thereof) evaluates to exactly
+zero at the initial point of `state.sys`, see [`evaluate_at_initial_point`](@ref). An
+array coefficient counts as zero if any of its entries is. Coefficients that cannot be
+evaluated are not zero.
+"""
+function is_zero_at_initial_point(state::TearingState, coeff)
+    if coeff isa AbstractArray
+        return any(Base.Fix1(is_zero_at_initial_point, state), coeff)
+    end
+    val = evaluate_at_initial_point(state, coeff)
+    return val !== nothing && iszero(val)
+end
+
 const _SUPPORTS_NEED_REMAINDER = isdefined(Symbolics, :SUPPORTS_LINEAR_EXPANDER_NEED_REMAINDER)
 
 function StateSelection.find_eq_solvables!(state::TearingState, ieq, to_rm = Int[], coeffs = nothing;
@@ -359,6 +375,11 @@ function StateSelection.find_eq_solvables!(state::TearingState, ieq, to_rm = Int
             if !_check_allow_symbolic_parameter(state, a, allow_symbolic, allow_parameter; fullvars_set)
                 continue
             end
+            # A coefficient that vanishes at the initial point (e.g. `sin(ω*t)` at `t = 0`)
+            # must not be divided by, like an expression in `maybe_zeros`.
+            if !allow_symbolic && is_zero_at_initial_point(state, a)
+                continue
+            end
             add_edge!(solvable_graph, ieq, j)
             continue
         end

lib/ModelingToolkitTearing/src/tearingstate.jl

@@ -112,12 +112,81 @@ mutable struct TearingState <: StateSelection.TransformationState{System}
     and put into `additional_observed`.
     """
     analytical_derivatives::Dict{SymbolicT, SymbolicT}
+    """
+    Lazily built substituter evaluating expressions at the initial point of `sys`, see
+    [`initial_point_substituter`](@ref). `nothing` until first used.
+    """
+    initial_point::Base.RefValue{Any}
 end
 
 function Base.show(io::IO, state::TearingState)
     print(io, "TearingState of ", typeof(state.sys))
 end
 
+"""
+    $TYPEDSIGNATURES
+
+Build a substituter that evaluates expressions at the initial point of `sys`: parameters
+take their bindings, unknowns their initial conditions, and the independent variable the
+start of `get_tspan(sys)` when the system has a `tspan`. Values given in `initial_point`
+(an iterable of `variable => value` pairs or a dict, e.g. the `u0` map later passed to the
+problem constructor) take precedence over all of these. Guesses are deliberately not
+used: they are starting points for the initialization solver, not the initial point.
+Variables without a value stay symbolic, so an expression depending on them does not
+evaluate to a number.
+"""
+function initial_point_substituter(sys::System; initial_point = nothing)
+    defs = copy(parent(bindings(sys)))
+    MTKBase.left_merge!(defs, initial_conditions(sys))
+    iv = MTKBase.get_iv(sys)
+    tspan = MTKBase.get_tspan(sys)
+    if iv !== nothing && tspan !== nothing && first(tspan) isa Real
+        defs[iv] = BSImpl.Const{VartypeT}(first(tspan))
+    end
+    if initial_point !== nothing
+        for (k, v) in initial_point
+            defs[unwrap(k)] = v isa SymbolicT ? v : BSImpl.Const{VartypeT}(v)
+        end
+    end
+    filter!(Base.Fix2(!==, MTKBase.COMMON_MISSING) ∘ last, defs)
+    return Symbolics.FixpointSubstituter{true}(
+        MTKBase.AADSubWrapper(defs); maxiters = clamp(length(defs), 10, 1000),
+        warn_maxiters = false
+    )
+end
+
+"""
+    $TYPEDSIGNATURES
+
+Evaluate `ex` at the initial point of `state.sys` (see [`initial_point_substituter`](@ref)).
+Returns the value as a `Float64`, or `nothing` if `ex` does not reduce to a finite real
+number there, including when evaluating it throws (e.g. a `DomainError`).
+
+The substituter is built on first use and cached in `state.initial_point`. This relies on
+the bindings, initial conditions and `tspan` of `state.sys` not changing during structural
+simplification, which holds for all passes that replace `state.sys`.
+"""
+function evaluate_at_initial_point(state::TearingState, ex)
+    ex = unwrap(ex)
+    if !(ex isa SymbolicT)
+        return ex isa Real ? Float64(ex) : nothing
+    end
+    subber = state.initial_point[]
+    if subber === nothing
+        subber = state.initial_point[] = initial_point_substituter(state.sys)
+    end
+    val = try
+        subber(ex)
+    catch err
+        @debug "Evaluating an expression at the initial point failed" ex err
+        return nothing
+    end
+    SU.isconst(val) || return nothing
+    val = SU.unwrap_const(val)
+    val isa Real || return nothing
+    return Float64(val)
+end
+
 StateSelection.has_equations(::TearingState) = true
 StateSelection.equations(ts::TearingState) = equations(ts)
 
@@ -510,7 +579,8 @@ function TearingState(sys::System, source_info::Union{Nothing, MTKBase.EquationS
                                 canonical_ranks, false)
     return TearingState(sys, fullvars, structure, Equation[], param_derivative_map,
                         no_deriv_params, original_eqs, Equation[], falses(length(fullvars)),
-                        typeof(sys)[], sources, nothing, Dict{SymbolicT, SymbolicT}())
+                        typeof(sys)[], sources, nothing, Dict{SymbolicT, SymbolicT}(),
+                        Base.RefValue{Any}(nothing))
 end
 
 """

lib/ModelingToolkitTearing/test/runtests.jl

@@ -794,3 +794,52 @@ end
     @test unwrap(xd) in vars
     @test !(unwrap(xd(k - 1)) in vars)
 end
+
+@testset "`evaluate_at_initial_point`" begin
+    @variables x(t) [guess = 2.0] y(t) z(t)
+    @parameters p = 3.0
+    @named sys = System([D(x) ~ y + z, D(z) ~ p * x], t; initial_conditions = [y => 4.0])
+    ts = TearingState(sys)
+    @test MTKTearing.evaluate_at_initial_point(ts, y * p) == 12.0
+    @test MTKTearing.evaluate_at_initial_point(ts, y - 4) == 0.0
+    @test MTKTearing.evaluate_at_initial_point(ts, 1.5) == 1.5
+    # guesses are not the initial point, and `z` has no value at all
+    @test MTKTearing.evaluate_at_initial_point(ts, x + p) === nothing
+    @test MTKTearing.evaluate_at_initial_point(ts, z + 1) === nothing
+    # without a `tspan` the initial time is unknown
+    @test MTKTearing.evaluate_at_initial_point(ts, sin(t)) === nothing
+    @test MTKTearing.is_zero_at_initial_point(ts, y - 4)
+    @test MTKTearing.is_zero_at_initial_point(ts, [p, y - 4])
+    @test !MTKTearing.is_zero_at_initial_point(ts, z)
+    # evaluation errors are not zero either
+    @test !MTKTearing.is_zero_at_initial_point(ts, sqrt(-y))
+    # the start of `tspan` is the initial time
+    @named sys = System([D(x) ~ y + z, D(z) ~ p * x], t; tspan = (1.0, 2.0))
+    ts = TearingState(sys)
+    @test MTKTearing.evaluate_at_initial_point(ts, 2t) == 2.0
+    # an explicit initial point takes precedence over everything
+    ts.initial_point[] = MTKTearing.initial_point_substituter(
+        ts.sys; initial_point = [x => 1.0, z => 2.0, t => 0.0]
+    )
+    @test MTKTearing.evaluate_at_initial_point(ts, x + z + t) == 3.0
+end
+
+@testset "coefficients that vanish at the initial point are not divided by" begin
+    @variables x(t) y(t)
+    @parameters p = 1.0 q = 1.0
+    # `sin(t)` is zero at `t = 0`, so `y` cannot be solved for from the second equation
+    @mtkcompile sys = System([D(x) ~ y, sin(t) * y ~ x], t; tspan = (0.0, 2.0))
+    @test issetequal(unknowns(sys), [x, y])
+    @test isempty(observables(sys))
+    # ... but it is not zero at `t = 1`, and unknown without a `tspan`
+    @mtkcompile sys = System([D(x) ~ y, sin(t) * y ~ x], t; tspan = (1.0, 2.0))
+    @test issetequal(unknowns(sys), [x])
+    @test issetequal(observables(sys), [y])
+    @mtkcompile sys = System([D(x) ~ y, sin(t) * y ~ x], t)
+    @test issetequal(unknowns(sys), [x])
+    # a combination of parameters that vanishes at their default values
+    @mtkcompile sys = System([D(x) ~ y, (p - q) * y ~ x], t)
+    @test issetequal(unknowns(sys), [x, y])
+    @mtkcompile sys = System([D(x) ~ y, (p + q) * y ~ x], t)
+    @test issetequal(unknowns(sys), [x])
+end