Add dae_order_lowering

Author: YingboMa

src/ModelingToolkit.jl

@@ -175,7 +175,7 @@ export ControlSystem
 export alias_elimination, flatten
 export connect, @connector, Connection, Flow, Stream, instream
 export isinput, isoutput, getbounds, hasbounds, isdisturbance, istunable, getdist, hasdist, tunable_parameters
-export ode_order_lowering, liouville_transform
+export ode_order_lowering, dae_order_lowering, liouville_transform
 export runge_kutta_discretize
 export PDESystem
 export Differential, expand_derivatives, @derivatives

src/structural_transformation/bareiss.jl

@@ -222,7 +222,7 @@ function nullspace(A; col_order=nothing)
             @swap(col_order[i],col_order[cp])
         end
     end
-    
+
     N = ModelingToolkit.reduced_echelon_nullspace(rank, B)
     apply_inv_pivot_rows!(N, column_pivots)
 end

src/systems/diffeqs/first_order_transform.jl

@@ -12,6 +12,14 @@ function ode_order_lowering(sys::ODESystem)
     return sys
 end
 
+function dae_order_lowering(sys::ODESystem)
+    iv = get_iv(sys)
+    eqs_lowered, new_vars = dae_order_lowering(equations(sys), iv, states(sys))
+    @set! sys.eqs = eqs_lowered
+    @set! sys.states = new_vars
+    return sys
+end
+
 function ode_order_lowering(eqs, iv, states)
     var_order = OrderedDict{Any,Int}()
     D = Differential(iv)
@@ -47,3 +55,51 @@ function ode_order_lowering(eqs, iv, states)
     # we want to order the equations and variables to be `(diff, alge)`
     return (vcat(diff_eqs, alge_eqs), vcat(diff_vars, setdiff(states, diff_vars)))
 end
+
+function dae_order_lowering(eqs, iv, states)
+    var_order = OrderedDict{Any,Int}()
+    D = Differential(iv)
+    diff_eqs = Equation[]
+    diff_vars = OrderedSet()
+    alge_eqs = Equation[]
+    vars = Set()
+    subs = Dict()
+
+    for (i, eq) ∈ enumerate(eqs)
+        vars!(vars, eq)
+        n_diffvars = 0
+        for vv in vars
+            isdifferential(vv) || continue
+            var, maxorder = var_from_nested_derivative(vv)
+            isparameter(var) && continue
+            n_diffvars += 1
+            order = get(var_order, var, nothing)
+            seen = order !== nothing
+            if !seen
+                order = 1
+            end
+            maxorder > order && (var_order[var] = maxorder)
+            var′ = lower_varname(var, iv, maxorder - 1)
+            subs[vv] = D(var′)
+            if !seen
+                push!(diff_vars, var′)
+            end
+        end
+        n_diffvars == 0 && push!(alge_eqs, eq)
+        empty!(vars)
+    end
+
+    for (var, order) ∈ var_order
+        for o in (order-1):-1:1
+            lvar = lower_varname(var, iv, o-1)
+            rvar = lower_varname(var, iv, o)
+            push!(diff_vars, lvar)
+
+            rhs = rvar
+            eq = Differential(iv)(lvar) ~ rhs
+            push!(diff_eqs, eq)
+        end
+    end
+
+    return ([diff_eqs; substitute.(eqs, (subs,))], vcat(collect(diff_vars), setdiff(states, diff_vars)))
+end