Merge branch 'master' into controlsystem

Author: YingboMa

Project.toml

@@ -79,9 +79,13 @@ Unitful = "1.1"
 julia = "1.6"
 
 [extras]
+AmplNLWriter = "7c4d4715-977e-5154-bfe0-e096adeac482"
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9"
+Ipopt_jll = "9cc047cb-c261-5740-88fc-0cf96f7bdcc7"
 Optimization = "7f7a1694-90dd-40f0-9382-eb1efda571ba"
+OptimizationMOI = "fd9f6733-72f4-499f-8506-86b2bdd0dea1"
 OptimizationOptimJL = "36348300-93cb-4f02-beb5-3c3902f8871e"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
@@ -93,4 +97,4 @@ Sundials = "c3572dad-4567-51f8-b174-8c6c989267f4"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["BenchmarkTools", "ForwardDiff", "Optimization", "OptimizationOptimJL", "OrdinaryDiffEq", "Random", "ReferenceTests", "SafeTestsets", "SteadyStateDiffEq", "Test", "StochasticDiffEq", "Sundials"]
+test = ["AmplNLWriter", "BenchmarkTools", "ForwardDiff", "Ipopt", "Ipopt_jll", "Optimization", "OptimizationOptimJL", "OptimizationMOI", "OrdinaryDiffEq", "Random", "ReferenceTests", "SafeTestsets", "SteadyStateDiffEq", "Test", "StochasticDiffEq", "Sundials"]

src/bipartite_graph.jl

@@ -184,6 +184,53 @@ function complete(g::BipartiteGraph{I}) where {I}
     BipartiteGraph(g.ne, g.fadjlist, badjlist)
 end
 
+# Matrix whose only purpose is to pretty-print the bipartite graph
+struct BipartiteAdjacencyList
+    u::Union{Vector{Int}, Nothing}
+end
+function Base.show(io::IO, l::BipartiteAdjacencyList)
+    if l.u === nothing
+        printstyled(io, '⋅', color = :light_black)
+    elseif isempty(l.u)
+        printstyled(io, '∅', color = :light_black)
+    else
+        print(io, l.u)
+    end
+end
+
+struct Label
+    s::String
+end
+Base.show(io::IO, l::Label) = print(io, l.s)
+
+struct BipartiteGraphPrintMatrix <:
+       AbstractMatrix{Union{Label, Int, BipartiteAdjacencyList}}
+    bpg::BipartiteGraph
+end
+Base.size(bgpm::BipartiteGraphPrintMatrix) = (max(nsrcs(bgpm.bpg), ndsts(bgpm.bpg)) + 1, 3)
+function Base.getindex(bgpm::BipartiteGraphPrintMatrix, i::Integer, j::Integer)
+    checkbounds(bgpm, i, j)
+    if i == 1
+        return (Label.(("#", "src", "dst")))[j]
+    elseif j == 1
+        return i - 1
+    elseif j == 2
+        return BipartiteAdjacencyList(i - 1 <= nsrcs(bgpm.bpg) ?
+                                      𝑠neighbors(bgpm.bpg, i - 1) : nothing)
+    elseif j == 3
+        return BipartiteAdjacencyList(i - 1 <= ndsts(bgpm.bpg) ?
+                                      𝑑neighbors(bgpm.bpg, i - 1) : nothing)
+    else
+        @assert false
+    end
+end
+
+function Base.show(io::IO, b::BipartiteGraph)
+    print(io, "BipartiteGraph with (", length(b.fadjlist), ", ",
+          isa(b.badjlist, Int) ? b.badjlist : length(b.badjlist), ") (𝑠,𝑑)-vertices\n")
+    Base.print_matrix(io, BipartiteGraphPrintMatrix(b))
+end
+
 """
 ```julia
 Base.isequal(bg1::BipartiteGraph{T}, bg2::BipartiteGraph{T}) where {T<:Integer}

src/structural_transformation/bipartite_tearing/modia_tearing.jl

@@ -37,17 +37,57 @@ end
 
 function tear_graph_modia(structure::SystemStructure; varfilter = v -> true,
                           eqfilter = eq -> true)
+    # It would be possible here to simply iterate over all variables and attempt to
+    # use tearEquations! to produce a matching that greedily selects the minimal
+    # number of torn variables. However, we can do this process faster if we first
+    # compute the strongly connected components. In the absence of cycles and
+    # non-solvability, a maximal matching on the original graph will give us an
+    # optimal assignment. However, even with cycles, we can use the maximal matching
+    # to give us a good starting point for a good matching and then proceed to
+    # reverse edges in each scc to improve the solution. Note that it is possible
+    # to have optimal solutions that cannot be found by this process. We will not
+    # find them here [TODO: It would be good to have an explicit example of this.]
+
     @unpack graph, solvable_graph = structure
     var_eq_matching = complete(maximal_matching(graph, eqfilter, varfilter))
     var_sccs::Vector{Union{Vector{Int}, Int}} = find_var_sccs(graph, var_eq_matching)
 
+    # Here, we're using a maximal matching on the post-pantelides system to find
+    # the strongly connected components of the system (of variables that depend
+    # on each other). The strongly connected components are unique, however, the
+    # maximal matching itself is not. Every maximal matching gives rise to the
+    # same set of strongly connected components, but the associated equations need
+    # not be the same. In the absence of solvability constraints, this may be a
+    # small issue, but here it is possible that an equation got assigned to an
+    # scc that cannot actually use it for solving a variable, but still precludes
+    # another scc from using it. To avoid this, we delete any assignments that
+    # are not in the solvable graph and extend the set of considered eqauations
+    # below.
+    for var in ndsts(solvable_graph)
+        var_eq_matching[var] === unassigned && continue
+        if !(BipartiteEdge(var, var_eq_matching[var]) in solvable_graph)
+            var_eq_matching[var] = unassigned
+        end
+    end
+
     for vars in var_sccs
         filtered_vars = filter(varfilter, vars)
         ieqs = Int[var_eq_matching[v]
                    for v in filtered_vars if var_eq_matching[v] !== unassigned]
         for var in vars
             var_eq_matching[var] = unassigned
         end
+        for var in filtered_vars
+            # Add any equations that we may not have been able to use earlier to see
+            # if a different matching may have been possible.
+            for eq′ in 𝑑neighbors(solvable_graph, var)
+                eqfilter(eq′) || continue
+                eq′ in ieqs && continue
+                if invview(var_eq_matching)[eq′] === unassigned
+                    push!(ieqs, eq′)
+                end
+            end
+        end
         tear_graph_block_modia!(var_eq_matching, graph, solvable_graph, ieqs, filtered_vars)
     end
 

src/systems/alias_elimination.jl

@@ -57,7 +57,13 @@ function alias_elimination(sys)
 
     newstates = []
     for j in eachindex(fullvars)
-        if !(j in keys(ag))
+        if j in keys(ag)
+            # Put back equations for alias eliminated dervars
+            if isdervar(state.structure, j) &&
+               !(invview(state.structure.var_to_diff)[j] in keys(ag))
+                push!(eqs, fullvars[j] ~ subs[fullvars[j]])
+            end
+        else
             isdervar(state.structure, j) || push!(newstates, fullvars[j])
         end
     end
@@ -192,6 +198,7 @@ struct AliasGraphKeySet <: AbstractSet{Int}
 end
 Base.keys(ag::AliasGraph) = AliasGraphKeySet(ag)
 Base.iterate(agk::AliasGraphKeySet, state...) = Base.iterate(agk.ag.eliminated, state...)
+Base.length(agk::AliasGraphKeySet) = Base.length(agk.ag.eliminated)
 function Base.in(i::Int, agk::AliasGraphKeySet)
     aliasto = agk.ag.aliasto
     1 <= i <= length(aliasto) && aliasto[i] !== nothing
@@ -210,9 +217,11 @@ function aag_bareiss!(graph, var_to_diff, mm_orig::SparseMatrixCLIL)
         is_linear_equations[e] = true
     end
 
-    # Variables that are highest order differentiated cannot be states of an ODE
-    is_not_potential_state = isnothing.(var_to_diff)
-    is_linear_variables = copy(is_not_potential_state)
+    # For now, only consider variables linear that are not differentiated.
+    # We could potentially apply the same logic to variables whose derivative
+    # is also linear, but that's a TODO.
+    diff_to_var = invview(var_to_diff)
+    is_linear_variables = .&(isnothing.(var_to_diff), isnothing.(diff_to_var))
     for i in 𝑠vertices(graph)
         is_linear_equations[i] && continue
         for j in 𝑠neighbors(graph, i)
@@ -230,11 +239,9 @@ function aag_bareiss!(graph, var_to_diff, mm_orig::SparseMatrixCLIL)
                 r !== nothing && return r
                 rank1 = k - 1
             end
-            if rank2 === nothing
-                r = find_masked_pivot(is_not_potential_state, M, k)
-                r !== nothing && return r
-                rank2 = k - 1
-            end
+            # TODO: It would be better to sort the variables by
+            # derivative order here to enable more elimination
+            # opportunities.
             return find_masked_pivot(nothing, M, k)
         end
         function find_and_record_pivot(M, k)
@@ -249,10 +256,9 @@ function aag_bareiss!(graph, var_to_diff, mm_orig::SparseMatrixCLIL)
         end
         bareiss_ops = ((M, i, j) -> nothing, myswaprows!,
                        bareiss_update_virtual_colswap_mtk!, bareiss_zero!)
-        rank3, = bareiss!(M, bareiss_ops; find_pivot = find_and_record_pivot)
-        rank1 = something(rank1, rank3)
-        rank2 = something(rank2, rank3)
-        (rank1, rank2, rank3, pivots)
+        rank2, = bareiss!(M, bareiss_ops; find_pivot = find_and_record_pivot)
+        rank1 = something(rank1, rank2)
+        (rank1, rank2, pivots)
     end
 
     return mm, solvable_variables, do_bareiss!(mm, mm_orig)
@@ -266,16 +272,27 @@ function alias_eliminate_graph!(graph, var_to_diff, mm_orig::SparseMatrixCLIL)
     # variables`.
     #
     # `do_bareiss` conceptually gives us this system:
-    # rank1 | [ M₁₁  M₁₂ | M₁₃  M₁₄ ]   [v₁] = [0]
-    # rank2 | [ 0    M₂₂ | M₂₃  M₂₄ ] P [v₂] = [0]
+    # rank1 | [ M₁₁  M₁₂ | M₁₃ ]   [v₁] = [0]
+    # rank2 | [ 0    M₂₂ | M₂₃ ] P [v₂] = [0]
     # -------------------|------------------------
-    # rank3 | [ 0    0   | M₃₃  M₃₄ ]   [v₃] = [0]
-    #         [ 0    0   | 0    0   ]   [v₄] = [0]
-    mm, solvable_variables, (rank1, rank2, rank3, pivots) = aag_bareiss!(graph, var_to_diff,
-                                                                         mm_orig)
+    #         [ 0    0   | 0   ]   [v₃] = [0]
+    #
+    # Where `v₁` are the purely linear variables (i.e. those that only appear in linear equations),
+    # `v₂` are the variables that may be potentially solved by the linear system and v₃ are the variables
+    # that contribute to the equations, but are not solved by the linear system. Note
+    # that the complete system may be larger than the linear subsystem and include variables
+    # that do not appear here.
+    mm, solvable_variables, (rank1, rank2, pivots) = aag_bareiss!(graph, var_to_diff,
+                                                                  mm_orig)
 
     # Step 2: Simplify the system using the Bareiss factorization
+
     ag = AliasGraph(size(mm, 2))
+
+    # First, eliminate variables that only appear in linear equations and were removed
+    # completely from the coefficient matrix. These are technically singularities in
+    # the matrix, but assigning them to 0 is a feasible assignment and works well in
+    # practice.
     for v in setdiff(solvable_variables, @view pivots[1:rank1])
         ag[v] = 0
     end
@@ -287,11 +304,7 @@ function alias_eliminate_graph!(graph, var_to_diff, mm_orig::SparseMatrixCLIL)
     function lss!(ei::Integer)
         vi = pivots[ei]
         may_eliminate = true
-        for v in 𝑠neighbors(graph, mm.nzrows[ei])
-            # the differentiated variable cannot be eliminated
-            may_eliminate &= isnothing(diff_to_var[v]) && isnothing(var_to_diff[v])
-        end
-        locally_structure_simplify!((@view mm[ei, :]), vi, ag, may_eliminate)
+        locally_structure_simplify!((@view mm[ei, :]), vi, ag, var_to_diff)
     end
 
     # Step 2.1: Go backwards, collecting eliminated variables and substituting
@@ -333,7 +346,7 @@ function exactdiv(a::Integer, b)
     return d
 end
 
-function locally_structure_simplify!(adj_row, pivot_col, ag, may_eliminate)
+function locally_structure_simplify!(adj_row, pivot_col, ag, var_to_diff)
     pivot_val = adj_row[pivot_col]
     iszero(pivot_val) && return false
 
@@ -375,21 +388,36 @@ function locally_structure_simplify!(adj_row, pivot_col, ag, may_eliminate)
         end
     end
 
-    if may_eliminate && nirreducible <= 1
+    if nirreducible <= 1
         # There were only one or two terms left in the equation (including the
         # pivot variable). We can eliminate the pivot variable.
         #
         # Note that when `nirreducible <= 1`, `alias_candidate` is uniquely
         # determined.
         if alias_candidate !== 0
+            # Verify that the derivative depth of the variable is at least
+            # as deep as that of the alias, otherwise, we can't eliminate.
+            pivot_var = pivot_col
+            alias_var = alias_candidate[2]
+            while (pivot_var = var_to_diff[pivot_col]) !== nothing
+                alias_var = var_to_diff[alias_var]
+                alias_var === nothing && return false
+            end
             d, r = divrem(alias_candidate[1], pivot_val)
             if r == 0 && (d == 1 || d == -1)
                 alias_candidate = -d => alias_candidate[2]
             else
                 return false
             end
         end
-        ag[pivot_col] = alias_candidate
+        diff_alias_candidate(ac) = ac === 0 ? 0 : ac[1] => var_to_diff[ac[2]]
+        while true
+            @assert !haskey(ag, pivot_col)
+            ag[pivot_col] = alias_candidate
+            pivot_col = var_to_diff[pivot_col]
+            pivot_col === nothing && break
+            alias_candidate = diff_alias_candidate(alias_candidate)
+        end
         zero!(adj_row)
         return true
     end