The Pantelides algorithm for delay differential-algebraic equations

We present a graph-theoretical approach that can detect which equations of a delay differential-algebraic equation (DDAE) need to be differentiated or shifted to construct a solution of the DDAE. Our approach exploits the observation that differentiation and shifting are very similar from a structural point of view, which allows us to generalize the Pantelides algorithm for differential-algebraic equations to the DDAE setting. The primary tool for the extension is the introduction of equivalence classes in the graph of the DDAE, which also allows us to derive a necessary and sufficient criterion for the termination of the new algorithm.


Introduction
We study delay differential-algebraic equations (DDAEs) of the form F (t, x(t), ẋ(t), x(t − τ )) = 0 (1.1) on the time interval [0, T ) = [0, T ] with τ > 0. The DDAE (1.1) is equipped with an (functional) initial condition of the form Differential equations with delay arise in various applications, where the current rate of change does not only depend on the current state but also on past information, for instance, population dynamics, infection disease models and laser dynamics.We refer to [8] and the references therein.Besides, if a dynamical system is controlled based on the current state of the system -a so-called feedback control law -then such a controller often requires some time to measure the current state and to compute the control action, which again results in a delay equation.On the other hand, modern modeling packages such as modelica1 or Matlab/Simulink2 automatically generate dynamical systems with constraints that arise due to interface conditions or conservation laws.As a consequence, we have to deal with equations of the form (1.1), which combine the features of delay differential equations (DDEs) [2,3,11,16] and differential-algebraic equations (DAEs) [20,24].Note that DDAEs may also arise from data-driven approaches [9,30,31], where the time delay is used to realize a transport phenomenon.Let us mention that our framework is not restricted to a single time delay, since multiple commensurate delays can be rewritten in the form (1.1) by introducing new variables [12].
One of the main difficulties for DDAEs is the fact that solutions may depend on derivatives of the function F in (1.1) as well as on evaluations of F at future time points [6,13], i.e., one has to study the interplay of the differentiation operator d/dt and the (time) shift operator ∆ τ [15], defined via (∆ τ x)(t) = x(t + τ ).In particular, it is important to understand which equations in (1.1) need to be differentiated and which need to be shifted.For linear systems, this is investigated for instance in [6, 13-15, 32, 33].For nonlinear systems a detailed analysis is not available -we instead refer to [1] for an analysis of (1.1) in a particular Hessenberg form.
In this paper, we propose a graph-theoretical approach to determine which equations need to be differentiated and which equations need to be shifted.To this end, we revisit the Pantelides algorithm [23] for DAEs in Section 2, a methodology that uses the information which variable appears in which equation to determine which equations need to be differentiated.Our main contributions are the following: (a) We introduce equivalence classes (cf.Definition 2.1) in the graph of an equation, which allows us to derive a criterion -see Theorem 2.10 -when the Pantelides algorithm terminates.We show that this criterion is equivalent to the criterion presented in [23] and thus provides another interpretation, that enables us to extend the algorithm to the DDAE case.
(b) In Section 3.1 we illustrate that, from a structural point of view, the operations differentiation and shifting are similar, such that we can essentially use the same procedure to determine the equations that need to be shifted or differentiated, respectively.
(c) We introduce a shifting and differentiation graph in Definition 3.5 and show that differentiation does not affect the shifting graph (Proposition 3.8).As a consequence, we follow the strategy in [6,32], and first determine the equations that need to be shifted.It turns out that during this procedure some of the equations need to be differentiated and we propose a linear integer program to determine which equations need to be differentiated how many times.Theorem 3.12 details that the integer program is equivalent to a standard linear program that can be solved with standard methods.
(d) Our methodology is summarized in Algorithm 4. We conclude this manuscript with a detailed analysis of Algorithm 4 in Section 4: We show in Theorem 4.1 that Algorithm 4 terminates if and only if the DDAE is structurally nonsingular with respect to a certain equivalence relation.Moreover Theorem 4.4 details in which cases Algorithm 4 concludes that no equations need to be shifted and that in this situation it produces the same result as if we apply the original Pantelides algorithm to the DDAE (1.1) where we replace ∆ −τ x with a function parameter λ.
We are convinced that our methodology can be combined with a dummy derivative approach [21] or with algebraic regularization techniques as proposed in [28,29], and thus serves as the first step to a numerical method for general nonlinear DDAEs.We emphasize that -similar to the original Pantelides algorithm [23] -our method not always determines the correct number of differentiations and shifts.For DAEs this fact is accounted for with a success check [25,26] and an algebraic analysis [28].We foresee that a similar validation of the results is also possible for DDAEs.

Notation
The natural numbers, the nonnegative integers, and the reals are denoted by N, N 0 , and R, respectively.For a set X the power set of X is denoted by P(X).For a differentiable function f : I → R n we use the notation ḟ := d dt f to denote the derivative from the right3 with respect to the (time) variable t.Similarly, the shift operator ∆ τ is defined as (∆ τ f )(t) = f (t + τ ).As a consequence, the DDAE (1.1) can be conveniently written as The partial derivative of F with respect to x, ẋ, and ∆ −τ x is denoted by ∂F ∂x , ∂F ∂ ẋ , and

Preliminary Results
A standard approach to solve differential equations with delay is the method of steps [2].
In the context of DDAEs, this requires to solve a sequence of initial value problems with a DAE, see for instance [33].More precisely, we study the DAE where One of the main differences of DAEs compared to ordinary differential equations (ODEs) is that the partial derivative of F with respect to ẋ, denoted via ∂F ∂ ẋ , may be singular.In fact, only if n = m the DAE (2.1a) is well-posed [20], and thus for the remainder of this manuscript we assume n = m.It is well-known that in general solutions of (2.1) depend on derivatives of (2.1a), see for instance [4,20,24].Following [5], we introduce the so-called derivative array of level ℓ . . .
, which can be used to determine the solution of (2.1) for large enough ℓ.The strangeness index [19,20] for instance uses D ℓ to determine matrix functions, which can be used to construct a reformulated version of (2.1a) that is suitable for numerical time integration.Although this procedure is tailored to numerical methods, it suffers from two issues: First, in a large-scale setting, the computation of the matrix functions may become computationally expensive, since all equations in (2.1a) are differentiated instead of only the required equations.Second, the computation of the matrix functions requires several numerical challenging rank decisions.
To overcome these issues, [28,29] propose to combine such algebraic regularization techniques with so-called structural analysis tools, such as the Σ-method [25,26] and the Pantelides algorithm [23].Hereby, structure is understood as the information which variable appears in which equations.In other words, the structure defines a bipartite graph G = (V E ∪V V , E) where the set of vertices V E ∪V V is given by the set of equations V E and the set of variables V V .An edge {F, x} ∈ E implies that the variable x ∈ V V appears in the equation F ∈ V E (see Definition 2.1).Note that we need to distinguish between a variable and the name of the variable in the definition of the bipartite graph.We, therefore, introduce the language which accounts for differentiation and also for shifting, which is required for the upcoming analysis of the DDAE (1.1).For notational convenience we write ξ j instead of ξ To map a vector of variable or equation names (as defined in L) to a set containing these elements, we define set : L N → P(L), (a 1 , . . ., a N ) → {a 1 , . . ., a N }.
For the remainder of this paper we identify a variable by its name and do not distinguish between both.
In the study of DDAEs it is essential to group several variables into one single vertex (see the forthcoming Section 3).We, therefore, introduce the following notation.

Definition 2.1. Consider the equation
(a) The set is called the set of variables of the equation (2.2).
(b) Let Θ denote the set of variables of (2.2) and Θ ⊆ Θ.The graph of (2.2) over the equivalence relation R ⊆ Θ × Θ is defined to be the bipartite graph where V E := set(F ) is the set of equations, V V := Θ/R denotes the set of variables and we denote the ith equation by 4) given by is visualized in Figure 2.1a (where we represent the equivalence classes by listing all its elements).By standard abuse of notation we do not distinguish between a graph and its visualization.
The main idea of the Pantelides algorithm [23]  for ℓ > 0 does not occur in any equation.For instance, in an ODE of the form ẋ = f (t, x) we can simply assign the ith equation to ẋi .
In the language of graph theory, we are thus interested in a matching (in the literature also called assignment) M ⊆ E of a graph (V, E), which is a subset of edges, such that for all vertices v ∈ V with v ∈ e 1 ∈ M and v ∈ e 2 ∈ M we have e 1 = e 2 , i.e., no vertex appears in more than one edge.An edge e = {v 1 , v 2 } ∈ M is called a matching edge and we say that v 1 and v 2 are matched.A vertex which does not occur in any matching edge is called an exposed vertex with respect to the matching.In a bipartite graph G = (V E ∪V V , E) a maximal matching is a matching such that no F ∈ V E is exposed.
Example 2.3.In the graph in Example 2.2 the variables x 1 , x 2 are not the highest derivative such that we can (temporarily) delete these variables from the graph, which we visualize by using gray edges in Figure 2.1b.A matching is for instance given by M = {{F 2 , { ẋ1 }}, {F 3 , { ẋ2 }}}, which we depict in Figure 2.1b via thick blue line.We immediately notice that F 1 is an exposed vertex (independent of any possible matching and indicated by a black vertex) and thus no maximal matching is possible.Remark 2.4.Let Θ denote the set of variables for the DAE (2.1a) (cf.Definition 2.1).Define and R equal := {(θ, θ) | θ ∈ Θ}.Then the graph of the DAE that is obtained after deleting all variables that are not highest derivative is the graph of (2.1a) over R equal .♦  Suppose that, after deleting variables that are no highest derivatives, we cannot find a maximal matching and thus have an exposed vertex F i .Then one of the following applies: The function F i does not depend on any variable, implying that in general, the DAE is not regular (see [19] for further details).Otherwise, the function F i determines a variable that is not a highest derivative, and thus F i needs to be differentiated to construct a solution.
More precisely, we replace F i by d dt F i (t, ϑ, θ) = 0 and add possible new variables.Note that all variables that appear in F i appear with one additional derivative in d dt F i (t, ϑ, θ) = 0 and thus the original variables cannot be highest derivatives and can thus be ignored in the following.
Example 2.5.Continuing with Example 2.3, we replace F 1 with the differentiated equation d dt F 1 (t, ϑ, θ) = 0, which we denote by Ḟ1 .Note that there is a (gray) edge between Ḟ1 and x 1 , since in our structural approach we do not work with the actual equation but only with the information that F 1 depends on x 1 .The chain rule then implies that Ḟ1 depends on x 1 and ẋ1 .A possible matching is given by M = {{ Ḟ1 , { ẋ1 }}, {F 3 , { ẋ2 }}}, such that F 2 is an exposed vertex with respect to M .The resulting graph is depicted in Figure 2.1c.Continuing with our strategy for Example 2.5 implies that we need to differentiate equation F 2 .This time, the situation is slightly different than before, since there exists a coupling between Ḟ1 and F 2 via the variable ẋ1 indicated by a black vertex, i.e., we have a path from Ḟ1 to F 2 via ẋ1 .
A path is a sequence of edges (e 1 , . . ., e n ) in a graph G = (V, E) such that e i = {v i−1 , v i } ∈ E, where v 0 , . . ., v n ∈ V are distinct.An alternating path with respect to a matching is a path with alternating non-matching and matching edges that starts with a non-matching edge.An augmenting path with respect to a matching is an alternating path that ends with a non-matching edge.It is easy to see that whenever an augmenting path exists, we can find another matching that includes more equations than the previous matching and hence we cannot determine (at this point) that the equation needs to be differentiated.Hence we only have to differentiate whenever we cannot find an augmenting path.In this case, we have to differentiate all equations that are connected with the exposed equation via an alternating path.For further details, we refer to [23].
Example 2.6.In Example 2.5 there exists no augmenting path that starts in F 2 (cf. Figure 2.1c).Thus, we have to differentiate F 2 and all equations F j with an alternating path from F 2 to F j .In Figure 2.1c these equations are denoted by a black dot.In this case the path (F 2 , { ẋ1 }, Ḟ1 ) implies that we have to differentiate Ḟ1 and F 2 .Note that in the resulting graph, given in Figure 2.1d, we find the maximal matching and can thus stop.
Before we present an algorithmic summary of the outline methodology, let us recall the main idea of the Pantelides algorithm, namely to match each equation in (2.1a) with a highest derivative.In more detail, we have considered the highest derivative of each variable (for the corresponding equivalence relation see Remark 2.4) to determine which equations need to be differentiated.Let {k 1 , . . ., k m } ⊆ {1, . . ., M } denote the equations that need to be differentiated.Then we define a subset of equations of (2.2) via . . .
with the understanding that ϑ contains only the subset of variables of ϑ that satisfies It Remark 2.8.The notion of structural singularity of a subset of equations is defined in [23].In contrast to [23] we restrict ourselves to the set of variables that are actually contained within this subset of equations.♦ Proposition 2.9.Consider the equation (2.2), let Θ denote the set of variables of (2.2).
For a subset Θ ⊆ Θ with equivalence relation R ⊆ Θ × Θ consider the graph G of (2.2) over Θ/R with a matching M .Let F j be an exposed node.Then one of the following is true: (a) There exists an augmenting path starting in F j .
(b) The set of equations F (t, ϑ) = 0 associated with Proof.Since all equations in C F j have an alternating path to F j we conclude that F j is the only exposed equation in C F j .This implies Suppose that there is no augmenting path that starts in F j .Assume that the inequality in (2.8) is strict.Then there exists at least one equivalence class in ( Θ ∩ set( ϑ))/R that is exposed.By definition of ( Θ ∩ set( ϑ))/R and C F j this implies that there exists an augmenting path that starts in F j , a contradiction.We conclude that we have equality in (2.8) and thus C F j is structurally singular.
Consider a proper subset CF j C F j .If F j ∈ CF j then all F k ∈ CF j are contained in a matching edge.Thus the set is structurally nonsingular.If F j ∈ CF j there exists an edge between a variable x k ∈ ( Θ ∩ set( ϑ))/R and an equation F k ∈ CF j which is not a matching edge.Thus CF j is structurally nonsingular.This completes the proof.Proposition 2.9 allows us to implement the outlined strategy as follows.Suppose that we have already found a matching for the first j − 1 equations.For the jth equation, we try to find an augmenting path.If there exists an augmenting path, then we can directly generate a matching that includes the first j equations.If not, then we differentiate all equations within the set C F j , see (2.7).Note that, to determine that there is no augmenting path, we have to check all alternating paths and hence the set C F j is automatically generated when we check for an augmenting path.One way to achieve this is via a recursive, depth-first search algorithm [7,34].The details are given in Algorithm 1.Here, we use the notation x i , F j on the one hand as vertices in V V and V E and on the other hand as indices referring to those vertices.The complete algorithm, which is known as the Pantelides algorithm [23], is presented in Algorithm 2. Note that in contrast to the original algorithm [23] we replace equations that are differentiated and do not add the equations to the graph.This has the advantage that we do not need to keep track which equation was obtained by differentiating another equation, and similarly for the variables.In Algorithm 2 we use a similar notation as in Algorithm 1 for vertices in V V and V E .The vertex x i,k represents the variable x (k) i .Since Algorithm 2 is an iterative algorithm we have to answer the questions whether Algorithm 2 terminates.We immediately observe that if Algorithm 2 terminates, then we

Algorithm 1 Augmentpath
⊲ vertices needed to be differentiated if pathfound== false Output: pathfound, assign, colorV, colorE, path pathfound ← true (pathfound, assign, colorV, colorE) ← Augmentpath(G, F k , assign, colorV, colorE) end if 15: end for have no exposed equation and thus each equation is matched with a highest derivative.We conclude that a necessary condition for termination of Algorithm 2 is that the DAE (2.1a) is structurally nonsingular with respect to Θ/R equal with Θ denoting the set of variables of (2.1a) and R equal the equivalence relation The following result shows that in fact the structural nonsingularity is also a sufficient condition.
Theorem 2.10.Consider the DAE (2.1a) and let Θ denote its set of variables and R equal the equivalence relation in (2.9).Then Algorithm 2 applied to the DAE (2.1a) terminates if and only if (2.1a) is structurally nonsingular with respect to Θ/R equal .
Proof.We show that the structural nonsingularity is equivalent to the termination criterion that is presented in the original paper [23].To this end let us rewrite (2.1a) by splitting the variable x into set(x) = set(y) ∪ set(z) such that • for all ξ ∈ set(y) it holds that ∂F ∂ ξ ≡ 0 and Using the new variables we can rewrite (2.1a) as Algorithm 2 Pantelides Algorithm for DAEs (p, q) ← (r, 0) assign, colorV, colorE) if pathfound == false then 11: for each (j, ℓ) s. t. colorE(F j,ℓ ) == true do 13: F j,ℓ ← F j,ℓ+1 14: end for 16: for each i with colorV(x i,k ) == true and F j,ℓ == assign(x i,k ) do with set(η, γ, µ) ⊆ set(y) and set( z) ⊆ set(z).Let G(γ, γ) = 0 denote the subset of (2.10b) that corresponds to the variable γ.Since the extended DAE is structurally nonsingular, we obtain and thus (2.10a) is structurally nonsingular with respect to Θ/R equal .
Even if Algorithm 2 terminates this is not a guarantee that the correct number of differentiations for each equation is determined.In fact, the Pantelides algorithm relies heavily on the so-called sparsity pattern of the DAE, i.e., each equation depends only on few variables.In particular, Algorithm 2 is not invariant under bijective transformations of the DAE (cf.[28,Remark 5]).
Example 2.11.The Pantelides algorithm applied to the DAE determines that the second equation needs to be differentiated one time.If we however add the first equation to the second equation, then we can match the first equation with ẋ1 and the second equation with ẋ2 and thus neither of the equations is differentiated.
Also, it may be the case that the Pantelides algorithm concludes that equations need to be differentiated several times, although this is not necessary [27].One possible reason for this is that the DAE is close to a high-index DAE and thus the numerical solution may be more difficult to obtain than the (potentially) small index suggests (see [27] for further details).In any case, it is possible to check that the equations are differentiated sufficiently often by applying a success check [25,26] and an algebraic criterion [28].

The Pantelides Algorithm for DDAEs
It is well-known that for DDAEs of the form (1.1) it is not sufficient to differentiate equations to obtain solutions, but in addition, some equations need to be shifted [6,13,14,32].As a direct consequence, a structural analysis for DDAEs should reveal which equations need to be shifted and which equations need to be differentiated.We thus have to understand the similarities and differences between differentiating and shifting from a structural point of view.

Differentiation vs. Shifting
To understand the effect of the two operators d dt and ∆ −τ from a structural point of view, we start with the following simple example.
Example 3.2.Let N ∈ R n×n be a nilpotent matrix with N ν = 0 and N ν−1 = 0 for some ν ∈ N and consider the two equations ) which can be understood as special cases of the Weierstraß canonical form [10] for a regular matrix pencil.It is easy to see that the solutions are of the form Apart from the operator applied to the inhomogeneity both solutions are the same.Thus the variable ẋ has in the DAE the same role as the variable y in the difference equation (3.2b).In other words, in the DAE, the variable with the highest derivative is ẋ, while in the delay equation, the variable y = ∆ 0 y is the variable with the highest shift, in the sense that 0 > −1.
Motivated by Example 3.2 we can use the Pantelides algorithm (Algorithm 2) for difference equations by replacing differentiation by shifting in Algorithm 2. In a difference equation we want to assign each equation to a different variable with highest shift, meaning that ∆ kτ x with k ≥ 0 occurs in some equation but ∆ (k+ℓ)τ x for ℓ > 0 does not.Note that by our definition the variable ∆ −τ x can never be a variable with highest shift, since this would imply that we can solve directly for this variable without shifting.This can also be seen in the solution formula (3.2c) for the difference equation, which requires an additional shift of the inhomogeneity g.
Example 3.3.The graph of the scalar difference equation is visualized in Figure 3.1a.By replacing differentiating with shifting in Algorithm 2, we shift F 1 once (cf. Figure 3.1b).Note that there is no edge between ∆ τ F 1 and ∆ −τ x 1 , since there is no chain rule for the shift operator.

The Shifting and Differentiation Graph
As outlined in Example 3.1, some equations in the DDAE (1.1) need to be shifted and some need to be differentiated.Motivated by the observation (cf.Section 3.1) that shifting and differentiating are quite similar from a structural point of view, we construct two graphs: One to determine which equations need to be shifted and one to determine which equations need to be differentiated.In the shifting graph, we do not want to match equations with delayed variables, since we cannot (directly) solve for these variables.More precisely, we do not want to match equations with variables that are not highest shifts and thus delete these variables in the shifting graph.On the other hand, we need to make sure that we do not match one equation with a variable and another equation with the derivative of that variable.This is illustrated in the following example.
Example 3.4.Consider again the DDAE (3.1).Deleting the variable ∆ τ x 2 results in the graph Figure 3.2a with maximal matching M = {{F 1 , {x 1 }}, {F 2 , { ẋ1 }}}.Note that the maximal matching is possible, since we match both equations with x 1 respectively its derivative ẋ1 .If we instead merge the variables x 1 and ẋ1 into a single vertex we obtain To ensure that we do not match different equations with the same variable (with different differentiation levels), we need to group the corresponding variables as in Example 3.4.To formalize this, consider the equation which may be understood as a generalization of (1.1) that appears after shifting and/or differentiation of some of the equations.
• The variables θ, θ ∈ Θ are called shifting similar if there exists ξ ∈ set(x) and integers k ∈ Z, p, q ∈ N 0 such that θ = ∆ kτ ξ (p) and θ = ∆ kτ ξ (q) .Shifting similarity introduces an equivalence relation on Θ, which we denote by R shift .The graph of (3.4) over R shift is called shifting graph.
• The variables θ, θ ∈ Θ are called differential similar if there exists ξ ∈ set(x) and integers p, k, ℓ ∈ N 0 such that θ = ∆ kτ ξ (p) and θ = ∆ ℓτ ξ (p) .Differential similarity introduces an equivalence relation on Θ with  A question that arises is how differentiation of equations affects the shifting graph and how shifting affects the differentiation graph (cf.[15] for an example where the order of shifting and differentiation results in different smoothness requirements).We immediately observe that the differentiation graph is affected by shifting, see the following example.
Example 3.7.The differentiation graph for the scalar DDAE does not contain any variable vertices and thus in particular no maximal matching.If we shift (3.5) we obtain 0 such that the differentiation graph is given by V = {F 1 , {x 1 }} and E = {{F 1 , {x 1 }}}.
(3.6) by replacing the ith equation (i ∈ {1, . . ., M }) in (3.4) with its total derivative with respect to t.The shifting graphs of (3.4) and (3.6) are denoted by G s := (V s E ∪V s V , E s ) and Ǧs := ( V s E ∪ V s V , Ěs ), respectively.Let us define the bijection φ : Then it is easy to see, that {u, v} ∈ E s if and only if {φ(u), φ(v)} ∈ Ěs , that is the graphs G s and Ǧs are isomorphic [17, p. 21].In particular, we have shown the following result.Since shifting affects the differentiation graph but differentiation does not affect the shifting graph, we conclude that we first determine which equations need to be shifted and afterwards determine which equations should be differentiated.
Remark 3.9.The strategy to shift first can be understood as replacing the DDAE (1.1) with a transformed DDAE that can be solved with the method of steps.For linear timeinvariant DDAEs such a strategy was investigated in [6,32] in a polynomial matrix framework.There the authors use a sequence of row compressions to determine which equations need to be shifted.It should be noted that the row-compression requires differentiation of some of the equations to enable the shift afterward.Fore more details we refer to [6,32] and the upcoming subsection.♦

Differentiating during the Shifting Step
If we find (for instance via Algorithm 2) an exposed equation in the shifting graph, we conclude that this equation and all equations that are connected with this equation by an alternating path need to be shifted.However, simply shifting all these equations may not be sufficient, since the connection may exist only implicitly via the equivalence class.We illustrate this with the following example.
Example 3.10.The shifting graph of Example 3.4, i.e., the shifting graph for the DDAE (3.1), is given in Figure 3.2b.If we match the first equation F 1 with the equivalence class {x 1 , ẋ1 } then F 2 is an exposed equation that is connected with F 1 via {x 1 , ẋ1 }.Accordingly, we have to shift both equations.If we check the graph of (3.1) (presented in Figure 3.2a), then we observe that there is no path from F 1 to F 2 , i.e., we only obtain a path if we differentiate F 2 .
To resolve an (possibly) implicit connection, we have to ensure that there also exists a path in the graph of the DDAE and not only in the shifting graph.This can be achieved by differentiating the equation that does not depend on the highest derivative in the equivalence class.Let F i k for k = 2, . . ., K denote the equations which have an alternating path to the exposed equation F i 1 and let ν k denote the number how often we have to differentiate the i k th equation.For each path5 (F i k j , v j , F i k j+1 ) in the shifting graph (j = 1, . . ., J) with v j ∈ V s V define b j ∈ Z to be the number how many more times we have to differentiate F i k j than F i k j+1 such that there exists a direct connection to a highest derivative of both equations in the graph of the DDAE.This results in the linear system with shifting graph given in Figure 3.4.Since there exists no maximal matching, we choose F 3 as exposed vertex and obtain that all equations need to be shifted.Since F 1 depends is solvable.In this case, the minimizing vector ν ⋆ of (3.10) is the unique minimizer of (3.8). Proof.
(a) Recall that J is the number of equations in the linear system Aν = b and K is the number of unknowns.We immediately observe J ≥ K − 1 (the equations and thus we conclude that rank(A) = K − 1. Suppose that the set (3.11) is not empty, i.e., there exists ν Then the feasible set of (3.10) is We conclude that ν ⋆ = ν + α ⋆ e is the unique minimizer of (3.10).
(b) Using (a) it suffices to show that (3.8) is solvable whenever (3.10) is solvable.Let ν ⋆ = [ν ⋆ k ] denote the unique minimizer of (3.10).Then the proof of (a) implies min{ν ⋆ k | k = 1, . . ., K} = 0. Since the equations (3.7) imply that the difference between two entries of ν ⋆ is an integer, we conclude ν ⋆ ∈ N K 0 .The result follows from the observation that the feasible set of (3.8) is given by Theorem 3.12 implies that we can compute the solution (provided it exists) of the linear integer program (3.8) with standard methods such as the simplex algorithm or interiorpoint methods [18,22].Alternatively, we can exploit the structure of the feasible set (3.12) and simply compute one solution of the linear system Aν = b and shift the solution accordingly.
Unfortunately, the following examples show that we cannot guarantee that the feasible set is non-empty and hence that the linear integer program (3.8) is solvable.
which can be obtained from (3.9) by adding the first equation to the second equation and renaming of f 2 .However, the additional variable in the second equation implies that we have to solve (3.8) with The first and the third equation imply that ν 1 = ν 2 = ν 3 , which poses a contradiction to the second equation and thus (3.8) is not solvable.
The reason that the linear integer program (3.8) is not always solvable is due to the fact that in order to resolve the implicit connection via the equivalence we need (at least) one equation with different levels of differentiation (see Example 3.13).Since our main goal is to determine how often we have to shift and differentiate each equation, we are interested only in the maximal number of differentiations that is required.This number may be obtained by determining the equations that are linearly dependent and then solve (3.8) by removing some of the linear dependent equations until (3.8) is solvable.Repeating this procedure with every possible combination of linear dependent equations yields a family of solutions and we can simply take the maximal value for each entry of ν.
Example 3.14.For the linear system in Example 3.13 we observe that the three equations are linearly dependent and that it is sufficient to remove one of them.If we remove the first equation, we obtain ν = (0, 0, 1), if we remove the second equation, we obtain ν = (0, 0, 0) and if we remove the third equation we obtain ν = (0, 1, 1).As a consequence, we have to differentiate the second and the third equation (in agreement with Example 3.11).
We conjecture that if the linear system Aν = b cannot be solved, then one can determine a-priori, which equations need to be removed to obtain the correct solution.This would allow reducing the computational cost significantly.

The Pantelides Algorithm for DDAEs
Before we summarize our findings in form of an algorithm, let us briefly recap the results of this section.
• We illustrated that we can use the original Pantelides algorithm [23] (see Algorithm 2) with a different equivalence class to determine which equation needs to be shifted.
Recall that we may have to differentiate some equations during the shifting step to resolve implicit connections via equivalence classes (cf.Section 3.3).To account for the additional differentiation, we present a slight modification of Algorithm 2 in Algorithm 3. We emphasize that the notation ξ i,k for ξ ∈ {F, x} has a different meaning for shifting and differentiation.In the case of shifting, the k represents the order of shift while in the case of differentiation k represents the order of differentiation.This enables us to shift in the shifting graph and differentiate in the differentiation graph with the same notation.
• Due to the structure of the equivalence classes, we shift first and then differentiate, see Proposition 3.8.Note that we do not have to resolve implicit connections within the differentiation step, since we only determine the maximal number that each equation needs to be shifted and differentiated.Since a possible implicit connection in the differentiation step results from equations that are already shifted, the implicit connection is resolved automatically by using the same equation with a smaller shift than the maximal shift that was determined in the shifting step.
The complete methodology is summarized in Algorithm 4.
is given in Figure 3.5a and in the notation of Algorithm 4 by the sets where each x i,j denotes one equivalence class.We observe that there is no maximal matching in the shifting graph and we may choose either F 1 or F 2 as exposed (cf. Figure 3.5b).Since both are connected via a matching, we shift F 1 and F 2 .Note, that F 3 has a path to F 1 and F 2 but the path is not an alternating path.Thus, we do not shift F 3 .The equations F 1 and F 2 are only implicitly connected via the equivalence class {x 1 , ẋ1 }.Constructing Algorithm 3 Generalization of Algorithm 2 for shifting and differentiation (p, q) ← (r, 0) (pathfound, assign, colorV, colorE) ← Augmentpath(( Ṽt E ∪V t V , Ẽt ), F p,q , assign, colorV, colorE) G ← updateGraph(G, type, boolFull ← true, colorV, colorE) 16: for each i with colorV(x i,k ) == true and F j,ℓ == assign(x i,k ) do Step 2: the linear integer program as described in Section 3.3 yields with the solution ν 1 = 1, ν 2 = 0. Consequently, we differentiate F 1 .In the resulting graph Figure 3.5c a possible matching is giving by M = {{F 1,1 , x 1,1 }, {F 2,1 , x 2,0 }, {F 4,0 , x 4,0 }}.The matching is not maximal and the exposed vertex F 3,1 is connected through an alternating path to F 1,1 and F 2,1 .Thus we need to shift these three equations.The linear program for these three equations with solution ν 1 = ν 2 = 1 and ν 3 = 0 yields that we differentiate the first two equations.
The resulting shifting graph Figure 3.5d provides a maximal matching.Now we construct the differentiation graph for the equations ∆ 2τ F1 , ∆ 2τ Ḟ2 , ∆ τ F 3 and F 4 given by or respectively Figure 3.6a.Note that the variable ∆ −τ x 1 does not occur in the graph since it is shifted.The vertex ∆ τ F 3 is exposed and connected through an alternating path to ∆ 2τ F1 and ∆ 2τ Ḟ2 .Thus we differentiate all three equations which results in Figure 3.6b.The existence of a maximal matching in Figure 3.6a leads to the termination of Algorithm 4. We conclude that we shift the first and second equation twice and that

Termination of the new algorithm
Similarly as for the Pantelides algorithm we have to answer the question under which circumstances Algorithm 4 applied to the DDAE (1.1) terminates.As in the DAE case (cf.Theorem 2.10), we partition the set of variables Θ of (1.1) in such a way, that variables that only differ by the level of differentiation or shifting are grouped together.In more detail, we consider the equivalence relation It is easy to see that if Algorithm 4 applied to the DDAE (1.1) terminates, then each equation is matched to a different variable and thus (1.1) is structurally nonsingular with respect to Θ/R equal .The following result shows that this is indeed also a sufficient condition.
(c) Shifting graph after differentiating F 1 and shifting Ḟ1 and F 2 with assignment showing that ∆ τ F is structurally nonsingular with respect to the equivalence relation R shift .
Note that this argument is easily extended to subsets of the general equation (3.4) and we conclude that whenever a MSS subset of equations is shifted during the shifting step it becomes structurally nonsingular.In Algorithm 3 we do not only replace F with ∆ τ F but may have to differentiate some of the equations.However, Proposition 3.8 implies that this does not affect the shifting graph and hence does not effect the structural nonsingularity.
In order to show that the shifting step terminates it suffices to show that any subset set(F ) ⊆ set(F ) that is disjoint from set ( F ) and is structurally nonsingular before we shift (4.2) remains structurally nonsingular after shifting.The only possibility for set(F ) to become structurally singular is if a subset of set(F ) has a matching edge to some vertex in set(x, ẋ)/R s .This implies that there is an alternating path from the exposed equation in F to one of the equations in set(F ), a contradiction to set( F )∩set(F ) = ∅.We conclude that the shifting step terminates.
Conversely, assume that (1.1) is structurally singular with respect to Θ/R equal .Then at some point Algorithm 4 identifies a subset (4.2) of (1.1) that is structurally singular with respect to Θ/R equal .In this case, we have set( x 2 ) = ∅ and thus In this case, shifting does not increase the set of equivalence classes and thus the set (4.2) is still structurally singular with respect to the relation R shift and thus results in an infinite loop.
To show that also the differentiation step in Algorithm 4 terminates, we want to apply Theorem 2.10.For this we have to show that the set of equations that we obtain after applying the shifting step is structurally nonsingular with respect to Θ/R equal with We conclude our analysis by investigating in which cases Algorithm 4 determines that no equation needs to be shifted.The following result shows that this is the case if the DDAE (1.1) is structurally nonsingular with respect to the set of equivalence classes that are obtained by using the relation R equal with a restricted set of variables.We conclude our analysis by emphasizing that the Pantelides algorithm for DDAEs suffers from the same problems as the original Pantelides algorithm (see for instance Example 2.11 and the discussion thereafter).In particular, there is no guarantee that the correct number of differentiations and shifts is identified.

Summary
We have presented a method (Algorithm 4) to determine which equations of the DDAE (1.1) need to be shifted and which equations need to be differentiated.The algorithm extends the Pantelides algorithm [23] for DAEs to the DDAE case.The main idea that enables the generalization is the introduction of equivalence classes in the bipartite graph associated with the DDAE (1.1).For further details see Definition 2.1.The Pantelides algorithm for DDAEs is divided into the shifting step and the differentiation step.We prove (Proposition 3.8) that differentiation does not affect the shifting graph and hence start with the shifting step.It turns out that already in the shifting step we have to differentiate some equations.To obtain the required number of differentiations for each equation, we may have to solve a sequence of linear systems, see Section 3.3 for further details.We presented a necessary and sufficient condition for the termination of Algorithm 4 in Theorem 4.1.We foresee that we can combine our framework with a dummy derivative approach [21] and an algebraic approach [28] such that Algorithm 4 can be used with numerical time integration methods.
Since our framework builds upon the Pantelides algorithm, it has similar strengths and weaknesses.One of the big advantages of our algorithm (in contrast to the algebraic index reduction procedure in [14]) is that the underlying graph-theoretical tools do not suffer from numerical rank decisions and can be computed efficiently also in a large-scale setting.
On the other hand, our algorithm is not invariant under equivalence transformations and thus may fail to determine the correct number of shifts and differentiations.In the DAE case, this fact is accounted for in the Σ-method [25,26] by including a success check.It is ongoing research to develop such a success check for DDAEs.
respectively.The union of two sets A and B is denoted by A ∪ B. If, moreover, the intersection A ∩ B of A and B is empty, then we write A ∪B.The number of elements that are contained in the set A is denoted by |A|.For an equivalence relation R ⊆ A × A we denote the set of equivalence classes by A/R.For a subset A ⊆ A we define A/R := A/ R with R := {(x, y) ∈ R | x, y ∈ A}.

Example 3 . 6 .
which we denote by R diff .The graph of (3.4) over R diff is called differentiation graph.The shifting graph for the DDAE (3.1) is given in Figure3.2b.For the differentiation graph we have Θ = {{x 1 }, { ẋ1 }} , such that the differentiation graph is given in Figure3.3.Note that {x 2 } ∈ Θ.

. 7 )( 3 . 8 ) 3 . 11 .
which we compactly write as Aν = b with ν = [ν k ] ∈ N K 0 , b = [b j ] ∈ Z J and A ∈ Z J×K .Since we do not want to differentiate equations more than necessary, we have to solve the linear integer program min K k=1 ν k , such that Aν = b, ν ∈ N 0 .Example Consider the DDAE ẋ1

Remark 3 . 15 .Example 3 . 16 .
During the shifting and differentiation step, we have to update the graph and the corresponding shifting and differentiation graphs.To unify this procedure, the necessary updates of the graph are summarized in Algorithm 5.This algorithm is called in the modified Pantelides algorithm (Algorithm 3).♦The initial shifting graph of the DDAE 0

Theorem 4 . 4 .
If the set of equations (1.1) is structurally nonsingular with respect to Θ/R equal with Θ := ξ ∈ set (x, ẋ) ∂F ∂ξ ≡ 0 then Algorithm 4 applied to (1.1) terminates and determines that no equation is shifted.In this case also Algorithm 2 applied to (1.1) where we replace ∆ −τ x with a function parameter λ terminates and the resulting graphs of both algorithms are isomorphic.Proof.First observe that the structural nonsingularity of (1.1) with respect to Θ/R equal implies structural nonsingularity with respect to Θ/R equal .Thus Theorem 4.1 ensures that Algorithm 4 terminates if applied to (1.1).Consider a subset F (t, η, γ, γ, μ, ∆ −τ x) = 0 of the DDAE (1.1) with set(η) ∪ set(γ) ∪ set(µ) ⊆ set(x) and set( x) ⊆ set(x).The structural nonsingularity with respect to Θ/R equal implies set( F ) ≤ |set (η, γ, γ, μ) /R equal | = |set (η, γ, γ, μ) /R shift |which is the structural nonsingularity of F after deleting variables with lower shifts in the shifting graph.Thus no equation is shifted.We conclude that Algorithm 4 can be reduced to the differentiation step.Within the differentiation step, Algorithm 2 and Algorithm 3 coincide and thus Algorithm 2 terminates with the same result as Algorithm 4.
Graph for (2.4) after differentiating F 1 twice and F 2 once [23]e Algorithm 4 is divided into two parts, namely the shifting step and the differentiation step, we also split the proof of Theorem 4.1 in two parts.The shifting and differentiation step are very similar to the original Pantelides algorithm.Thus we can extend the termination proof presented in[23]to our setting.The main idea in[23]is to show that differentiating a MSS subset (see Definition 2.7) renders it structurally nonsingular and a structurally nonsingular subset remains structurally nonsingular after differentiation.Let us first show that under the assumptions of Theorem 4.1 the shifting step terminates.Consider the DDAE (1.1) with set of variables Θ and equivalence relation R equal as defined in (4.1).The shifting step of Algorithm 4 applied to (1.1) terminates if and only if (1.1) is structurally nonsingular with respect to Θ/R equal .Proof.Recall that in the shifting step we group variables that are shifting similar, i.e., we work with the equivalence relation R s from Definition 3.5.By virtue of Proposition 2.9, any subset of equations of the DDAE (1.1) that Algorithm 4, or more precisely Algorithm 3, identifies to be shifted is MSS with respect to Θ/R s .Suppose that the subset are not necessarily disjoint.More precisely, decompose x = ( x 1 , x 2 ) such that set( x) = set( x 1 ) ∪ set( x 2 ), set( x 2 ) ∩ set(x, x) = ∅, and set( x 1 ) ⊆ set(x, x).Let us assume first that the DDAE (1.1) is structurally nonsingular with respect to Θ/R equal .This implies that set( x 2 ) = ∅.Shifting (4.2) results in the equation ∆ τ F given by F (t + τ, ∆ τ x, ∆ τ ẋ, x) = 0.