feat: enable structural_simplify to automatically handle non-fully-determined systems

Uses the result of `check_consistency`

Author: AayushSabharwal

src/structural_transformation/utils.jl

@@ -58,7 +58,41 @@ end
 ###
 ### Structural check
 ###
-function check_consistency(state::TransformationState, orig_inputs)
+
+"""
+    $(TYPEDSIGNATURES)
+
+Check if the `state` represents a singular system, and return the unmatched variables.
+"""
+function singular_check(state::TransformationState)
+    @unpack graph, var_to_diff = state.structure
+    fullvars = get_fullvars(state)
+    # This is defined to check if Pantelides algorithm terminates. For more
+    # details, check the equation (15) of the original paper.
+    extended_graph = (@set graph.fadjlist = Vector{Int}[graph.fadjlist;
+                                                        map(collect, edges(var_to_diff))])
+    extended_var_eq_matching = maximal_matching(extended_graph)
+
+    nvars = ndsts(graph)
+    unassigned_var = []
+    for (vj, eq) in enumerate(extended_var_eq_matching)
+        vj > nvars && break
+        if eq === unassigned && !isempty(𝑑neighbors(graph, vj))
+            push!(unassigned_var, fullvars[vj])
+        end
+    end
+    return unassigned_var
+end
+
+"""
+    $(TYPEDSIGNATURES)
+
+Check the consistency of `state`, given the inputs `orig_inputs`. If `nothrow == false`,
+throws an error if the system is under-/over-determined or singular. In this case, if the
+function returns it will return `true`. If `nothrow == true`, it will return `false`
+instead of throwing an error. The singular case will print a warning.
+"""
+function check_consistency(state::TransformationState, orig_inputs; nothrow = false)
     fullvars = get_fullvars(state)
     neqs = n_concrete_eqs(state)
     @unpack graph, var_to_diff = state.structure
@@ -72,6 +106,7 @@ function check_consistency(state::TransformationState, orig_inputs)
     is_balanced = n_highest_vars == neqs
 
     if neqs > 0 && !is_balanced
+        nothrow && return false
         varwhitelist = var_to_diff .== nothing
         var_eq_matching = maximal_matching(graph, eq -> true, v -> varwhitelist[v]) # not assigned
         # Just use `error_reporting` to do conditional
@@ -85,20 +120,7 @@ function check_consistency(state::TransformationState, orig_inputs)
         error_reporting(state, bad_idxs, n_highest_vars, iseqs, orig_inputs)
     end
 
-    # This is defined to check if Pantelides algorithm terminates. For more
-    # details, check the equation (15) of the original paper.
-    extended_graph = (@set graph.fadjlist = Vector{Int}[graph.fadjlist;
-                                                        map(collect, edges(var_to_diff))])
-    extended_var_eq_matching = maximal_matching(extended_graph)
-
-    nvars = ndsts(graph)
-    unassigned_var = []
-    for (vj, eq) in enumerate(extended_var_eq_matching)
-        vj > nvars && break
-        if eq === unassigned && !isempty(𝑑neighbors(graph, vj))
-            push!(unassigned_var, fullvars[vj])
-        end
-    end
+    unassigned_var = singular_check(state)
 
     if !isempty(unassigned_var) || !is_balanced
         io = IOBuffer()
@@ -107,10 +129,14 @@ function check_consistency(state::TransformationState, orig_inputs)
         errmsg = "The system is structurally singular! " *
                  "Here are the problematic variables: \n" *
                  unassigned_var_str
+        if nothrow
+            @warn errmsg
+            return false
+        end
         throw(InvalidSystemException(errmsg))
     end
 
-    return nothing
+    return true
 end
 
 ###

src/systems/systemstructure.jl

@@ -677,7 +677,11 @@ function _structural_simplify!(state::TearingState, io; simplify = false,
         check_consistency = true, fully_determined = true, warn_initialize_determined = false,
         dummy_derivative = true,
         kwargs...)
-    check_consistency &= fully_determined
+    if fully_determined isa Bool
+        check_consistency &= fully_determined
+    else
+        check_consistency = true
+    end
     has_io = io !== nothing
     orig_inputs = Set()
     if has_io
@@ -690,7 +694,8 @@ function _structural_simplify!(state::TearingState, io; simplify = false,
     end
     sys, mm = ModelingToolkit.alias_elimination!(state; kwargs...)
     if check_consistency
-        ModelingToolkit.check_consistency(state, orig_inputs)
+        fully_determined = ModelingToolkit.check_consistency(
+            state, orig_inputs; nothrow = fully_determined === nothing)
     end
     if fully_determined && dummy_derivative
         sys = ModelingToolkit.dummy_derivative(