Pantelides algorithm

src/ModelingToolkit.jl

@@ -83,6 +83,7 @@ include("systems/diffeqs/odesystem.jl")
 include("systems/diffeqs/sdesystem.jl")
 include("systems/diffeqs/abstractodesystem.jl")
 include("systems/diffeqs/first_order_transform.jl")
+include("systems/diffeqs/index_reduction.jl")
 include("systems/diffeqs/modelingtoolkitize.jl")
 include("systems/diffeqs/validation.jl")
 

src/systems/diffeqs/first_order_transform.jl

@@ -19,8 +19,8 @@ Takes a Nth order ODESystem and returns a new ODESystem written in first order
 form by defining new variables which represent the N-1 derivatives.
 """
 function ode_order_lowering(sys::ODESystem)
-    eqs_lowered, _ = ode_order_lowering(sys.eqs, sys.iv)
-    ODESystem(eqs_lowered, sys.iv, [var_from_nested_derivative(eq.lhs)[1] for eq in eqs_lowered], sys.ps)
+    eqs_lowered, new_vars = ode_order_lowering(sys.eqs, sys.iv)
+    ODESystem(eqs_lowered, sys.iv, vcat(new_vars, states(sys)), sys.ps)
 end
 
 function ode_order_lowering(eqs, iv)
@@ -30,14 +30,18 @@ function ode_order_lowering(eqs, iv)
     new_vars = Variable[]
 
     for (i, eq) ∈ enumerate(eqs)
-        var, maxorder = var_from_nested_derivative(eq.lhs)
-        if maxorder > get(var_order, var, 0)
-            var_order[var] = maxorder
-            any(isequal(var), vars) || push!(vars, var)
+        if isequal(eq.lhs, Constant(0))
+            push!(new_eqs, eq)
+        else
+            var, maxorder = var_from_nested_derivative(eq.lhs)
+            if maxorder > get(var_order, var, 0)
+                var_order[var] = maxorder
+                any(isequal(var), vars) || push!(vars, var)
+            end
+            var′ = lower_varname(var, iv, maxorder - 1)
+            rhs′ = rename_lower_order(eq.rhs)
+            push!(new_eqs,Differential(iv())(var′(iv())) ~ rhs′)
         end
-        var′ = lower_varname(var, iv, maxorder - 1)
-        rhs′ = rename_lower_order(eq.rhs)
-        push!(new_eqs,Differential(iv())(var′(iv())) ~ rhs′)
     end
 
     for var ∈ vars

src/systems/diffeqs/index_reduction.jl

@@ -0,0 +1,162 @@
+# V-nodes `[x_1, x_2, x_3, ..., dx_1, dx_2, ..., y_1, y_2, ...]` where `x`s are
+# differential variables and `y`s are algebraic variables.
+function get_vnodes(sys)
+    dxvars = Operation[]
+    edges = map(_->Int[], 1:length(sys.eqs))
+    for (i, eq) in enumerate(sys.eqs)
+        if !(eq.lhs isa Constant)
+            # Make sure that the LHS is a first order derivative of a var.
+            @assert eq.lhs.op isa Differential
+            @assert !(eq.lhs.args[1] isa Differential) # first order
+
+            push!(dxvars, eq.lhs)
+            # For efficiency we note down the diff edges here
+            push!(edges[i], length(dxvars))
+        end
+    end
+
+    xvars = (first ∘ var_from_nested_derivative).(dxvars)
+    algvars  = setdiff(states(sys), xvars)
+    return xvars, dxvars, edges, algvars
+end
+
+function sys2bigraph(sys)
+    xvars, dxvars, edges, algvars = get_vnodes(sys)
+    xvar_offset = length(xvars)
+    algvar_offset = 2xvar_offset
+    for edge in edges
+        isempty(edge) || (edge .+= xvar_offset)
+    end
+
+    for (i, eq) in enumerate(sys.eqs)
+        # T or D(x):
+        # We assume no derivatives appear on the RHS at this point
+        vs = vars(eq.rhs)
+        for v in vs
+            for (j, target_v) in enumerate(xvars)
+                if v == target_v
+                    push!(edges[i], j)
+                end
+            end
+            for (j, target_v) in enumerate(algvars)
+                if v == target_v
+                    push!(edges[i], j+algvar_offset)
+                end
+            end
+        end
+    end
+
+    fullvars = [xvars; dxvars; algvars] # full list of variables
+    vars_asso = [(1:xvar_offset) .+ xvar_offset; zeros(Int, length(fullvars) - xvar_offset)] # variable association list
+    return edges, fullvars, vars_asso
+end
+
+print_bigraph(sys, vars, edges) = print_bigraph(stdout, sys, vars, edges)
+function print_bigraph(io::IO, sys, vars, edges)
+    println(io, "Equations:")
+    foreach(x->println(io, x), [i => sys.eqs[i] for i in 1:length(sys.eqs)])
+    for (i, edge) in enumerate(edges)
+        println(io, "\nEq $i has:")
+        print(io, '[')
+        for e in edge
+            print(io, "$(vars[e]), ")
+        end
+        print(io, ']')
+    end
+    return nothing
+end
+
+function match_equation!(edges, i, assign, active, vcolor=falses(length(active)), ecolor=falses(length(edges)))
+    # `edge[active]` are active edges
+    # i: equations
+    # j: variables
+    # assign: assign[j] == i means (i-j) is assigned
+    #
+    # color the equation
+    ecolor[i] = true
+    # if a V-node j exists s.t. edge (i-j) exists and assign[j] == 0
+    for j in edges[i]
+        if active[j] && assign[j] == 0
+            assign[j] = i
+            return true
+        end
+    end
+    # for every j such that edge (i-j) exists and j is uncolored
+    for j in edges[i]
+        (active[j] && !vcolor[j]) || continue
+        # color the variable
+        vcolor[j] = true
+        if match_equation!(edges, assign[j], assign, active, vcolor, ecolor)
+            assign[j] = i
+            return true
+        end
+    end
+    return false
+end
+
+function matching(edges, nvars, active=trues(nvars))
+    assign = zeros(Int, nvars)
+    for i in 1:length(edges)
+        match_equation!(edges, i, assign, active)
+    end
+    return assign
+end
+
+function pantelides(sys::ODESystem; kwargs...)
+    edges, fullvars, vars_asso = sys2bigraph(sys)
+    return pantelides!(edges, fullvars, vars_asso; kwargs...)
+end
+
+function pantelides!(edges, vars, vars_asso; maxiter = 8000)
+    neqs = length(edges)
+    nvars = length(vars)
+    assign = zeros(Int, nvars)
+    eqs_asso = fill(0, neqs)
+    neqs′ = neqs
+    for k in 1:neqs′
+        i = k
+        pathfound = false
+        # In practice, `maxiter=8000` should never be reached, otherwise, the
+        # index would be on the order of thousands.
+        for _ in 1:maxiter
+            # run matching on (dx, y) variables
+            active = vars_asso .== 0
+            vcolor = falses(nvars)
+            ecolor = falses(neqs)
+            pathfound = match_equation!(edges, i, assign, active, vcolor, ecolor)
+            pathfound && break # terminating condition
+            # for every colored V-node j
+            for j in eachindex(vcolor); vcolor[j] || continue
+                # introduce a new variable
+                nvars += 1
+                push!(vars_asso, 0)
+                vars_asso[j] = nvars
+                push!(assign, 0)
+            end
+
+            # for every colored E-node l
+            for l in eachindex(ecolor); ecolor[l] || continue
+                neqs += 1
+                # create new E-node
+                push!(edges, copy(edges[l]))
+                # create edges from E-node `neqs` to all V-nodes `j` and
+                # `vars_asso[j]` s.t. edge `(l-j)` exists
+                for j in edges[l]
+                    if !(vars_asso[j] in edges[neqs])
+                        push!(edges[neqs], vars_asso[j])
+                    end
+                end
+                push!(eqs_asso, 0)
+                eqs_asso[l] = neqs
+            end
+
+            # for every colored V-node j
+            for j in eachindex(vcolor); vcolor[j] || continue
+                assign[vars_asso[j]] = eqs_asso[assign[j]]
+            end
+            i = eqs_asso[i]
+        end # for _ in 1:maxiter
+        pathfound || error("maxiter=$maxiter reached! File a bug report if your system has a reasonable index (<100), and you are using the default `maxiter`. Try to increase the maxiter by `pantelides(sys::ODESystem; maxiter=1_000_000)` if your system has an incredibly high index and it is truly extremely large.")
+    end # for k in 1:neqs′
+    return edges, assign, vars_asso, eqs_asso
+end

src/systems/diffeqs/odesystem.jl

@@ -75,15 +75,17 @@ function ODESystem(deqs::AbstractVector{<:Equation}, iv, dvs, ps;
 end
 
 var_from_nested_derivative(x) = var_from_nested_derivative(x,0)
+var_from_nested_derivative(x::Constant) = (missing, missing)
 var_from_nested_derivative(x,i) = x.op isa Differential ? var_from_nested_derivative(x.args[1],i+1) : (x.op,i)
 iv_from_nested_derivative(x) = x.op isa Differential ? iv_from_nested_derivative(x.args[1]) : x.args[1].op
+iv_from_nested_derivative(x::Constant) = missing
 
 function ODESystem(eqs; kwargs...)
-    ivs = unique(iv_from_nested_derivative(eq.lhs) for eq ∈ eqs)
+    ivs = unique(skipmissing(iv_from_nested_derivative(eq.lhs) for eq ∈ eqs))
     length(ivs) == 1 || throw(ArgumentError("one independent variable currently supported"))
     iv = first(ivs)
 
-    dvs = unique(var_from_nested_derivative(eq.lhs)[1] for eq ∈ eqs)
+    dvs = unique(skipmissing(var_from_nested_derivative(eq.lhs)[1] for eq ∈ eqs))
     ps = filter(vars(eq.rhs for eq ∈ eqs)) do x
         isparameter(x) & !isequal(x, iv)
     end |> collect

src/utils.jl

@@ -67,6 +67,7 @@ Base.occursin(t::Expression, x::Expression) = isequal(x, t)
 vars(exprs) = foldl(vars!, exprs; init = Set{Variable}())
 function vars!(vars, O)
     isa(O, Operation) || return vars
+    O.op isa Variable && push!(vars, O.op)
     for arg ∈ O.args
         if isa(arg, Operation)
             isa(arg.op, Variable) && push!(vars, arg.op)

test/index_reduction.jl

@@ -0,0 +1,63 @@
+using ModelingToolkit
+using ModelingToolkit: sys2bigraph
+using DiffEqBase
+using Test
+
+# Define some variables
+@parameters t L g
+@variables x(t) y(t) w(t) z(t) T(t) xˍt(t) yˍt(t)
+@derivatives D'~t
+
+eqs2 = [D(D(x)) ~ T*x,
+        D(D(y)) ~ T*y - g,
+        0 ~ x^2 + y^2 - L^2]
+pendulum2 = ODESystem(eqs2, t, [x, y, T], [L, g], name=:pendulum)
+lowered_sys = ModelingToolkit.ode_order_lowering(pendulum2)
+
+lowered_eqs = [D(xˍt) ~ T*x,
+               D(yˍt) ~ T*y - g,
+               0 ~ x^2 + y^2 - L^2,
+               D(x) ~ xˍt,
+               D(y) ~ yˍt]
+@test ODESystem(lowered_eqs, t, [xˍt, yˍt, x, y, T], [L, g]) == lowered_sys
+@test isequal(lowered_sys.eqs, lowered_eqs)
+
+# Simple pendulum in cartesian coordinates
+eqs = [D(x) ~ w,
+       D(y) ~ z,
+       D(w) ~ T*x,
+       D(z) ~ T*y - g,
+       0 ~ x^2 + y^2 - L^2]
+pendulum = ODESystem(eqs, t, [x, y, w, z, T], [L, g], name=:pendulum)
+
+edges, vars, vars_asso = sys2bigraph(pendulum)
+@test ModelingToolkit.matching(edges, length(vars), vars_asso .== 0) == [0, 0, 0, 0, 1, 2, 3, 4, 0]
+
+edges, assign, vars_asso, eqs_asso = ModelingToolkit.pantelides(pendulum)
+
+@test sort.(edges) == sort.([
+ [5, 3],               # 1
+ [6, 4],               # 2
+ [7, 9, 1],            # 3
+ [8, 9, 2],            # 4
+ [2, 1],               # 5
+ [2, 1, 6, 5],         # 6
+ [5, 3, 10, 7],        # 7
+ [6, 4, 11, 8],        # 8
+ [2, 1, 6, 5, 11, 10], # 9
+])
+#                  [1, 2, 3, 4, 5,  6,  7,  8,  9, 10,   11]
+#                  [x, y, w, z, x', y', w', z', T, x'', y'']
+@test vars_asso == [5, 6, 7, 8, 10, 11, 0,  0,  0,  0,    0]
+#1: D(x) ~ w
+#2: D(y) ~ z
+#3: D(w) ~ T*x
+#4: D(z) ~ T*y - g
+#5: 0 ~ x^2 + y^2 - L^2
+# ----
+#6: D(5) -> 0 ~ 2xx'+ 2yy'
+#7: D(1) -> D(D(x)) ~ D(w)
+#8: D(2) -> D(D(y)) ~ D(z)
+#9: D(6) -> 0 ~ 2xx'' + 2x'x' + 2yy'' + 2y'y'
+#                 [1, 2, 3, 4, 5, 6, 7, 8, 9]
+@test eqs_asso == [7, 8, 0, 0, 6, 9, 0, 0, 0]