Numerical and Structural Characterizations for DAEs

Numerical solvers cannot solve all DAEs. Consider the DAE system

\[F(\{x', y', z'\}, \{x, y, z\}, t) = \begin{pmatrix} x + y - \sin(t) \\ z - \sin(t) \\ z' - \cos(t) \end{pmatrix} = 0.\]

Note that the second equation and the third equation are equivalent, and there are not enough constraints to uniquely determine $x$ and $y$. Let's see how numerical solvers and ModelingToolkit behave.

using DifferentialEquations, Sundials, ModelingToolkit, Plots, LinearAlgebra
using ModelingToolkit: t_nounits as t, D_nounits as D

function f!(out, du, u, p, t)
    # u[1]: x, du[1]: x'
    # u[2]: y, du[2]: y'
    # u[3]: z, du[3]: z'
    out[1] = u[1] + u[2] - sin(t)
    out[2] = u[3] - sin(t)
    out[3] = du[3] - cos(t)
end
prob = DAEProblem(f!, [0.0, 0.0, 1.0], [0.0, 0.0, 0.0], (0, 100.0), differential_vars=[false, false, true])
sol1 = solve(prob, IDA())
sol2 = solve(prob, DFBDF())
println(sol1.retcode, "\t", sol2.retcode)

[IDAS ERROR]  IDASolve
  At t = 0 and h = 8.25906e-13, the corrector convergence failed repeatedly or with |h| = hmin.

ConvergenceFailure	Success

As expected, numerical solvers exit early and they cannot take more than one step. We call such systems non-integrable or singular. Let's try to use ModelingToolkit to analyze this problem.

@variables x(t) y(t) z(t)

eqs = [
    x + y ~ sin(t)
    z ~ sin(t)
    D(z) ~ cos(t)
]
@named sys = ODESystem(eqs, t)
sys = complete(sys);
try structural_simplify(sys) catch e println(e) end

ModelingToolkit.InvalidSystemException("The system is structurally singular! Here are the problematic variables: \n y(t)")

We can see that ModelingToolkit also errors, correctly identifying that the system is singular.

Let's consider another DAE system

\[F(\{x', y'\}, \{x, y\}, t) = \begin{pmatrix} x - \sin(t) \\ x' - y \end{pmatrix} = 0.\]

If we plug the first equation into the second equation, we have

\[y(t) = (\sin(t))' = \cos(t).\]

Superficially, we can have the initial condition

\[x'(0) = x(0) = y(0) = y'(0) = 0.\]

However, if we differentiate the first equation once, we get the hidden constraint

\[x'(t) - \cos(t) = 0.\]

Thus, the true consistent initial condition is

\[x'(0) = y(0) = 1, x(0) = y'(0) = 0.\]

Let's try to solve this simple DAE using a numerical solver.

function f!(out, du, u, p, t)
    # u[1]: x, du[1]: x'
    # u[2]: y, du[2]: y'
    out[1] = u[1] - sin(t)
    out[2] = du[1] - u[2]
end
prob = DAEProblem(f!, [1, 0.0], [0.0, 1.0], (0, 100pi), differential_vars=[true, false])
sol1 = solve(prob, IDA())
sol2 = solve(prob, DFBDF())
println("[sol1: ", sol1.retcode, "]",
"\n", "[sol2 ", sol2.retcode, ": y(100pi)=", sol2[2, end], " steps: ", length(sol2.t), "]")

[IDAS ERROR]  IDASolve
  At t = 16.335 and h = 8.49321e-06, the error test failed repeatedly or with |h| = hmin.

[sol1: Unstable]
[sol2 Success: y(100pi)=1.0001207192054915 steps: 3086]

To better understand the numerical behavior, let's analyze the variable step size behavior of the implicit Euler method of the original system. The local truncation error for $y$ is

\[\begin{align} &\text{lte} = \frac{\frac{y(t_n) - y(t_{n-1})}{h_n} - \frac{y(t_{n-1}) - y(t_{n-2})}{h_{n-1}}} {h_{n} + h_{n-1}} (h_{n} + h_{n-1}) h_n \\ =& \left(\frac{y(t_n) - y(t_{n-1})}{h_n} - \frac{y(t_{n-1}) - y(t_{n-2})}{h_{n-1}}\right) h_n \\ =& h_n\left(\frac{\frac{\sin(t_{n}) - \sin(t_{n-1})}{h_n} - \frac{\sin(t_{n-1}) - \sin(t_{n-2})}{h_{n-1}}}{h_n} - \frac{\frac{\sin(t_{n-1}) - \sin(t_{n-2})}{h_{n-1}} - \frac{\sin(t_{n-2}) - \sin(t_{n-3})}{h_{n-2}}}{h_{n-1}}\right) \\ =& \frac{\sin(t_{n}) - \sin(t_{n-1})}{h_n} - \frac{\sin(t_{n-1}) - \sin(t_{n-2})}{h_{n-1}} - h_n \frac{\frac{\sin(t_{n-1}) - \sin(t_{n-2})}{h_{n-1}} - \frac{\sin(t_{n-2}) - \sin(t_{n-3})}{h_{n-2}}}{h_{n-1}} \end{align}\]

Note that when $h_{n} \to 0$, the local truncation error becomes

\[\lim_{h_{n}\to 0} \text{lte} = \frac{\sin(t_{n}) - \sin(t_{n-1})}{h_n} - \frac{\sin(t_{n-1}) - \sin(t_{n-2})}{h_{n-1}} = \cos(t_{n}) - \frac{\sin(t_{n-1}) - \sin(t_{n-2})}{h_{n-1}}\]

Note that this in general does not converge to $0$. Thus, numerical solvers could have difficulties in solving this system. Let's also confirm this numerical behavior by setting the maximum order to $1$.

solve(prob, DFBDF(max_order=Val(1)))

retcode: MaxIters
Interpolation: 1st order linear
t: 999997-element Vector{Float64}:
   0.0
   5.9122259226118514e-8
   3.482873403416318e-7
   9.878837822595348e-7
   1.939419135851393e-6
   3.1005754683177864e-6
   4.38401403708843e-6
   5.734204194967433e-6
   7.119991498544613e-6
   8.52489468131624e-6
   ⋮
 168.01102271497368
 168.0111152896093
 168.01120773239938
 168.0113000435422
 168.01139222321862
 168.01148427162423
 168.01157618895002
 168.01166797537934
 168.01175963109773
u: 999997-element Vector{Vector{Float64}}:
 [0.0, 1.0]
 [5.912225922611848e-8, 0.9999999999999996]
 [3.4828734034162474e-7, 0.9999999999999758]
 [9.878837822593741e-7, 0.9999999999997597]
 [1.939419135850177e-6, 0.9999999999988911]
 [3.1005754683128185e-6, 0.999999999996769]
 [4.384014037074387e-6, 0.9999999999929289]
 [5.7342041949360085e-6, 0.9999999999871271]
 [7.119991498484456e-6, 0.9999999999792661]
 [8.524894681212984e-6, 0.9999999999693221]
 ⋮
 [-0.9979408979309449, -0.06418644985168312]
 [-0.9979468314096368, -0.06409399995040455]
 [-0.9979527479036147, -0.06400168118153492]
 [-0.9979586474619813, -0.06390949332473912]
 [-0.9979645301326023, -0.06381743622441657]
 [-0.9979703959641649, -0.06372550968654046]
 [-0.9979762450049405, -0.0636337135112199]
 [-0.9979820773025937, -0.06354204751854926]
 [-0.99798789290481, -0.06345051151707824]

Let's replace the first equation with the differentiated equation and solve it numerically,

function g!(out, du, u, p, t)
    # u[1]: x, du[1]: x'
    # u[2]: y, du[2]: y'
    out[1] = du[1] - cos(t)
    out[2] = du[1] - u[2]
end
prob = DAEProblem(g!, [1, 0.0], [0.0, 1.0], (0, 100pi), differential_vars=[true, false])
sol1_diff = solve(prob, IDA())
sol2_diff = solve(prob, DFBDF())
println("[sol1_diff: ", sol1_diff.retcode, ": y(100pi)=", sol1_diff[2, end], " steps: ", length(sol1_diff.t), "]",
"\n", "[sol2_diff ", sol2_diff.retcode, ": y(100pi)=", sol2_diff[2, end], " steps: ", length(sol2_diff.t), "]")

[sol1_diff: Success: y(100pi)=1.0000000000000002 steps: 1837]
[sol2_diff Success: y(100pi)=1.0 steps: 1130]

We can see that it takes far fewer iterations to solve the system, and the numerical solution is much closer to the analytical solution $y(100\pi) = 1$. If we check the residual of the original constraint, we get

plot(sol1_diff.t, sol1_diff[1, :] - sin.(sol1_diff.t), lab = "IDA")
plot!(sol2_diff.t, sol2_diff[1, :] - sin.(sol2_diff.t), lab = "DFBDF")

We see a significant numerical drift from the original constraint for both DAE solvers. Again, let's see how ModelingToolkit does.

@variables x(t) y(t)
eqs = [
    x ~ sin(t)
    D(x) ~ y
]
@named sys = ODESystem(eqs, t)
sys = complete(sys)
model = structural_simplify(sys)
prob = ODEProblem(model, [x=>0.0, y=>1.0, D(x)=>0.0, D(y)=>1.0], (0, 100pi))
sol = solve(prob, Rodas5P())
println("[sol: ", sol.retcode, ": y(100pi)=", sol[y, end], " steps: ", length(sol.t), "]")

[sol: Success: y(100pi)=1.0 steps: 2]

norm(sol[x, :] - sin.(sol.t))

0.0

We can see that ModelingToolkit handles this DAE system perfectly. $y(t)$ is completely accurate, the original constraints are satisfied without drift, and the numerical solver needs to take minimal steps. In my following lectures, we will study how ModelingToolkit handles this example system.

Numerical Integrability Criterion for DAEs

There are both differential equations and algebraic equations in acausal models. Thus, a generic acausal model is a system of differential-algebraic equations (DAEs). In general, additional processing steps are required to simulate DAEs. To see this more clearly, we will analyze the implicit Euler algorithm which is the most basic form of both Runge-Kutta and linear multistep methods for DAEs.

Given the DAE of the form

\[\begin{equation} 0 = F(u'(t), u(t), p, t), \end{equation}\]

where $F: (\mathbb{R}^n, \mathbb{R}^n, \mathbb{R}^m, \mathbb{R}) \rightarrow \mathbb{R}^n$. The implicit Euler solves for $u(t+h)$ from

\[\begin{equation} 0 = F\left(\frac{\hat{u}(t+h) - u(t)}{h},\; \hat{u}(t+h),\; p,\; t+h\right), \end{equation}\]

with fixed $t, h, p$, and $u(t)$.

Numerically, we use Newton's method to solve potentially nonlinear equations by iteratively solving the best approximating linear equations (i.e., the Jacobian of the nonlinear function with respect to the unknowns) to refine an initial guess. By the chain rule, we have

\[\begin{equation} \frac{\partial F}{\partial u} = \frac{1}{h}F_{u'} + F_{u}. \end{equation}\]

The Newton iteration is then

\[\begin{align} \frac{\partial F}{\partial u} \Delta^{[i]} &= F\left(\frac{\hat{u}^{[i]}(t+h) - u(t)}{h}, \hat{u}^{[i]}(t+h), p, t+h\right) \\ \hat{u}^{[i+1]} &= \hat{u}^{[i]} - \Delta^{[i]}, \end{align}\]

where $\{\cdot\}^{[i]}$ denotes the iteration variable at the $i$-th iteration. Thus, for the implicit Euler algorithm to work, $\lambda F_{u'} + F_{u}$ has to be non-singular for sufficiently small $\lambda\in\mathbb{R}$. In fact, this conclusion holds for all linear multistep methods and Runge-Kutta methods, that is, they all need to solve an iteration matrix of the form $\lambda F_{u'} + F_{u}$. Curious readers can check this development documentation on how the Newton iteration is set up for all the other cases. Formally, the numerical integrability criterion is

\[\begin{equation} \forall u, u', \exists \lambda > 0, \det(\lambda F_{u'} + F_{u}) \ne 0. \end{equation}\]

Structural Analysis

The above criterion is too hard to verify at compile time, since we do not yet know the exact values to solve the DAE. To make it computationally easier to check, we can change the $\forall$ to $\exists$, which is

\[\begin{equation} \exists u, u', \lambda > 0, \det(\lambda F_{u'} + F_{u}) \ne 0, \end{equation}\]

so that we only need to validate a single instance. To make the criterion even easier to check, we introduce the following definitions.

Although structurally non-singular is a strictly weaker condition than numerically non-singular, checking it seems to be more difficult on the surface because we need to come up with an example that is numerically non-singular. However, the following powerful theorem gives us a dramatically simpler way of checking if a sparse matrix is structurally non-singular.

The above condition is equivalent with checking the existence of a perfect matching on the induced bipartite graph of the incidence matrix, and it can be efficiently solved by using the augmenting path algorithm to find the maximum cardinality matching.

function construct_augmenting_path!(matching::Matching, g::BipartiteGraph, vsrc, dstfilter,
        dcolor = falses(ndsts(g)), scolor = nothing)
    scolor === nothing || (scolor[vsrc] = true)

    # if a `vdst` is unassigned and the edge `vsrc <=> vdst` exists
    for vdst in 𝑠neighbors(g, vsrc)
        if dstfilter(vdst) && matching[vdst] === unassigned
            matching[vdst] = vsrc
            return true
        end
    end

    # for every `vsrc` such that edge `vsrc <=> vdst` exists and `vdst` is uncolored
    for vdst in 𝑠neighbors(g, vsrc)
        (dstfilter(vdst) && !dcolor[vdst]) || continue
        dcolor[vdst] = true
        if construct_augmenting_path!(matching, g, matching[vdst], dstfilter, dcolor,
            scolor)
            matching[vdst] = vsrc
            return true
        end
    end
    return false
end

function maximal_matching(g::BipartiteGraph, srcfilter = vsrc -> true,
        dstfilter = vdst -> true, ::Type{U} = Unassigned) where {U}
    matching = Matching{U}(ndsts(g))
    foreach(Iterators.filter(srcfilter, 𝑠vertices(g))) do vsrc
        construct_augmenting_path!(matching, g, vsrc, dstfilter)
    end
    return matching
end

Structural Integrability Criterion for DAEs

To utilize the structural analysis framework, we need weaken the integrability criterion further from

\[\begin{equation} \exists u, u', \lambda > 0, \det(\lambda F_{u'} + F_{u}) \ne 0, \end{equation}\]

to the structural integrability criterion

\[\begin{equation} \mathfrak{I}(\mathfrak{I}(F, \{u_i'\}) + \mathfrak{I}(F, \{u_i\})) \text{ is structurally non-singular.} \end{equation}\]

Consistency Solvability Criterion for DAEs

Consider the DAE system

\[\begin{equation} F(\{x', y'\}, \{x, y\}, t) = \begin{pmatrix} f_1(x, t) \\ f_2(x', y', t) \end{pmatrix} = 0, \end{equation}\]

where $f_1$ and $f_2$ are some arbitrary smooth functions. We have

\[\begin{equation} \mathfrak{I}(F, \{x', y'\}) = \begin{pmatrix} 0 & 0 \\ 1 & 1 \end{pmatrix},\quad \mathfrak{I}(F, \{x, y\}) = \begin{pmatrix} 1 & 0 \\ 0 & 0 \end{pmatrix}. \end{equation}\]

Thus,

\[\begin{equation} \mathfrak{I}(\mathfrak{I}(F, \{u_i'\}) + \mathfrak{I}(F, \{u_i\})) = \begin{pmatrix} 1 & 0 \\ 1 & 1 \end{pmatrix} \end{equation}\]

is structurally non-singular, which means that the DAE system is structurally integrable.

However, solving for a consistent initial condition $u(t_0)$ and $u'(t_0)$ is not as simple as simply solving for $F(u'(t_0), u(t_0), t_0) = 0$. Because given a general DAE in the form of $F(u', u, p, t) = 0$, all its total time derivatives are also valid constraints, i.e.

\[\begin{align} F(u', u, p, t) &= 0 \\ F'(u'', u', u, p, t) &= 0 \\ F''(u''', u'', u', u, p, t) &= 0 \\ &\vdots \nonumber \end{align}\]

For the above example, we can differentiate the $f_1$ equation once and get

\[\begin{equation} \begin{pmatrix} f'_1(x', t) \\ f_2(x', y', t) \end{pmatrix} = 0. \end{equation}\]

Note that differentiating $f'_1$ and $f_2$ further is not necessary because we will not get additional constraints for the states ($u(t)$ and $u'(t)$) of the DAE system.

We need a more systematic way of knowing when differentiating $F$ does not add new "information" into the system. First, let's develop a characterization on the variables. Let $z=\{z_i\}$ be the set of the highest order derivative variables, and let $\lambda = \{\lambda_i\}$ contain the rest of the variables. Note that $z$ and $\lambda$ must be disjoint. Therefore, DAEs can then be written as

\[0 = F(z, \lambda, p, t) \\\]

Differentiating the above equation gives us

\[F_z z' + F_\lambda \lambda' + F_t = 0\]

Note that $z' = \{z'_i\}$ contains all the new terms generated by the differentiation, as it contains variables with higher order derivatives than before. Rearranging terms, we get

\[F_z z' = -F_\lambda \lambda' - F_t.\]

When $F_z$ is non-singular, we can use old terms to explicitly solve for $z'$, so we don't generate genuinely new equations, and this is the numerical consistency solvability criterion. Therefore, we only add new equations to the system if and only if $F_z$ is singular. It's wasteful to differentiate the entire system until the matrix is invertible; we can differentiate a minimal subset of equations to make $F_z$ non-singular.

Similarly, the structural consistency solvability criterion is then

\[\mathfrak{I}(F, z) \text{ is structurally non-singular.}\]

For the above system, the differentiated system has

\[\mathfrak{I}(F, \{x', y'\}) = \begin{pmatrix} 1 & 0 \\ 1 & 1 \end{pmatrix},\]

which is structurally non-singular. Thus, the differentiated system satisfies the structural consistency solvability criterion.

Note that the sparsity pattern of $\mathfrak{I}(\mathfrak{I}(F, \{u_i'\}) + \mathfrak{I}(F, \{u_i\}))$ is always a subset of $\mathfrak{I}(F, z)$ for arbitrary systems.

Pantelides Algorithm

We can use the Pantelides algorithm ^{[Pantelides1988]} to efficiently convert a DAE system that has structural integrability to a system with structural consistency solvability, even if it initially lacks this property.

The gist of the Pantelides algorithm is that, we can try to find an augmenting path for all equation (source) vertices on the sub-graph that only contains highest differentiated variable (destination) vertices, and if there is no augmenting path starting at an equation vertex, then we can differentiate all the equations and variables reached in the augmenting path search until there is an augmenting path starting at the differentiated equation vertex. Note that if we differentiate an equation in the form of $f(x, y, ...)$ we get

\[\frac{d}{dt}f(x, y, ...) = \frac{\partial f}{\partial x}x' + \frac{\partial f}{\partial y}y' + ...\]

Thus, the incidence of the differentiated equation is trivial to compute.

The following code comes directly from ModelingToolkit.

function pantelides!(state::TransformationState; finalize = true, maxiters = 8000)
    @unpack graph, solvable_graph, var_to_diff, eq_to_diff = state.structure
    neqs = nsrcs(graph)
    nvars = nv(var_to_diff)
    vcolor = falses(nvars)
    ecolor = falses(neqs)
    var_eq_matching = Matching(nvars)
    neqs′ = neqs
    nnonemptyeqs = count(eq -> !isempty(𝑠neighbors(graph, eq)) && eq_to_diff[eq] === nothing,
        1:neqs′)

    varwhitelist = computed_highest_diff_variables(state.structure)

    if nnonemptyeqs > count(varwhitelist)
        throw(InvalidSystemException("System is structurally singular"))
    end

    for k in 1:neqs′
        eq′ = k
        eq_to_diff[eq′] === nothing || continue
        isempty(𝑠neighbors(graph, eq′)) && continue
        pathfound = false
        # In practice, `maxiters=8000` should never be reached, otherwise, the
        # index would be on the order of thousands.
        for iii in 1:maxiters
            # run matching on (dx, y) variables
            #
            # the derivatives and algebraic variables are zeros in the variable
            # association list
            resize!(vcolor, nvars)
            fill!(vcolor, false)
            resize!(ecolor, neqs)
            fill!(ecolor, false)
            pathfound = construct_augmenting_path!(var_eq_matching, graph, eq′,
                v -> varwhitelist[v], vcolor, ecolor)
            pathfound && break # terminating condition
            if is_only_discrete(state.structure)
                error("The discrete system has high structural index. This is not supported.")
            end
            for var in eachindex(vcolor)
                vcolor[var] || continue
                if var_to_diff[var] === nothing
                    # introduce a new variable
                    nvars += 1
                    var_diff = var_derivative!(state, var)
                    push!(var_eq_matching, unassigned)
                    push!(varwhitelist, false)
                    @assert length(var_eq_matching) == var_diff
                end
                varwhitelist[var] = false
                varwhitelist[var_to_diff[var]] = true
            end

            for eq in eachindex(ecolor)
                ecolor[eq] || continue
                # introduce a new equation
                neqs += 1
                eq_derivative!(state, eq)
            end

            for var in eachindex(vcolor)
                vcolor[var] || continue
                # the newly introduced `var`s and `eq`s have the inherits
                # assignment
                var_eq_matching[var_to_diff[var]] = eq_to_diff[var_eq_matching[var]]
            end
            eq′ = eq_to_diff[eq′]
        end # for _ in 1:maxiters
        pathfound ||
            error("maxiters=$maxiters reached! File a bug report if your system has a reasonable index (<100), and you are using the default `maxiters`. Try to increase the maxiters by `pantelides(sys::ODESystem; maxiters=1_000_000)` if your system has an incredibly high index and it is truly extremely large.")
    end # for k in 1:neqs′

    finalize && for var in 1:ndsts(graph)
        varwhitelist[var] && continue
        var_eq_matching[var] = unassigned
    end
    return var_eq_matching
end