Merge pull request #1738 from SciML/myb/alias_irr

Make sure that irreducible variables don't get eliminated by accident

Author: YingboMa

src/systems/alias_elimination.jl

@@ -570,7 +570,9 @@ function simple_aliases!(ag, graph, var_to_diff, mm_orig, irreducibles = ())
                                                                          irreducibles)
 
     # Step 2: Simplify the system using the Bareiss factorization
-    for v in setdiff(solvable_variables, @view pivots[1:rank1])
+    rk1vars = BitSet(@view pivots[1:rank1])
+    for v in solvable_variables
+        v in rk1vars && continue
         ag[v] = 0
     end
 
@@ -699,6 +701,16 @@ function alias_eliminate_graph!(graph, var_to_diff, mm_orig::SparseMatrixCLIL)
         end
         for (i, v) in enumerate(level_to_var)
             _alias = get(ag, v, nothing)
+            v_eqs = 𝑑neighbors(graph, v)
+            # if an irreducible appears in only one equation, we need to make
+            # sure that the other variables don't get eliminated
+            if length(v_eqs) == 1
+                eq = v_eqs[1]
+                for av in 𝑠neighbors(graph, eq)
+                    push!(irreducibles, av)
+                end
+                ag[v] = nothing
+            end
             push!(irreducibles, v)
             if _alias !== nothing && iszero(_alias[1]) && i < length(level_to_var)
                 # we have `x = 0`
@@ -709,7 +721,8 @@ function alias_eliminate_graph!(graph, var_to_diff, mm_orig::SparseMatrixCLIL)
         end
         if nlevels < (new_nlevels = length(level_to_var))
             for i in (nlevels + 1):new_nlevels
-                var_to_diff[level_to_var[i - 1]] = level_to_var[i]
+                li = level_to_var[i]
+                var_to_diff[level_to_var[i - 1]] = li
                 push!(updated_diff_vars, level_to_var[i - 1])
             end
         end

test/odesystem.jl

@@ -839,8 +839,10 @@ let
 
     sys_alias = alias_elimination(sys_con)
     D = Differential(t)
-    true_eqs = [0 ~ sys.v - D(sys.x)
-                0 ~ ctrl.kv * D(sys.x) + ctrl.kx * sys.x - D(sys.v)]
+    true_eqs = [0 ~ D(sys.v) - sys.u
+                0 ~ sys.x - ctrl.x
+                0 ~ sys.v - D(sys.x)
+                0 ~ ctrl.kv * D(sys.x) + ctrl.kx * ctrl.x - D(sys.v)]
     @test isequal(full_equations(sys_alias), true_eqs)
 
     sys_simp = structural_simplify(sys_con)

test/reduction.jl

@@ -244,3 +244,14 @@ eqs = [D(x) ~ σ * (y - x)
 lorenz1 = ODESystem(eqs, t, name = :lorenz1)
 lorenz1_reduced = structural_simplify(lorenz1)
 @test z in Set(parameters(lorenz1_reduced))
+
+# MWE for #1722
+@variables t
+vars = @variables a(t) w(t) phi(t)
+eqs = [a ~ D(w)
+       w ~ D(phi)
+       w ~ sin(t)]
+@named sys = ODESystem(eqs, t, vars, [])
+ss = alias_elimination(sys)
+@test equations(ss) == [0 ~ D(D(phi)) - a, 0 ~ sin(t) - D(phi)]
+@test observed(ss) == [w ~ D(phi)]