Merge pull request #2187 from SciML/aliasfix

Add back rank2 (highest diffed variables) in alias elimination

Author: YingboMa

examples/serial_inductor.jl

@@ -14,3 +14,20 @@ eqs = [connect(source.p, resistor.p)
 
 @named ll_model = ODESystem(eqs, t)
 ll_model = compose(ll_model, [source, resistor, inductor1, inductor2, ground])
+
+@named source = ConstantVoltage(V = 10.0)
+@named resistor1 = Resistor(R = 1.0)
+@named resistor2 = Resistor(R = 1.0)
+@named inductor1 = Inductor(L = 1.0e-2)
+@named inductor2 = Inductor(L = 2.0e-2)
+@named ground = Ground()
+
+eqs = [connect(source.p, inductor1.p)
+    connect(inductor1.n, resistor1.p)
+    connect(inductor1.n, resistor2.p)
+    connect(resistor1.n, resistor2.n)
+    connect(resistor2.n, inductor2.p)
+    connect(source.n, inductor2.n)
+    connect(inductor2.n, ground.g)]
+@named ll2_model = ODESystem(eqs, t)
+ll2_model = compose(ll2_model, [source, resistor1, resistor2, inductor1, inductor2, ground])

src/structural_transformation/StructuralTransformations.jl

@@ -54,6 +54,7 @@ export tearing_assignments, tearing_substitution
 export torn_system_jacobian_sparsity
 export full_equations
 export but_ordered_incidence, lowest_order_variable_mask, highest_order_variable_mask
+export computed_highest_diff_variables
 
 include("utils.jl")
 include("pantelides.jl")

src/structural_transformation/pantelides.jl

@@ -70,7 +70,7 @@ function pantelides_reassemble(state::TearingState, var_eq_matching)
 end
 
 """
-    computed_highest_diff_variables(var_to_diff)
+    computed_highest_diff_variables(structure)
 
 Computes which variables are the "highest-differentiated" for purposes of
 pantelides. Ordinarily this is relatively straightforward. However, in our

src/systems/alias_elimination.jl

@@ -264,7 +264,8 @@ function find_linear_variables(graph, linear_equations, var_to_diff, irreducible
     return linear_variables
 end
 
-function aag_bareiss!(graph, var_to_diff, mm_orig::SparseMatrixCLIL{T, Ti}) where {T, Ti}
+function aag_bareiss!(structure, mm_orig::SparseMatrixCLIL{T, Ti}) where {T, Ti}
+    @unpack graph, var_to_diff = structure
     mm = copy(mm_orig)
     linear_equations_set = BitSet(mm_orig.nzrows)
 
@@ -279,6 +280,7 @@ function aag_bareiss!(graph, var_to_diff, mm_orig::SparseMatrixCLIL{T, Ti}) wher
         v -> var_to_diff[v] === nothing === invview(var_to_diff)[v]
     end
     is_linear_variables = is_algebraic.(1:length(var_to_diff))
+    is_highest_diff = computed_highest_diff_variables(structure)
     for i in 𝑠vertices(graph)
         # only consider linear algebraic equations
         (i in linear_equations_set && all(is_algebraic, 𝑠neighbors(graph, i))) &&
@@ -291,25 +293,31 @@ function aag_bareiss!(graph, var_to_diff, mm_orig::SparseMatrixCLIL{T, Ti}) wher
 
     local bar
     try
-        bar = do_bareiss!(mm, mm_orig, is_linear_variables)
+        bar = do_bareiss!(mm, mm_orig, is_linear_variables, is_highest_diff)
     catch e
         e isa OverflowError || rethrow(e)
         mm = convert(SparseMatrixCLIL{BigInt, Ti}, mm_orig)
-        bar = do_bareiss!(mm, mm_orig, is_linear_variables)
+        bar = do_bareiss!(mm, mm_orig, is_linear_variables, is_highest_diff)
     end
 
     return mm, solvable_variables, bar
 end
 
-function do_bareiss!(M, Mold, is_linear_variables)
+function do_bareiss!(M, Mold, is_linear_variables, is_highest_diff)
     rank1r = Ref{Union{Nothing, Int}}(nothing)
+    rank2r = Ref{Union{Nothing, Int}}(nothing)
     find_pivot = let rank1r = rank1r
         (M, k) -> begin
             if rank1r[] === nothing
                 r = find_masked_pivot(is_linear_variables, M, k)
                 r !== nothing && return r
                 rank1r[] = k - 1
             end
+            if rank2r[] === nothing
+                r = find_masked_pivot(is_highest_diff, M, k)
+                r !== nothing && return r
+                rank2r[] = k - 1
+            end
             # TODO: It would be better to sort the variables by
             # derivative order here to enable more elimination
             # opportunities.
@@ -334,15 +342,19 @@ function do_bareiss!(M, Mold, is_linear_variables)
     bareiss_ops = ((M, i, j) -> nothing, myswaprows!,
         bareiss_update_virtual_colswap_mtk!, bareiss_zero!)
 
-    rank2, = bareiss!(M, bareiss_ops; find_pivot = find_and_record_pivot)
+    rank3, = bareiss!(M, bareiss_ops; find_pivot = find_and_record_pivot)
+    rank2 = something(rank2r[], rank3)
     rank1 = something(rank1r[], rank2)
-    (rank1, rank2, pivots)
+    (rank1, rank2, rank3, pivots)
 end
 
-function simple_aliases!(ils, graph, solvable_graph, eq_to_diff, var_to_diff)
-    ils, solvable_variables, (rank1, rank2, pivots) = aag_bareiss!(graph,
-        var_to_diff,
-        ils)
+function alias_eliminate_graph!(state::TransformationState, ils::SparseMatrixCLIL)
+    @unpack structure = state
+    @unpack graph, solvable_graph, var_to_diff, eq_to_diff = state.structure
+    # Step 1: Perform Bareiss factorization on the adjacency matrix of the linear
+    #         subsystem of the system we're interested in.
+    #
+    ils, solvable_variables, (rank1, rank2, rank3, pivots) = aag_bareiss!(structure, ils)
 
     ## Step 2: Simplify the system using the Bareiss factorization
     rk1vars = BitSet(@view pivots[1:rank1])
@@ -362,14 +374,6 @@ function simple_aliases!(ils, graph, solvable_graph, eq_to_diff, var_to_diff)
     return ils
 end
 
-function alias_eliminate_graph!(state::TransformationState, ils::SparseMatrixCLIL)
-    @unpack graph, solvable_graph, var_to_diff, eq_to_diff = state.structure
-    # Step 1: Perform Bareiss factorization on the adjacency matrix of the linear
-    #         subsystem of the system we're interested in.
-    #
-    return simple_aliases!(ils, graph, solvable_graph, eq_to_diff, var_to_diff)
-end
-
 function exactdiv(a::Integer, b)
     d, r = divrem(a, b)
     @assert r == 0

test/clock.jl

@@ -67,7 +67,7 @@ ci, varmap = infer_clocks(sys)
 eqmap = ci.eq_domain
 tss, inputs = ModelingToolkit.split_system(deepcopy(ci))
 sss, = SystemStructures._structural_simplify!(deepcopy(tss[1]), (inputs[1], ()))
-@test equations(sss) == [D(x) ~ u - y]
+@test equations(sss) == [D(x) ~ u - x]
 sss, = SystemStructures._structural_simplify!(deepcopy(tss[2]), (inputs[2], ()))
 @test isempty(equations(sss))
 @test observed(sss) == [yd ~ Sample(t, dt)(y); r ~ 1.0; ud ~ kp * (r - yd)]

test/components.jl

@@ -159,6 +159,12 @@ u0 = states(sys) .=> 0
 prob = DAEProblem(sys, Differential(t).(states(sys)) .=> 0, u0, (0, 0.5))
 @test_nowarn sol = solve(prob, DFBDF())
 
+sys2 = structural_simplify(ll2_model)
+@test length(equations(sys2)) == 3
+u0 = states(sys) .=> 0
+prob = ODEProblem(sys, u0, (0, 10.0))
+@test_nowarn sol = solve(prob, FBDF())
+
 @variables t x1(t) x2(t) x3(t) x4(t)
 D = Differential(t)
 @named sys1_inner = ODESystem([D(x1) ~ x1], t)

test/linearize.jl

@@ -207,12 +207,12 @@ if VERSION >= v"1.8"
 
     @named model = ODESystem(eqs, t, [], []; systems = [link1, cart, force, fixed])
     def = ModelingToolkit.defaults(model)
-    def[link1.y1] = 0
+    for s in states(model)
+        def[s] = 0
+    end
     def[link1.x1] = 10
+    def[link1.fy1] = -def[link1.g] * def[link1.m]
     def[link1.A] = -pi / 2
-    def[link1.dA] = 0
-    def[cart.s] = 0
-    def[force.flange.v] = 0
     lin_outputs = [cart.s, cart.v, link1.A, link1.dA]
     lin_inputs = [force.f.u]
 
@@ -237,7 +237,9 @@ if VERSION >= v"1.8"
     def = merge(def, Dict(x => 0.0 for x in dummyder))
 
     @test substitute(lsyss.A, def) ≈ lsys.A
-    @test substitute(lsyss.B, def) ≈ lsys.B
+    # We cannot pivot symbolically, so the part where a linear solve is required
+    # is not reliable.
+    @test substitute(lsyss.B, def)[1:6, :] ≈ lsys.B[1:6, :]
     @test substitute(lsyss.C, def) ≈ lsys.C
     @test substitute(lsyss.D, def) ≈ lsys.D
 end

test/reduction.jl

@@ -294,7 +294,7 @@ ss = structural_simplify(sys)
 @test isempty(equations(ss))
 @test sort(string.(observed(ss))) == ["x(t) ~ 0.0"
     "xˍt(t) ~ 0.0"
-    "y(t) ~ xˍt(t)"]
+    "y(t) ~ xˍt(t) - x(t)"]
 
 eqs = [D(D(x)) ~ -x]
 @named sys = ODESystem(eqs, t, [x], [])

test/structural_transformation/tearing.jl

@@ -138,11 +138,6 @@ eqs = [
     0 ~ x + z,
 ]
 @named nlsys = NonlinearSystem(eqs, [x, y, z], [])
-let (mm, _, _) = ModelingToolkit.aag_bareiss(nlsys)
-    @test mm == [-1 1 0;
-        0 -1 -1;
-        0 0 0]
-end
 
 newsys = tearing(nlsys)
 @test length(equations(newsys)) <= 1