feat: add allow_algebraic, default allow_symbolic to false

Author: AayushSabharwal

src/structural_transformation/partial_state_selection.jl

@@ -170,11 +170,13 @@ 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.structure.solvable_graph === nothing && find_solvables!(state; kwargs...)
+        state_priority = nothing, log = Val(false), allow_symbolic = false, kwargs...)
+    state.structure.solvable_graph === nothing &&
+        find_solvables!(state; allow_symbolic, kwargs...)
     complete!(state.structure)
-    var_eq_matching = complete(pantelides!(state; kwargs...))
-    dummy_derivative_graph!(state.structure, var_eq_matching, jac, state_priority, log)
+    var_eq_matching = complete(pantelides!(state; allow_symbolic, kwargs...))
+    dummy_derivative_graph!(
+        state.structure, var_eq_matching, jac, state_priority, log; allow_symbolic)
 end
 
 struct DummyDerivativeSummary
@@ -184,7 +186,7 @@ end
 
 function dummy_derivative_graph!(
         structure::SystemStructure, var_eq_matching, jac = nothing,
-        state_priority = nothing, ::Val{log} = Val(false)) where {log}
+        state_priority = nothing, ::Val{log} = Val(false); allow_symbolic = false) where {log}
     @unpack eq_to_diff, var_to_diff, graph = structure
     diff_to_eq = invview(eq_to_diff)
     diff_to_var = invview(var_to_diff)
@@ -342,7 +344,12 @@ function dummy_derivative_graph!(
         @warn "The number of dummy derivatives ($n_dummys) does not match the number of differentiated equations ($n_diff_eqs)."
     end
 
-    ret = tearing_with_dummy_derivatives(structure, BitSet(dummy_derivatives))
+    dummy_derivatives_set = BitSet(dummy_derivatives)
+    if !allow_symbolic
+        make_differential_denominators_unsolvable!(structure, dummy_derivatives_set)
+    end
+
+    ret = tearing_with_dummy_derivatives(structure, dummy_derivatives_set)
     if log
         (ret..., DummyDerivativeSummary(var_dummy_scc, var_state_priority))
     else

src/structural_transformation/symbolics_tearing.jl

@@ -992,8 +992,12 @@ end
     ndims = ndims(arr)
 end
 
-function tearing(state::TearingState; kwargs...)
-    state.structure.solvable_graph === nothing && find_solvables!(state; kwargs...)
+function tearing(state::TearingState; allow_symbolic = false, kwargs...)
+    state.structure.solvable_graph === nothing &&
+        find_solvables!(state; allow_symbolic, kwargs...)
+    if !allow_symbolic
+        make_differential_denominators_unsolvable!(state.structure)
+    end
     complete!(state.structure)
     tearing_with_dummy_derivatives(state.structure, ())
 end

src/structural_transformation/utils.jl

@@ -275,6 +275,22 @@ function find_eq_solvables!(state::TearingState, ieq::Int, to_rm = Int[], coeffs
     all_int_vars, term
 end
 
+"""
+    $(TYPEDSIGNATURES)
+
+Remove edges in `structure.solvable_graph` that require differential variables in the
+denominator to solve. `additional_algevars` is a collection of integers corresponding to
+differential variables that should be considered as algebraic for the purpose of this
+transformation.
+"""
+function make_differential_denominators_unsolvable!(
+        structure::SystemStructure, additional_algevars = ())
+    for ((eqi, vari), denoms) in structure.denominators
+        all(i -> isalgvar(structure, i) || i in additional_algevars, denoms) && continue
+        rem_edge!(structure.solvable_graph, eqi, vari)
+    end
+end
+
 function find_solvables!(state::TearingState; kwargs...)
     @assert state.structure.solvable_graph === nothing
     eqs = equations(state)

src/systems/systems.jl

@@ -18,21 +18,30 @@ Structurally simplify algebraic equations in a system and compute the
 topological sort of the observed equations in `sys`.
 
 ### Optional Arguments:
-+ optional argument `io` may take a tuple `(inputs, outputs)`. This will convert all `inputs` to parameters and allow them to be unconnected, i.e., simplification will allow models where `n_unknowns = n_equations - n_inputs`.
++ optional argument `io` may take a tuple `(inputs, outputs)`. This will convert all
+  `inputs` to parameters and allow them to be unconnected, i.e., simplification will
+  allow models where `n_unknowns = n_equations - n_inputs`.
 
 ### Optional Keyword Arguments:
-+ When `simplify=true`, the `simplify` function will be applied during the tearing process.
-+ `allow_symbolic=true`, `allow_parameter=true`, and `conservative=false` limit 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``.
-+ `fully_determined=true` controls whether or not an error will be thrown if the number of equations don't match the number of inputs, outputs, and equations.
++ When `simplify=true`, the `simplify` function will be applied during the tearing
+  process.
++ `allow_symbolic=false`, `allow_algebraic=true`, `allow_parameter=true`, and
+  `conservative=false` limit 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``. `allow_symbolic` allows arbitrary
+  symbolic coefficients. If it is false, `allow_algebraic` allows symbolic coefficients
+  involving only algebraic variables and parameters. Otherwise, `allow_parameter` only
+  allows coefficients containing parameters.
++ `fully_determined=true` controls whether or not an error will be thrown if the number
+  of equations don't match the number of inputs, outputs, and equations.
 """
 function structural_simplify(
         sys::AbstractSystem, io = nothing; additional_passes = [], simplify = false, split = true,
-        allow_symbolic = true, allow_parameter = true, conservative = false, fully_determined = true,
-        kwargs...)
+        allow_symbolic = false, allow_algebraic = true, allow_parameter = true, conservative = false,
+        fully_determined = true, kwargs...)
     isscheduled(sys) && throw(RepeatedStructuralSimplificationError())
-    newsys′ = __structural_simplify(sys, io; simplify,
-        allow_symbolic, allow_parameter, conservative, fully_determined,
-        kwargs...)
+    newsys′ = __structural_simplify(sys, io; simplify, allow_symbolic, allow_algebraic,
+        allow_parameter, conservative, fully_determined, kwargs...)
     if newsys′ isa Tuple
         @assert length(newsys′) == 2
         newsys = newsys′[1]

src/systems/systemstructure.jl

@@ -693,7 +693,7 @@ end
 
 function _structural_simplify!(state::TearingState, io; simplify = false,
         check_consistency = true, fully_determined = true, warn_initialize_determined = false,
-        dummy_derivative = true,
+        dummy_derivative = true, allow_symbolic = false,
         kwargs...)
     if fully_determined isa Bool
         check_consistency &= fully_determined
@@ -710,24 +710,27 @@ function _structural_simplify!(state::TearingState, io; simplify = false,
     else
         input_idxs = 0:-1 # Empty range
     end
-    sys, mm = ModelingToolkit.alias_elimination!(state; kwargs...)
+    sys, mm = ModelingToolkit.alias_elimination!(state; allow_symbolic, kwargs...)
     if check_consistency
         fully_determined = ModelingToolkit.check_consistency(
             state, orig_inputs; nothrow = fully_determined === nothing)
     end
     if fully_determined && dummy_derivative
         sys = ModelingToolkit.dummy_derivative(
-            sys, state; simplify, mm, check_consistency, kwargs...)
+            sys, state; simplify, mm, check_consistency, allow_symbolic, kwargs...)
     elseif fully_determined
-        var_eq_matching = pantelides!(state; finalize = false, kwargs...)
+        var_eq_matching = pantelides!(state; finalize = false, allow_symbolic, kwargs...)
+        if !allow_symbolic
+            StructuralTransformations.make_differential_denominators_unsolvable!(state.structure)
+        end
         sys = pantelides_reassemble(state, var_eq_matching)
         state = TearingState(sys)
-        sys, mm = ModelingToolkit.alias_elimination!(state; kwargs...)
+        sys, mm = ModelingToolkit.alias_elimination!(state; allow_symbolic, kwargs...)
         sys = ModelingToolkit.dummy_derivative(
-            sys, state; simplify, mm, check_consistency, kwargs...)
+            sys, state; simplify, mm, check_consistency, allow_symbolic, kwargs...)
     else
         sys = ModelingToolkit.tearing(
-            sys, state; simplify, mm, check_consistency, kwargs...)
+            sys, state; simplify, mm, check_consistency, allow_symbolic, kwargs...)
     end
     fullunknowns = [map(eq -> eq.lhs, observed(sys)); unknowns(sys)]
     @set! sys.observed = ModelingToolkit.topsort_equations(observed(sys), fullunknowns)