Add dummy_derivative (WIP)

Author: YingboMa

src/bipartite_graph.jl

@@ -234,13 +234,13 @@ Base.in(edge::BipartiteEdge, g::BipartiteGraph) = Graphs.has_edge(g, edge)
 
 ### Maximal matching
 """
-    construct_augmenting_path!(m::Matching, g::BipartiteGraph, vsrc, dstfilter, vcolor=falses(ndsts(g)), ecolor=falses(nsrcs(g))) -> path_found::Bool
+    construct_augmenting_path!(m::Matching, g::BipartiteGraph, vsrc, dstfilter, vcolor=falses(ndsts(g)), ecolor=nothing) -> path_found::Bool
 
 Try to construct an augmenting path in matching and if such a path is found,
 update the matching accordingly.
 """
-function construct_augmenting_path!(matching::Matching, g::BipartiteGraph, vsrc, dstfilter, dcolor=falses(ndsts(g)), scolor=falses(nsrcs(g)))
-    scolor[vsrc] = true
+function construct_augmenting_path!(matching::Matching, g::BipartiteGraph, vsrc, dstfilter, dcolor=falses(ndsts(g)), scolor=nothing)
+    scolor === nothing || (scolor[vsrc] = true)
 
     # if a `vdst` is unassigned and the edge `vsrc <=> vdst` exists
     for vdst in 𝑠neighbors(g, vsrc)

src/structural_transformation/StructuralTransformations.jl

@@ -27,7 +27,7 @@ using ModelingToolkit.BipartiteGraphs
 import .BipartiteGraphs: invview
 using Graphs
 using ModelingToolkit.SystemStructures
-using ModelingToolkit.SystemStructures: algeqs
+using ModelingToolkit.SystemStructures: algeqs, EquationsView
 
 using ModelingToolkit.DiffEqBase
 using ModelingToolkit.StaticArrays
@@ -40,6 +40,7 @@ using SparseArrays
 using NonlinearSolve
 
 export tearing, partial_state_selection, dae_index_lowering, check_consistency
+export dummy_derivative
 export build_torn_function, build_observed_function, ODAEProblem
 export sorted_incidence_matrix, pantelides!, tearing_reassemble, find_solvables!
 export tearing_assignments, tearing_substitution

src/structural_transformation/partial_state_selection.jl

@@ -139,3 +139,80 @@ function partial_state_selection_graph!(structure::SystemStructure, var_eq_match
 
     var_eq_matching
 end
+
+function dummy_derivative_graph!(state::TransformationState)
+    var_eq_matching = complete(pantelides!(state))
+    complete!(state.structure)
+    dummy_derivative_graph!(state.structure, var_eq_matching)
+end
+
+function dummy_derivative_graph!(structure::SystemStructure, var_eq_matching)
+    @unpack eq_to_diff, var_to_diff, graph = structure
+    diff_to_eq = invview(eq_to_diff)
+    diff_to_var = invview(var_to_diff)
+    invgraph = invview(graph)
+
+    neqs = nsrcs(graph)
+    eqlevel = zeros(Int, neqs)
+    maxlevel = 0
+    for i in 1:neqs
+        level = 0
+        eq = i
+        while diff_to_eq[eq] !== nothing
+            eq = diff_to_eq[eq]
+            level += 1
+        end
+        maxlevel = max(maxlevel, level)
+        eqlevel[i] = level
+    end
+
+    nvars = ndsts(graph)
+    varlevel = zeros(Int, nvars)
+    for i in 1:nvars
+        level = 0
+        var = i
+        while diff_to_var[var] !== nothing
+            var = diff_to_var[var]
+            level += 1
+        end
+        maxlevel = max(maxlevel, level)
+        varlevel[i] = level
+    end
+
+    var_sccs = find_var_sccs(graph, var_eq_matching)
+    eqcolor = falses(neqs)
+    dummy_derivatives = Int[]
+    for vars in var_sccs
+        eqs = [var_eq_matching[var] for var in vars if var_eq_matching[var] !== unassigned]
+        isempty(eqs) && continue
+        maxlevel = maximum(map(x->eqlevel[x], eqs))
+        iszero(maxlevel) && continue
+
+        rank_matching = Matching(nvars)
+        for level in maxlevel:-1:1
+            eqs = filter(eq->diff_to_eq[eq] !== nothing, eqs)
+            nrows = length(eqs)
+            iszero(nrows) && break
+            eqs_set = BitSet(eqs)
+
+            structural_rank = 0
+            for var in vars
+                pathfound = construct_augmenting_path!(rank_matching, invgraph, var, eq->eq in eqs_set, eqcolor)
+                pathfound || continue
+                push!(dummy_derivatives, var)
+                structural_rank += 1
+                structural_rank == nrows && break
+            end
+            if structural_rank != nrows
+                @warn "The DAE system is structurally singular!"
+            end
+            fill!(rank_matching, unassigned)
+
+            # prepare the next iteration
+            eqs = map(eq->diff_to_eq[eq], eqs)
+            vars = [diff_to_var[var] for var in vars if diff_to_var[var] !== nothing]
+        end
+    end
+
+    dummy_derivatives
+end

src/structural_transformation/symbolics_tearing.jl

@@ -282,3 +282,15 @@ function partial_state_selection(sys; simplify=false)
 
     tearing_reassemble(state, var_eq_matching; simplify=simplify)
 end
+
+"""
+    dummy_derivative(sys)
+
+Perform index reduction and use the dummy derivative techinque to ensure that
+the system is balanced.
+"""
+function dummy_derivative(sys)
+    state = TearingState(sys)
+    dds = dummy_derivative_graph!(state)
+    EquationsView(state), state.fullvars[dds]
+end