Merge pull request #1804 from Keno/kf/pantelidesrobust

Make pantelides robust against unused derivative vars

Author: YingboMa

src/inputoutput.jl

@@ -269,8 +269,8 @@ function inputs_to_parameters!(state::TransformationState, io)
         @assert new_v > 0
         new_var_to_diff[new_i] = new_v
     end
-    @set! structure.var_to_diff = new_var_to_diff
-    @set! structure.graph = new_graph
+    @set! structure.var_to_diff = complete(new_var_to_diff)
+    @set! structure.graph = complete(new_graph)
 
     @set! sys.eqs = map(Base.Fix2(substitute, input_to_parameters), equations(sys))
     @set! sys.states = setdiff(states(sys), keys(input_to_parameters))

src/structural_transformation/StructuralTransformations.jl

@@ -21,7 +21,8 @@ using ModelingToolkit: ODESystem, AbstractSystem, var_from_nested_derivative, Di
                        ExtraVariablesSystemException,
                        get_postprocess_fbody, vars!,
                        IncrementalCycleTracker, add_edge_checked!, topological_sort,
-                       invalidate_cache!, Substitutions, get_or_construct_tearing_state
+                       invalidate_cache!, Substitutions, get_or_construct_tearing_state,
+                       AliasGraph
 
 using ModelingToolkit.BipartiteGraphs
 import .BipartiteGraphs: invview

src/structural_transformation/pantelides.jl

@@ -74,14 +74,37 @@ end
 
 Perform Pantelides algorithm.
 """
-function pantelides!(state::TransformationState; maxiters = 8000)
+function pantelides!(state::TransformationState, ag::Union{AliasGraph, Nothing} = nothing;
+                     maxiters = 8000)
     @unpack graph, solvable_graph, var_to_diff, eq_to_diff = state.structure
     neqs = nsrcs(graph)
     nvars = nv(var_to_diff)
     vcolor = falses(nvars)
     ecolor = falses(neqs)
     var_eq_matching = Matching(nvars)
     neqs′ = neqs
+    nnonemptyeqs = count(eq -> !isempty(𝑠neighbors(graph, eq)), 1:neqs′)
+
+    # Allow matching for the highest differentiated variable that
+    # currently appears in an equation (or implicit equation in a side ag)
+    varwhitelist = falses(nvars)
+    for var in 1:nvars
+        if var_to_diff[var] === nothing
+            while isempty(𝑑neighbors(graph, var)) && (ag === nothing || !haskey(ag, var))
+                var′ = invview(var_to_diff)[var]
+                var′ === nothing && break
+                var = var′
+            end
+            if !isempty(𝑑neighbors(graph, var)) || (ag !== nothing && haskey(ag, var))
+                varwhitelist[var] = true
+            end
+        end
+    end
+
+    if nnonemptyeqs > count(varwhitelist)
+        throw(InvalidSystemException("System is structurally singular"))
+    end
+
     for k in 1:neqs′
         eq′ = k
         isempty(𝑠neighbors(graph, eq′)) && continue
@@ -93,7 +116,6 @@ function pantelides!(state::TransformationState; maxiters = 8000)
             #
             # the derivatives and algebraic variables are zeros in the variable
             # association list
-            varwhitelist = var_to_diff .== nothing
             resize!(vcolor, nvars)
             fill!(vcolor, false)
             resize!(ecolor, neqs)
@@ -103,11 +125,16 @@ function pantelides!(state::TransformationState; maxiters = 8000)
             pathfound && break # terminating condition
             for var in eachindex(vcolor)
                 vcolor[var] || continue
-                # introduce a new variable
-                nvars += 1
-                var_diff = var_derivative!(state, var)
-                push!(var_eq_matching, unassigned)
-                @assert length(var_eq_matching) == var_diff
+                if var_to_diff[var] === nothing
+                    # introduce a new variable
+                    nvars += 1
+                    var_diff = var_derivative!(state, var)
+                    push!(var_eq_matching, unassigned)
+                    push!(varwhitelist, false)
+                    @assert length(var_eq_matching) == var_diff
+                end
+                varwhitelist[var] = false
+                varwhitelist[var_to_diff[var]] = true
             end
 
             for eq in eachindex(ecolor)

src/structural_transformation/partial_state_selection.jl

@@ -1,6 +1,7 @@
-function partial_state_selection_graph!(state::TransformationState)
+function partial_state_selection_graph!(state::TransformationState;
+                                        ag::Union{AliasGraph, Nothing} = nothing)
     find_solvables!(state; allow_symbolic = true)
-    var_eq_matching = complete(pantelides!(state))
+    var_eq_matching = complete(pantelides!(state, ag))
     complete!(state.structure)
     partial_state_selection_graph!(state.structure, var_eq_matching)
 end

src/structural_transformation/symbolics_tearing.jl

@@ -567,7 +567,7 @@ function tearing_reassemble(state::TearingState, var_eq_matching; simplify = fal
     eq_to_diff = new_eq_to_diff
     diff_to_var = invview(var_to_diff)
 
-    @set! state.structure.graph = graph
+    @set! state.structure.graph = complete(graph)
     @set! state.structure.var_to_diff = var_to_diff
     @set! state.structure.eq_to_diff = eq_to_diff
     @set! state.fullvars = fullvars = fullvars[invvarsperm]

src/systems/systemstructure.jl

@@ -350,7 +350,8 @@ function TearingState(sys; quick_cancel = false, check = true)
     eq_to_diff = DiffGraph(nsrcs(graph))
 
     return TearingState(sys, fullvars,
-                        SystemStructure(var_to_diff, eq_to_diff, graph, nothing), Any[])
+                        SystemStructure(complete(var_to_diff), complete(eq_to_diff),
+                                        complete(graph), nothing), Any[])
 end
 
 function linear_subsys_adjmat(state::TransformationState)

test/structural_transformation/utils.jl

@@ -21,7 +21,7 @@ StructuralTransformations.find_solvables!(state)
 sss = state.structure
 @unpack graph, solvable_graph, var_to_diff = sss
 @test graph.fadjlist == [[1, 7], [2, 8], [3, 5, 9], [4, 6, 9], [5, 6]]
-@test graph.badjlist == 9
+@test length(graph.badjlist) == 9
 @test ne(graph) == nnz(incidence_matrix(graph)) == 12
 @test nv(solvable_graph) == 9 + 5
 let N = nothing