Add conservative kwarg in structural_transformation

`conservative=true` limits tearing to only solve for trivial linear systems
where the coefficient has the absolute value of 1. This is useful for debugging
numerical stability issues after tearing. By default, we set
`conservative=false`.

Author: YingboMa

src/structural_transformation/pantelides.jl

@@ -123,7 +123,8 @@ end
 
 Perform Pantelides algorithm.
 """
-function pantelides!(state::TransformationState; finalize = true, maxiters = 8000)
+function pantelides!(
+        state::TransformationState; finalize = true, maxiters = 8000, kwargs...)
     @unpack graph, solvable_graph, var_to_diff, eq_to_diff = state.structure
     neqs = nsrcs(graph)
     nvars = nv(var_to_diff)
@@ -181,7 +182,7 @@ function pantelides!(state::TransformationState; finalize = true, maxiters = 800
                 ecolor[eq] || continue
                 # introduce a new equation
                 neqs += 1
-                eq_derivative!(state, eq)
+                eq_derivative!(state, eq; kwargs...)
             end
 
             for var in eachindex(vcolor)

src/structural_transformation/partial_state_selection.jl

@@ -173,7 +173,7 @@ function dummy_derivative_graph!(state::TransformationState, jac = nothing;
         state_priority = nothing, log = Val(false), kwargs...)
     state.structure.solvable_graph === nothing && find_solvables!(state; kwargs...)
     complete!(state.structure)
-    var_eq_matching = complete(pantelides!(state))
+    var_eq_matching = complete(pantelides!(state; kwargs...))
     dummy_derivative_graph!(state.structure, var_eq_matching, jac, state_priority, log)
 end
 

src/structural_transformation/symbolics_tearing.jl

@@ -56,7 +56,7 @@ function eq_derivative_graph!(s::SystemStructure, eq::Int)
     return eq_diff
 end
 
-function eq_derivative!(ts::TearingState{ODESystem}, ieq::Int)
+function eq_derivative!(ts::TearingState{ODESystem}, ieq::Int; kwargs...)
     s = ts.structure
 
     eq_diff = eq_derivative_graph!(s, ieq)
@@ -75,7 +75,8 @@ function eq_derivative!(ts::TearingState{ODESystem}, ieq::Int)
         add_edge!(s.graph, eq_diff, s.var_to_diff[var])
     end
     s.solvable_graph === nothing ||
-        find_eq_solvables!(ts, eq_diff; may_be_zero = true, allow_symbolic = false)
+        find_eq_solvables!(
+            ts, eq_diff; may_be_zero = true, allow_symbolic = false, kwargs...)
 
     return eq_diff
 end

src/structural_transformation/utils.jl

@@ -181,7 +181,9 @@ end
 
 function find_eq_solvables!(state::TearingState, ieq, to_rm = Int[], coeffs = nothing;
         may_be_zero = false,
-        allow_symbolic = false, allow_parameter = true, kwargs...)
+        allow_symbolic = false, allow_parameter = true,
+        conservative = false,
+        kwargs...)
     fullvars = state.fullvars
     @unpack graph, solvable_graph = state.structure
     eq = equations(state)[ieq]
@@ -220,6 +222,7 @@ function find_eq_solvables!(state::TearingState, ieq, to_rm = Int[], coeffs = no
             coeffs === nothing || push!(coeffs, convert(Int, a))
         else
             all_int_vars = false
+            conservative && continue
         end
         if a != 0
             add_edge!(solvable_graph, ieq, j)

src/systems/systems.jl

@@ -9,8 +9,10 @@ $(SIGNATURES)
 Structurally simplify algebraic equations in a system and compute the
 topological sort of the observed equations. When `simplify=true`, the `simplify`
 function will be applied during the tearing process. It also takes kwargs
-`allow_symbolic=false` and `allow_parameter=true` which limits the coefficient
-types during tearing.
+`allow_symbolic=false`, `allow_parameter=true`, and `conservative=false` which
+limits the coefficient types during tearing. In particular, `conservative=true`
+limits tearing to only solve for trivial linear systems where the coefficient
+has the absolute value of ``1``.
 
 The optional argument `io` may take a tuple `(inputs, outputs)`.
 This will convert all `inputs` to parameters and allow them to be unconnected, i.e.,

src/systems/systemstructure.jl

@@ -691,15 +691,18 @@ function _structural_simplify!(state::TearingState, io; simplify = false,
         ModelingToolkit.check_consistency(state, orig_inputs)
     end
     if fully_determined && dummy_derivative
-        sys = ModelingToolkit.dummy_derivative(sys, state; simplify, mm, check_consistency)
+        sys = ModelingToolkit.dummy_derivative(
+            sys, state; simplify, mm, check_consistency, kwargs...)
     elseif fully_determined
         var_eq_matching = pantelides!(state; finalize = false, kwargs...)
         sys = pantelides_reassemble(state, var_eq_matching)
         state = TearingState(sys)
         sys, mm = ModelingToolkit.alias_elimination!(state; kwargs...)
-        sys = ModelingToolkit.dummy_derivative(sys, state; simplify, mm, check_consistency)
+        sys = ModelingToolkit.dummy_derivative(
+            sys, state; simplify, mm, check_consistency, kwargs...)
     else
-        sys = ModelingToolkit.tearing(sys, state; simplify, mm, check_consistency)
+        sys = ModelingToolkit.tearing(
+            sys, state; simplify, mm, check_consistency, kwargs...)
     end
     fullunknowns = [map(eq -> eq.lhs, observed(sys)); unknowns(sys)]
     @set! sys.observed = ModelingToolkit.topsort_equations(observed(sys), fullunknowns)