diff --git a/Project.toml b/Project.toml index 00f22c3..e09b30b 100644 --- a/Project.toml +++ b/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" diff --git a/lib/ModelingToolkitTearing/Project.toml b/lib/ModelingToolkitTearing/Project.toml index 38907e6..a37e65e 100644 --- a/lib/ModelingToolkitTearing/Project.toml +++ b/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 "] [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" diff --git a/lib/ModelingToolkitTearing/src/stateselection_interface.jl b/lib/ModelingToolkitTearing/src/stateselection_interface.jl index e86f8ca..20d5d2e 100644 --- a/lib/ModelingToolkitTearing/src/stateselection_interface.jl +++ b/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 diff --git a/lib/ModelingToolkitTearing/src/tearingstate.jl b/lib/ModelingToolkitTearing/src/tearingstate.jl index beef0a8..413ac4a 100644 --- a/lib/ModelingToolkitTearing/src/tearingstate.jl +++ b/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 """ diff --git a/lib/ModelingToolkitTearing/test/runtests.jl b/lib/ModelingToolkitTearing/test/runtests.jl index 6e488a3..c7260c5 100644 --- a/lib/ModelingToolkitTearing/test/runtests.jl +++ b/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 diff --git a/src/partial_state_selection.jl b/src/partial_state_selection.jl index f7ed304..975cc42 100644 --- a/src/partial_state_selection.jl +++ b/src/partial_state_selection.jl @@ -1,4 +1,5 @@ using BipartiteGraphs: Unassigned, maximal_matching +using LinearAlgebra: norm function partial_state_selection_graph!(state::TransformationState) var_eq_matching = complete(pantelides!(state)) @@ -175,14 +176,15 @@ function partial_state_selection_graph!(structure::SystemStructure, var_eq_match end function dummy_derivative_graph!(state::TransformationState, jac = nothing; - state_priority = nothing, log = Val(false), kwargs...) + state_priority = nothing, log = Val(false), numjac = nothing, kwargs...) state.structure.solvable_graph === nothing && find_solvables!(state; kwargs...) complete!(state.structure) var_eq_matching = complete(pantelides!(state; kwargs...)) # NOTE: `get_mm` must be queried after `pantelides!`, which extends the # linear subsystem matrix with differentiated rows (`eq_derivative!`). dummy_derivative_graph!( - state.structure, var_eq_matching, jac, state_priority, log; mm = get_mm(state)) + state.structure, var_eq_matching, jac, state_priority, log; + mm = get_mm(state), numjac) end struct DummyDerivativeSummary @@ -190,6 +192,96 @@ struct DummyDerivativeSummary state_priority::Vector{Vector{Float64}} end +""" + $TYPEDSIGNATURES + +Greedily select, in the given order, the columns of `J` that are numerically linearly +independent of the previously selected ones. Rows are first scaled to unit maximum +magnitude, which does not change which column sets are singular but makes the test +independent of the units of the equations. A column is skipped when its norm is below +`rtol` times the largest column norm, and accepted when the norm of its component +orthogonal to the span of the accepted columns exceeds `rtol` times its own norm. At most +`size(J, 1)` columns are selected. +""" +function numerically_independent_columns( + J::AbstractMatrix{<:Real}, cols; rtol::Float64 = sqrt(eps(Float64)) + ) + m = size(J, 1) + rowscale = [(s = maximum(abs, @view J[i, :]); s > 0 ? inv(s) : 1.0) for i in 1:m] + maxnorm = maximum((norm(rowscale .* @view J[:, c]) for c in axes(J, 2)); init = 0.0) + basis = Vector{Vector{Float64}}() + accepted = Int[] + r = zeros(m) + for c in cols + r .= rowscale .* @view J[:, c] + nc = norm(r) + nc > rtol * maxnorm || continue + # two passes of modified Gram-Schmidt for numerical stability + for _ in 1:2, q in basis + r .-= (q' * r) .* q + end + nr = norm(r) + nr > rtol * nc || continue + push!(basis, r ./ nr) + push!(accepted, c) + length(accepted) == m && break + end + return accepted +end + +""" + $TYPEDSIGNATURES + +Select dummy derivatives among `vars` for the equations `eqs` by structural rank: walk +`vars` in order and accept a variable when an augmenting path to an unmatched equation +exists. Accepted variables are appended to `dummy_derivatives`. Returns the rank found. +""" +function structural_rank_selection!( + dummy_derivatives::Vector{Int}, vars::Vector{Int}, eqs::Vector{Int}, + rank_matching::Matching, invgraph, eqs_set::BitSet, eqcolor::BitVector, nrows::Int + ) + empty!(eqs_set) + union!(eqs_set, eqs) + rank = 0 + for var in vars + eqcolor .= false + # We need `invgraph` here because we are matching from + # variables to equations. + pathfound = construct_augmenting_path!(rank_matching, invgraph, var, + Base.Fix2(in, eqs_set), eqcolor) + pathfound || continue + push!(dummy_derivatives, var) + rank += 1 + rank == nrows && break + end + fill!(rank_matching, unassigned) + return rank +end + +""" + $TYPEDSIGNATURES + +Return the highest-derivative dummy-derivative candidates incident to the equations `eqs` +that are not in `vars` and that `var_eq_matching` leaves unassigned, i.e. that would become +differential states. Sorted by variable index. +""" +function unassigned_candidates( + structure::SystemStructure, var_eq_matching, eqs::Vector{Int}, vars::Vector{Int} + ) + (; var_to_diff, graph) = structure + diff_to_var = invview(var_to_diff) + extra = Int[] + for eq in eqs, var in 𝑠neighbors(graph, eq) + var_eq_matching[var] === unassigned || continue + var_to_diff[var] === nothing || continue + diff_to_var[var] !== nothing && is_present(structure, var) || continue + var in vars && continue + var in extra && continue + push!(extra, var) + end + return sort!(extra) +end + """ $TYPEDSIGNATURES @@ -202,12 +294,25 @@ Perform the dummy derivatives algorithm. return `nothing`. - `state_priority` is a function taking the index of a variable and returning its priority. Higher priority variables are more likely to be chosen as states. + +# Keyword Arguments + +- `numjac` is `nothing` or a function taking a list of equation and variable indices and + returning the jacobian for the same evaluated numerically at the initial point, as a + matrix of finite reals, or `nothing` if it cannot be evaluated. It is only consulted for + SCCs without an integer jacobian, after the structural selection: when the variables + selected as dummy derivatives have a numerically rank-deficient jacobian at the initial + point, the candidates are reordered so that a numerically independent set (found greedily + in priority order) is selected instead. If the candidates matched in the SCC admit no + such set, the highest derivatives incident to its equations that the matching left + unassigned are added to the candidates. When no such set exists at all, a warning is + emitted and the structural selection is kept. """ function dummy_derivative_graph!( structure::SystemStructure, var_eq_matching, jac = nothing, state_priority = nothing, ::Val{log} = Val(false); tearing_alg::TearingAlgorithm = DummyDerivativeTearing(), - mm = nothing, kwargs...) where {log} + mm = nothing, numjac = nothing, kwargs...) where {log} (; eq_to_diff, var_to_diff, graph) = structure diff_to_eq = invview(eq_to_diff) diff_to_var = invview(var_to_diff) @@ -314,21 +419,60 @@ function dummy_derivative_graph!( push!(dummy_derivatives, vars[col_order[i]]) end else - empty!(eqs_set) - union!(eqs_set, eqs) - rank = 0 - for var in vars - eqcolor .= false - # We need `invgraph` here because we are matching from - # variables to equations. - pathfound = construct_augmenting_path!(rank_matching, invgraph, var, - Base.Fix2(in, eqs_set), eqcolor) - pathfound || continue - push!(dummy_derivatives, var) - rank += 1 - rank == nrows && break + n_before = length(dummy_derivatives) + rank = structural_rank_selection!( + dummy_derivatives, vars, eqs, rank_matching, invgraph, eqs_set, + eqcolor, nrows) + # The structural choice may be singular at the initial point (e.g. a + # constraint jacobian column that vanishes there). Only then, reorder the + # candidates by a numerically independent set and select again. This is + # done at every differentiation level: the highest level also contains + # the differentiated kinematic equations, which can hide a singular + # choice of coordinates. The order of `vars` is inherited by the lower + # levels. + Jn = (numjac !== nothing && rank == nrows) ? numjac(eqs, vars) : nothing + if Jn !== nothing && size(Jn) == (nrows, length(vars)) + chosen_cols = Int[ + findfirst(==(v), vars) + for v in @view dummy_derivatives[(n_before + 1):end] + ] + if length(numerically_independent_columns(Jn, chosen_cols)) != nrows + pivots = numerically_independent_columns(Jn, eachindex(vars)) + if length(pivots) != nrows && isfirst + # The matched candidates alone are singular. Widen the pool + # with the highest derivatives incident to `eqs` that the + # matching left unassigned (i.e. would become states), so + # that one of them can be eliminated instead. + extra = unassigned_candidates(structure, var_eq_matching, eqs, vars) + sort!(extra; by = var -> ( + state_priority === nothing ? 0 : extended_sp(var), + cranks === nothing ? 0 : cranks[var])) + Jn_ext = isempty(extra) ? nothing : numjac(eqs, vcat(vars, extra)) + if Jn_ext !== nothing && size(Jn_ext) == (nrows, length(vars) + length(extra)) + pivots_ext = numerically_independent_columns(Jn_ext, axes(Jn_ext, 2)) + if length(pivots_ext) == nrows + append!(vars, extra) + pivots = pivots_ext + end + end + end + @debug "Dummy derivative selection singular at the initial point" eqs = repr(eqs) vars = repr(vars) chosen_cols = repr(chosen_cols) pivots = repr(pivots) + if length(pivots) == nrows + perm = vcat(pivots, setdiff(eachindex(vars), pivots)) + permute!(vars, perm) + if state_priority !== nothing && isfirst + var_dummy_scc[end] = copy(vars) + var_state_priority[end] = extended_sp.(vars) + end + resize!(dummy_derivatives, n_before) + rank = structural_rank_selection!( + dummy_derivatives, vars, eqs, rank_matching, invgraph, + eqs_set, eqcolor, nrows) + else + @warn "The dummy derivatives selected for a set of $(nrows) differentiated equations are numerically singular at the initial point, and no non-singular selection was found among the candidates at this differentiation level. Consider changing the initial point or the `state_priority` of the involved variables." + end + end end - fill!(rank_matching, unassigned) end if rank != nrows @warn "The DAE system is singular!" diff --git a/test/runtests.jl b/test/runtests.jl index 1a03238..f886a40 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -72,3 +72,21 @@ include("carpanzano_tearing.jl") @test mm2.row_cols == [[2]] @test mm2.row_vals == [[1]] end + +@testset "`numerically_independent_columns`" begin + J = [1.0 2.0 0.0; 0.0 0.0 1.0] + @test SSel.numerically_independent_columns(J, 1:3) == [1, 3] + @test SSel.numerically_independent_columns(J, [2, 1, 3]) == [2, 3] + @test SSel.numerically_independent_columns(J, [1, 2]) == [1] + # never more columns than rows + @test SSel.numerically_independent_columns([1.0 0.0 1.0], 1:3) == [1] + # zero columns are skipped + @test SSel.numerically_independent_columns([0.0 1.0; 0.0 1.0], 1:2) == [2] + # dependence below the tolerance counts as dependence + @test SSel.numerically_independent_columns([1.0 1.0 + 1e-12; 1.0 1.0], 1:2) == [1] + @test SSel.numerically_independent_columns([1.0 1.0 + 1e-12; 1.0 1.0], 1:2; rtol = 1e-14) == [1, 2] + # rows are scaled, so a badly scaled but regular matrix keeps its rank + @test SSel.numerically_independent_columns([1e8 1e8; 1.0 2.0], 1:2) == [1, 2] + # columns negligible against the largest one are skipped + @test SSel.numerically_independent_columns([1.0 1e-12; 1.0 0.0], 1:2) == [1] +end