Conversion Methods, Block Triangularization, and Structural Analysis of Differential-Algebraic Equation Systems
BBIT manuscript No. (will be inserted by the editor)
Conversion Methods, Block Triangularization, andStructural Analysis of Differential-Algebraic EquationSystems
Guangning Tan · Nedialko S. Nedialkov · John D. Pryce
Received: date / Accepted: date
Abstract
In a previous article, the authors developed two conversion methods to im-prove the Σ -method for structural analysis (SA) of differential-algebraic equations(DAEs). These methods reformulate a DAE on which the Σ -method fails into anequivalent problem on which this SA is more likely to succeed with a genericallynonsingular Jacobian. The basic version of these methods processes the DAE as awhole. This article presents the block version that exploits block triangularizationof a DAE. Using a block triangular form of a Jacobian sparsity pattern, we identifywhich diagonal blocks of the Jacobian are identically singular and then perform aconversion on each such block. This approach improves the efficiency of finding asuitable conversion for fixing SA’s failures. All of our conversion methods can be im-plemented in a computer algebra system so that every conversion can be automated. Keywords differential-algebraic equations · structural analysis · block triangularform · modeling · symbolic computation Mathematics Subject Classification (2000) · · · This article is a continuation of [13], in which we presented two conversion methodsfor improving the Σ -method [9] for structural analysis (SA) of DAEs. When thisSA fails on a DAE with an identically singular (but structurally nonsingular) System G. TanSchool of Computational Science and Engineering, McMaster University,1280 Main Street West, Hamilton, Ontario L8S 4K1, Canada, E-mail: [email protected]. S. NedialkovDepartment of Computing and Software, McMaster University,1280 Main Street West, , Hamilton, Ontario L8S 4K1, Canada, E-mail: [email protected]. D. PryceSchool of Mathematics, Cardiff University,Senghennydd Road, Cardiff CF24 4AG, Wales, UK., E-mail: [email protected] a r X i v : . [ c s . S C ] A ug Guangning Tan et al.
Jacobian, our conversion methods reformulate the DAE into an equivalent problemon which the SA is more likely to succeed with a generically nonsingular SystemJacobian [12, 13].These two conversion methods are the linear combination (LC) method and theexpression substitution (ES) method. The former is based on replacing an existingequation by a linear combination of some equations and derivatives of them. Thelatter is based on replacing some existing derivatives by expressions that containnewly introduced variables and derivatives of them. In the ES method, the equationsthat prescribe such replacements are also appended to the original DAE, so the result-ing system is an enlarged one. The main result of a conversion using either method isa strict decrease in the value of the signature matrix [13]. Based on our experience,we conjecture that such a decrease tends to give a better problem formulation of aDAE from SA perspective.Our works [6, 10, 11] show how to construct block triangular forms (BTFs) ofa DAE using the structural data obtained from the Σ -method. A BTF indicates howeach part of the DAE influences [resp. is influenced by] other parts. The interde-pendences between all pairs of blocks may be depicted by a fine-block graph [11].Exploiting the underlying structure of a DAE, we can compute the derivatives of itssolution in a blockwise fashion [7], or perform a dummy derivative index reductionalgorithm [4, 5].We refer the reader to the previous article [13] for a summary of the Σ -method,details of its failures, the basic conversion methods, and explanations of the equiv-alence of DAEs. By “basic” we mean that these methods do not exploit BTFs of aDAE. We shall follow the notation in [13].In this article, we combine our conversion methods with a block triangularizationof a DAE and derive our block conversion methods. When the System Jacobian isidentically singular, and the DAE has a nontrivial BTF—that is, having at least twodiagonal blocks—we can identify which blocks are identically singular and performa conversion on each such block. Now that we only deal with equations and variableswithin a block, which is usually of a smaller size compared to the whole DAE, theseblock methods require fewer symbolic computations and hence are expected to bemore efficient in finding a useful conversion for fixing SA’s failures.Section 2 reviews BTFs of a sparsity pattern and BTFs of a DAE. Section 3presents our block conversion methods and demonstrates their application on a DAEfrom [1]. Section 4 gives more examples, in which the two DAEs are obtained fromelectrical circuit analysis [3]. Section 5 gives concluding remarks. In § § Σ -methodtheory, such as a signature matrix ΣΣΣ = ( σ i j ) and its value Val ( ΣΣΣ ) , a highest-value Throughout this article, “derivatives of a variable” include the variable itself as its 0th derivative.onversion Methods, Block Triangularization, and Structural Analysis of DAEs 3 transversal (HVT) T of ΣΣΣ , a valid offset pair ( c ; d ) , a System Jacobian J ( c ; d ) = ( J i j ) ,and so forth. We refer the reader to [13] for details.Terms are in slanted font at their defining occurrence. We use bold font for matri-ces that may split into blocks, and for the sub-matrices. Individual entries of a matrixare in lowercase. For example, matrix A has sub-matrices A lm and entries a i j .2.1 Block triangular forms of a sparsity pattern.Let R = n be the set of indices of n rows (equations), and let C = n be the setof indices of n columns (variables). A sparsity pattern A is a subset of the Cartesianproduct R × C that contains row-column index pairs ( i , j ) . We can view A as its in-cidence matrix ( a i j ) , where a i j equals 1 if ( i , j ) ∈ A and 0 otherwise. A transversal of A is n positions in A with exactly one position in each row and each column. If A has some transversal, then it is structurally nonsingular . The union of all transver-sals of A comprise its essential sparsity pattern A ess [11]. Obviously, A is structurallynonsingular if and only if A ess is nonempty.Assume henceforth that A is structurally nonsingular. Let P and Q be two suitablepermutation matrices for A , such that the permuted incidence matrix A (cid:48) = P A Q canbe written in a p × p block form A (cid:48) = A A · · · A p A · · · A p . . . ... A pp , (2.1)where each diagonal block A qq , q = p , is square of positive size N q . We say theblock form (2.1) is a BTF of A . Blanks in (2.1) mean that a sub-matrix A kl below theblock diagonal with k > l is empty.A sparsity pattern is irreducible , if it cannot be permuted to the form (2.1) with p > reducible . A BTF is irreducible if each diagonal block is irreducible ; otherwise it is reducible [11]. Hence, if (2.1) is irreducible, then p is thelargest number of diagonal blocks among all possible BTFs of A (cid:48) .When we say block q of a matrix in a BTF, we shall refer to the q th diagonalblock submatrix. For q = p , we define for block q the index set B q = the set of indices i that belong to block q . Throughout this article, they are the indices of the permuted A (cid:48) , not those of theoriginal A .Another useful notation is blockOf ( i ) that denotes the block number q such thatindex i ∈ B q . Since each diagonal block is square, both B q and blockOf ( i ) notationapply to rows and columns equally. To summarize, for i ∈ n and q ∈ p ,blockOf ( i ) = q ⇐⇒ i ∈ B q ⇐⇒ q − ∑ m = N m + ≤ i ≤ q ∑ m = N m . The colon notation p : q for integers p , q denotes either the unordered set or the enumerated list ofintegers i with p ≤ i ≤ q , depending on context. Guangning Tan et al. × ×× ×× ×× ×× ×× × × ×× ×× ×× ×× ×× × q N q B q { } { } { } { , , } (a) (b) Fig. 2.1: (a) Two nontrivial BTFs of the same sparsity pattern. The left one is re-ducible with number of blocks p =
2. The right one is irreducible with p =
4. (b)Block information for the irreducible BTF.
Example 2.1
We illustrate in Figure 2.1 the above block notation with a sparsity pat-tern of two nontrivial BTFs.The following lemma connects the transversals of a sparsity pattern A and thetransversals of its diagonal blocks in some BTF. Lemma 2.1 [11, Lemma 2.4] Any transversal T of a sparsity pattern A is containedin the union of the diagonal blocks of any BTF of A , that is, T ⊆ A ∪ · · · ∪ A pp . Equivalently, the intersection of T with block q of A is a transversal T q of A qq .2.2 Block triangular forms of a DAE.The natural sparsity pattern of a DAE indicates if a variable x j occurs in an equation f i . Each such occurrence corresponds to a finite entry σ i j in ΣΣΣ , and hence we have S = (cid:8) ( i , j ) | σ i j > − ∞ (cid:9) (the sparsity pattern of ΣΣΣ ) . If S has some transversal, then ΣΣΣ has a transversal with finite σ i j ’s and a finite Val ( ΣΣΣ ) [9], so the DAE is structurally well posed (SWP) [10]; otherwise it is structurally illposed. Here, we shall deal with the SWP case only.A more informative BTF derives from the sparsity pattern S = S ( c ; d ) of aSystem Jacobian J = J ( c ; d ) as defined in [13, (2.6)]: S = (cid:8) ( i , j ) | d j − c i = σ i j (cid:9) (the sparsity pattern of J ) . (2.2)By [13, (2.2)], d j − c i = σ i j holds on a HVT T of ΣΣΣ , so T is also a transversal of S .A less obvious set contains the positions that contribute to det ( J ) : S ess = the union of all HVTs of ΣΣΣ (the essential sparsity pattern of
ΣΣΣ ) , which is also the essential sparsity pattern of S for any valid offset pair ( c ; d ) [11,Lemma 3.1].Since d j − c i = σ i j holds on each HVT and hence implies σ i j > − ∞ , we have S ess ⊆ S ⊆ S for any offset pair ( c ; d ) . onversion Methods, Block Triangularization, and Structural Analysis of DAEs 5 Our experience suggests that the (irreducible) BTF based on S can be significantlyfiner than that based on S . We refer to the former BTF as fine BTF , and to the latteras coarse BTF . We refer to the diagonal blocks in the fine BTF as fine blocks , andrefer to those in the coarse BTF as coarse blocks .Assume that S is permuted into a p × p BTF. Following this BTF, we apply thesame permutations on J and ΣΣΣ , and write them in p × p block forms: J = J J · · · J p · · · J p ... . . . . . . ... · · · pp and ΣΣΣ = ΣΣΣ ΣΣΣ · · · ΣΣΣ p ΣΣΣ ΣΣΣ · · · ΣΣΣ p ... . . . . . . ... ΣΣΣ p · · · · · · ΣΣΣ pp . (2.3)We call this procedure a block triangularization of the DAE, and note that the sparsitypattern S of ΣΣΣ may not be in the same BTF as S of J . That is, every σ i j below theblock diagonal of ΣΣΣ is not necessarily − ∞ , but must satisfy σ i j < d j − c i as J i j ≡ d j − c i (cid:40) > σ i j if blockOf ( j ) < blockOf ( i ) ≥ σ i j if blockOf ( j ) ≥ blockOf ( i ) . (2.4)We refer the reader to [6, 11] for more details on BTFs. Example 2.2
We illustrate the coarse and fine BTFs with the (artificially) modifieddouble pendula DAE in [7]. The state variables are x , y , λ , u , v , µ ; G is gravity, L > α is a constant.0 = f = x (cid:48)(cid:48) + x λ = f = u (cid:48)(cid:48) + u µ = f = y (cid:48)(cid:48) + y λ + ( x (cid:48) ) − g = f = ( v (cid:48)(cid:48)(cid:48) ) + v µ − G = f = x + y − L = f = u + v − ( L + αλ ) + λ (cid:48)(cid:48) ΣΣΣ = v µ u x y λ c i f • f • f • f • f • f • d j J = v (cid:48)(cid:48)(cid:48) µ u (cid:48)(cid:48) x ( ) y ( ) λ ( ) f v (cid:48)(cid:48)(cid:48) vf u f (cid:48)(cid:48) u f ( ) x yf ( ) yf ( ) x The row and column labels in J , showing equations and variables differentiated toorder c i and d j , aim to remind the reader of the formula for J in [13, (2.6)].There are two 3 × f , f , f and variables v , µ , u , can further decompose into three 1 × f , f , f and variables x , y , λ , is irreducible.Hence there are four blocks in the fine BTF. Guangning Tan et al.
The sparsity pattern S of J is exactly the one in Figure 2.1(a), so the fine BTFinformation is in Figure 2.1(b).If we state Lemma 2.1 in the context of a Jacobian sparsity pattern, then we havethe following lemma. Lemma 2.2 [11, Lemma 3.3] Assume that a Jacobian sparsity pattern S is in someBTF. Let ( ΣΣΣ qm ) q , m = p be the corresponding sub-matrices of ΣΣΣ . Then a HVT T of
ΣΣΣ is the union of HVTs T q of the diagonal blocks ΣΣΣ qq : T = T ∪ · · · ∪ T p . This lemma is not difficult to prove, given that a transversal T of S is the unionof transversals T q of the diagonal blocks of S .The following lemma is useful for proving the main Theorems 3.1 and 3.2 of theblock conversion methods in § Lemma 2.3
Assume that
ΣΣΣ has a finite
Val ( ΣΣΣ ) and is in a p × p block form as in(2.3). Let c and d be two nonnegative integer n-vectors. Assume also that(a) d j − c i > σ i j holds for all entries below the diagonal blocks of ΣΣΣ ,(b) d j − c i ≥ σ i j holds elsewhere, and(c) Val ( ΣΣΣ ) = ∑ j d j − ∑ i c i .Then(i) ( c ; d ) is a valid offset pair of ΣΣΣ ,(ii) the block form of
ΣΣΣ is a BTF of the Jacobian sparsity pattern S , and(iii) a HVT of ΣΣΣ is the union of HVTs T q of the diagonal blocks ΣΣΣ qq , for all q = p.Proof (i) We let T denote a HVT of ΣΣΣ . Since Val ( ΣΣΣ ) is finite, σ i j ≥ ( i , j ) ∈ T .For ( c ; d ) to be a valid offset of ΣΣΣ , d j − c i ≥ σ i j must hold for all i , j = n , withequalities for all ( i , j ) ∈ T [9].By (a) and (b), d j − c i ≥ σ i j holds everywhere. Summing these inequalities over T gives ∑ ( i , j ) ∈ T ( d j − c i ) ≥ ∑ ( i , j ) ∈ T σ i j . The left-hand side equals ∑ j d j − ∑ i c i , and the right-hand side equals Val ( ΣΣΣ ) bydefinition. By (c), these two values are equal, so d j − c i = σ i j holds for all ( i , j ) ∈ T ,and ( c ; d ) is valid for ΣΣΣ . (ii) By (a), the blocks below the block diagonal in S , derived from ΣΣΣ and ( c ; d ) using (2.2), are empty. By the definition of a BTF of a Jacobian sparsity pattern, S is in a BTF as described by the p × p block form. (iii) This follows immediately from (ii) and Lemma 2.2. (cid:117)(cid:116)
Following a p × p BTF based on S , we can write any valid offset pair ( c ; d ) of ΣΣΣ in a block form as ( c ; d ) , ( c ; d ) , . . . , ( c p ; d p ) , (2.5)where each of the sub-vectors c q and d q is of length N q , where q = p . onversion Methods, Block Triangularization, and Structural Analysis of DAEs 7 Lemma 2.4
Assume that a Jacobian pattern S , derived by ΣΣΣ and a valid offset pair ( c ; d ) , is in some BTF. If we write ( c ; d ) into block form as in (2.5), then ( c q ; d q ) is avalid offset pair of ΣΣΣ qq .Proof Let T be a HVT of ΣΣΣ . By Lemma 2.2, the intersection of T with block q is aHVT T q of ΣΣΣ qq . Then d j − c i = σ i j holds for all ( i , j ) ∈ T q ⊆ T . Since ( c ; d ) is valid for ΣΣΣ , d j − c i ≥ σ i j and c i ≥ ΣΣΣ qq , where i , j ∈ B q . Thus the offset pair ( c q ; d q ) matched to block q satisfies the conditions [13, (2.2)] for being valid for ΣΣΣ qq . (cid:117)(cid:116) From the view of Lemma 2.4, we can regard each diagonal block
ΣΣΣ qq as a sig-nature matrix in its own right. Equivalently, each block q , having N q equations in N q variables, can be viewed as a sub-DAE, with a signature matrix ΣΣΣ qq , a finite Val ( ΣΣΣ qq ) ,a local offset pair ( c q ; d q ) , and a sub-Jacobian J qq . Expressions that contribute to en-tries in an off-diagonal block ΣΣΣ qm , where q (cid:54) = m , can be considered as driving terms,or equivalently, the influence of variables in block m on those in block q . We refer to ( c ; d ) of ΣΣΣ as a global offset pair . The reader is referred to [11] for more theoreticalresults about block triangularization and global/local offset pairs.
They are suitable for improving the efficiency of finding a useful conversion for fixingSA’s failures. If J is identically singular, then by (2.3), det ( J ) = ∏ pq = det ( J qq ) ≡
0, soat least one J qq for some q ∈ p is identically singular. As discussed before, we canregard block q as a sub-DAE with a signature matrix ΣΣΣ qq . Then we wish to apply thebasic conversion methods on this sub-DAE to achieve a strict decrease in Val ( ΣΣΣ qq ) ,provided the conditions for applying these methods are satisfied for those variablesand equations within block q .However, what we should ensure is a strict decrease in the value of the whole signature matrix, namely Val ( ΣΣΣ ) < Val ( ΣΣΣ ) , where ΣΣΣ is the signature matrix of theresulting DAE. Proving this inequality from a decrease in Val ( ΣΣΣ qq ) is nontrivial,because a conversion on block q may affect blocks ΣΣΣ qm for m = , . . . , q − , q + , . . . , p . Especially in the ES method, ΣΣΣ qq and these blocks are enlarged. Hence,the conditions and the conversion process need to be carefully modified, so that theconversion methods can adapt to a BTF based on S .We give an introductory example in § § § S with some valid ( c ; d ) , but also to anyBTF of S . For example, the basic conversion methods consider a DAE in a (trivial)BTF of one n × n block. Guangning Tan et al. = f = x + x + h ( t ) = f = x + ( x (cid:48) + x (cid:48) ) x (cid:48) + h ( t ) = f = x (cid:48) + h ( t ) . (3.1) ΣΣΣ = x x x c i (cid:34) (cid:35) f • f • f • d j J = x (cid:48) x (cid:48) x (cid:48) (cid:34) (cid:35) f (cid:48) f x (cid:48) x (cid:48) x (cid:48) + x (cid:48) f h , h , h are driving functions.) The coarse BTF and the fine BTF are identical,both having two diagonal blocks.In the basic LC method, we can choose u = [ − x (cid:48) , , − x (cid:48) − x (cid:48) ] T ∈ coker ( J ) . Using[13, (4.3)], we have I = (cid:8) i | u i (cid:54)≡ (cid:9) = (cid:8) , , (cid:9) , c = min i ∈ I c i = , L = (cid:8) l ∈ I | c l = c (cid:9) = (cid:8) , (cid:9) . We let σ ( x j , u ) denote the order of the highest derivative to which x j occurs in u , or − ∞ if x j does not occur in u [13, (4.1)]. The LC condition [13, (4.4)] is violated since σ ( x j , u ) = (cid:54) < = d j − c for all j = . Not surprisingly, replacing either f or f by f = ∑ i ∈ I u i f ( c i − c ) i = − x (cid:48) h (cid:48) ( t ) + (cid:0) x + h ( t ) (cid:1) − ( x (cid:48) + x (cid:48) ) (cid:0) x (cid:48) + h ( t ) (cid:1) does not result in a decrease in Val ( ΣΣΣ ) ; verifying this is not difficult.Notice that only the sub-Jacobian of block 1, J = ∂ ( f (cid:48) , f ) / ∂ ( x (cid:48) , x (cid:48) ) , is sin-gular. Suppose we consider block 1, with B = (cid:8) , (cid:9) , as a sub-DAE, and choose u = [ − x (cid:48) , ] T ∈ coker ( J ) . Within block 1, the LC method derives I = (cid:8) i ∈ B | u i (cid:54)≡ (cid:9) = (cid:8) , (cid:9) , c = min i ∈ I c i = , L = (cid:8) l ∈ I | c l = c (cid:9) = (cid:8) (cid:9) . Now the LC condition [13, (4.4)] is satisfied for the column indices in block 1: σ ( x j , u ) = − ∞ < d j − c for j = , ∈ B . Replacing f by f = u f (cid:48) + u f = x + h ( t ) − x (cid:48) h (cid:48) ( t ) results in the DAE with thefollowing SA result. onversion Methods, Block Triangularization, and Structural Analysis of DAEs 9 ΣΣΣ = x x x c i (cid:34) (cid:35) f • f • f • d j J = x x x (cid:48) (cid:34) (cid:35) f f g (cid:48) ( t ) f J is nonsingular. The conversion results in a decrease in the valueof the signature matrix: Val ( ΣΣΣ ) = < = Val ( ΣΣΣ ) .The basic ES method can work on (3.1) by choosing v = [ , − , ] T ∈ ker ( J ) .It is simpler—though trivial for this example—to work on block 1 only. We find v = [ , − ] T ∈ ker ( J ) , and use [13, (4.10)] to derive J = (cid:8) l ∈ B | v l (cid:54)≡ (cid:9) = (cid:8) , (cid:9) , s = | J | = , M = (cid:8) , (cid:9) , c = max i ∈ M c i = . Since v is constant, it is not difficult to verify that the ES conditions [13, (4.11)] hold.We choose l = ∈ J and introduce for x a new variable y = x ( d − c ) − v v · x ( d − c ) = x + x . The ES method hence says: replace x by y − x in f , and replace x (cid:48) by y (cid:48) − x (cid:48) in f . Finally we append the equation g that prescribes such replacements, and obtain0 = f = y + h ( t ) = f = x (cid:48) + h ( t ) = f = x + y (cid:48) x (cid:48) + h ( t ) = g = − y + x + x . ΣΣΣ = x x y x c i g • f • f • f • d j J = x x y (cid:48) x (cid:48) g − f x (cid:48) y (cid:48) f (cid:48) f ( ΣΣΣ ) = < = Val ( ΣΣΣ ) , and the SA succeeds as det ( J ) = r denotethe zero column vector of size r . Assume that a J qq is identically singular. Let (cid:98) u ∈ coker ( J qq ) , where (cid:98) u (cid:54)≡ N q . Let also u = N + ··· + N q − (cid:98) u0 N q + + ··· + N p . Denote I = (cid:8) i | u i (cid:54)≡ (cid:9) ⊆ B q , c = min i ∈ I c i , L = (cid:8) l ∈ I | c l = c (cid:9) , and L = (cid:8) l ∈ L | u l is (nonzero) constant (cid:9) . (3.2)The set L is used to seek a conversion that guarantees equivalence between the orig-inal DAE and the converted one. The block LC method is based on the followingtheorem. Theorem 3.1 If σ ( x j , u ) < d j − c for all j ∈ B q (3.3) and we replace an equation f l , l ∈ L, byf = ∑ i ∈ I u i f ( c i − c ) i , then Val ( ΣΣΣ ) < Val ( ΣΣΣ ) , where ΣΣΣ = ( σ i j ) is the signature matrix of the resulting DAE. Before proving this theorem, we show how to apply the block LC method andprove a related lemma.
Example 3.1
We illustrate the block LC method with the Campbell-Griepentrog two-link robot arm DAE [1]. We slightly simplify the problem formulation to (3.4), allow-ing the first-order derivatives x (cid:48) , x (cid:48) , and x (cid:48) to occur implicitly in the equations. Thetwo state variables u and u in the original formulation are renamed x and x , re-spectively (and not to be confused with entries in a vector u in our notation).The equations of this problem are0 = A = x (cid:48)(cid:48) − (cid:104) c ( x )( x (cid:48) + x (cid:48) ) + x (cid:48) d ( x ) + ( x − x ) (cid:0) a ( x ) + b ( x ) (cid:1) + a ( x )( x − x ) (cid:105) = B = x (cid:48)(cid:48) − (cid:104) − c ( x )( x (cid:48) + x (cid:48) ) − x (cid:48) d ( x ) + ( x − x ) (cid:0) − a ( x ) − b ( x ) (cid:1) − a ( x ) x + (cid:0) a ( x ) + (cid:1) x (cid:105) = C = x (cid:48)(cid:48) − (cid:104) − c ( x )( x (cid:48) + x (cid:48) ) − x (cid:48) d ( x ) + ( x − x ) (cid:0) a ( x ) − b ( x ) (cid:1) − x (cid:48) c ( x ) − d ( x ) (cid:0) x (cid:48) + x (cid:48) (cid:1) − (cid:0) a ( x ) + b ( x ) (cid:1) ( x − x ) (cid:105) = D = cos x + cos ( x + x ) − p ( t ) = E = sin x + sin ( x + x ) − p ( t ) , (3.4)where a ( θ ) = / ( − cos θ ) c ( θ ) = sin θ / ( − cos θ ) p ( t ) = cos ( − e t ) + cos ( − t ) b ( θ ) = cos θ / ( − cos θ ) d ( θ ) = sin θ cos θ / ( − cos θ ) p ( t ) = sin ( − e t ) + sin ( − t ) . onversion Methods, Block Triangularization, and Structural Analysis of DAEs 11 ΣΣΣ = x x x x x c i f B • f C • f A • f D • f E • d j J = x (cid:48)(cid:48) x x x (cid:48)(cid:48) x (cid:48)(cid:48) B a − a − C a + b − a − b A − a a D (cid:48)(cid:48) ∂ D ∂ x ∂ D ∂ x E (cid:48)(cid:48) ∂ E ∂ x ∂ E ∂ x Here in J , a = a ( x ) = / ( − cos x ) b = b ( x ) = cos x / ( − cos x ) ∂ D / ∂ x = − sin x − sin ( x + x ) ∂ D / ∂ x = − sin ( x + x ) ∂ E / ∂ x = cos x + cos ( x + x ) ∂ E / ∂ x = cos ( x + x ) . The DAE (3.4) is of differentiation index 5, while the SA reports structural index ν S =
3. Hence this must be a failure case, because ν S is an upper bound for thedifferentiation index when the SA succeeds [9]. We can see that the sub-Jacobian J of block 2 is identically singular.Our method first computes (cid:98) u = [ , + cos x ] T ∈ coker ( J ) . Then u = [ , , + cos x , , ] T . Using (3.2), we have I = (cid:8) i | u i (cid:54)≡ (cid:9) = (cid:8) , (cid:9) , c = min i ∈ I c i = , L = { , } , and L = { } . The variables x and x in block 2 do not occur in u , so the condition (3.3) is satisfied.Considering equivalence, we pick l = ∈ L over l = ∈ L \ L , and replace f l = C by C = u C + u A = C + ( + cos x ) A . The SA results of the resulting DAE are asfollows. ΣΣΣ = x x x x x c i A • B • C • D • E • d j J = x x x (cid:48)(cid:48) x ( ) x ( ) A − a a B a − a − C (cid:48)(cid:48) ∂ C ∂ x + cos x D ( ) ∂ D ∂ x ∂ D ∂ x E ( ) ∂ E ∂ x ∂ E ∂ x Here ∂ C / ∂ x = ( a − a b + b )( − cos x ) . The SA reports ν S = ( J ) = ( a − a b + b ) sin x (cid:54) = . Now Val ( ΣΣΣ ) = < = Val ( ΣΣΣ ) . Lemma 3.1
Consider a BTF of a Jacobian pattern S derived from ΣΣΣ and ( c ; d ) . Ifwe perform the LC conversion as described in Theorem 3.1, then in the resulting ΣΣΣ ,d j − c i (cid:40) > σ i j if blockOf ( j ) < blockOf ( i ) ≥ σ i j if blockOf ( j ) ≥ blockOf ( i ) . (3.5) Proof
We only replace f l by f l = f in a conversion, so σ i j = σ i j for all i (cid:54) = l and all j . By (2.4), (3.5) holds for all i (cid:54) = l .When i = l , we consider two cases: (a) blockOf ( j ) < q and (b) blockOf ( j ) ≥ q .(a) blockOf ( j ) < q = blockOf ( l ) . By (2.4), σ l j < d j − c l . Then σ l j is σ (cid:0) x j , f l (cid:1) = σ (cid:16) x j , ∑ i ∈ I u i f ( c i − c ) i (cid:17) ≤ max (cid:110) σ ( x j , u ) , max i ∈ I σ (cid:0) x j , f ( c i − c ) i (cid:1)(cid:111) . (3.6)We use some simple derivations to obtain σ ( x j , u ) ≤ σ ( x j , J qq ) ≤ max i ∈ I σ ( x j , f i ) = max i ∈ I σ i j < max i ∈ I ( d j − c i ) = d j − min i ∈ I c i = d j − c l and (3.7a)max i ∈ I σ (cid:16) x j , f ( c i − c ) i (cid:17) = max i ∈ I ( σ i j + c i − c ) < d j − c = d j − c l . (3.7b)Substituting (3.7a) and (3.7b) in (3.6), we obtain σ l j = σ (cid:0) x j , f l (cid:1) < d j − c l .(b) blockOf ( j ) ≥ q = blockOf ( l ) . By (2.4), σ l j ≤ d j − c l . We can replace the two“ < ” in (3.7) by “ ≤ ”, and using these inequalities in (3.6), we have σ l j ≤ d j − c l . (cid:117)(cid:116) Using Lemma 3.1, we can now prove Theorem 3.1.
Proof
By Lemma 2.4, we can regard block q as a sub-DAE with ΣΣΣ qq and ( c q ; d q ) .The conversion described in Theorem 3.1 can be considered as an application ofthe basic LC method to this sub-DAE. Since the block LC condition (3.3) holds,that is, σ ( x j , u ) < d j − c for all j ∈ B q that belong to this sub-DAE, the basic LCcondition [13, (4.4)] also holds for the sub-DAE. Hence Val ( ΣΣΣ qq ) < Val ( ΣΣΣ qq ) .Let T be a HVT of ΣΣΣ . Using (3.5) in Lemma 3.1, we have Val ( ΣΣΣ ) = ∑ ( i , j ) ∈ T σ i j ≤ d j − c i = Val ( ΣΣΣ ) . Now we prove Val ( ΣΣΣ ) < Val ( ΣΣΣ ) by contradiction. First assume thatVal ( ΣΣΣ ) = ∑ j d j − ∑ i c i = Val ( ΣΣΣ ) ≥
0. With (3.5) , the three conditions in Lemma 2.3are satisfied. From this lemma, it follows that the Jacobian patterns S , derived from ΣΣΣ and ( c ; d ) , and S , derived from ΣΣΣ and ( c ; d ) , are in the same p × p BTF.By Lemma 2.2, T is the union of HVTs T m of all diagonal blocks ΣΣΣ mm , m = p .By the construction of ΣΣΣ , Val ( ΣΣΣ mm ) = Val ( ΣΣΣ mm ) for all m (cid:54) = q . Then a contradictionfollows fromVal ( ΣΣΣ ) = ∑ ( i , j ) ∈ T σ i j = p ∑ m = ∑ ( i , j ) ∈ T m σ i j = p ∑ m = Val ( ΣΣΣ mm )= ∑ m (cid:54) = q Val ( ΣΣΣ mm ) + Val ( ΣΣΣ qq ) < ∑ m (cid:54) = q Val ( ΣΣΣ mm ) + Val ( ΣΣΣ qq )= p ∑ m = Val ( ΣΣΣ mm ) = p ∑ m = ∑ ( i , j ) ∈ T m σ i j = ∑ ( i , j ) ∈ T σ i j = Val ( ΣΣΣ ) , (3.8)where T is a HVT of ΣΣΣ and T m are HVTs of its diagonal blocks ΣΣΣ mm . The assumptionVal ( ΣΣΣ ) =
Val ( ΣΣΣ ) is hence false, so Val ( ΣΣΣ ) < Val ( ΣΣΣ ) holds. onversion Methods, Block Triangularization, and Structural Analysis of DAEs 13 J qq is identically singular. Let (cid:98) v ∈ ker ( J qq ) , where (cid:98) v (cid:54)≡ N q .Similarly, we construct the column n -vector v = N + ··· + N q − (cid:98) v0 N q + + ··· + N p . We use notation similar to that used in the basic ES method (see [13, § J = (cid:8) j | v j (cid:54)≡ (cid:9) ⊆ B q , M = (cid:8) i ∈ B q | d j − c i = σ i j for some j ∈ J (cid:9) , s = | J | , c = max i ∈ M c i and J = (cid:8) l | v l is (nonzero) constant (cid:9) . (3.9)The set J is used to seek a conversion that guarantees equivalence between the origi-nal DAE and the converted one. The conditions for applying the block ES method are σ ( x j , v ) (cid:40) < d j − c if j ∈ J or blockOf ( j ) < q ≤ d j − c if j ∈ B q \ J or blockOf ( j ) > q , and d j − c ≥ j ∈ J . (3.10)We choose an l ∈ J , and introduce s − y j = x ( d j − c ) j − v j v l · x ( d l − c ) l for all j ∈ J \ (cid:8) l (cid:9) . (3.11)In each f i with i ∈ B q , wereplace each x ( σ ij ) j with d j − c i = σ i j and j ∈ J \ (cid:8) l (cid:9) by (cid:16) y j + v j v l · x ( d l − c ) l (cid:17) ( c − c i ) . (3.12)Note that because of M in (3.9), we actually perform replacements (equivalently re-ferred to as “expression substitutions”) in only f i ’s with i ∈ M ⊆ B q . Denote each new f i by f i , and let also f i = f i for the unchanged equations with i / ∈ M .By (3.11), we append s − = g j = − y j + x ( d j − c ) j − v j v l · x ( d l − c ) l for all j ∈ J \ (cid:8) l (cid:9) . (3.13)The block ES method is based on the following theorem. Theorem 3.2
Let J, s, M, and c be as defined in (3.9). Assume that the conditions(3.10) hold. For an l ∈ J, if we1) introduce s − new variables x j , j ∈ J \ (cid:8) l (cid:9) , as defined in (3.11),2) perform replacements in f i , for all i ∈ B q , as described in (3.12), and3) append s − equations g j , j ∈ J \ (cid:8) l (cid:9) , as defined in (3.13),then Val ( ΣΣΣ ) < Val ( ΣΣΣ ) , where ΣΣΣ is the signature matrix of the resulting DAE.
Before proving this theorem, we illustrate the block ES method with the robotarm DAE (3.4) and prove two related lemmas.
Example 3.2
The method finds first (cid:98) v = [ , ] T ∈ ker ( J ) . Then v = [ , , , , ] T .Using (3.9), we have J = J = (cid:8) j | v j (cid:54)≡ (cid:9) = { , } , s = | J | = , M = { , } , c = max i ∈ M c i = . Since v is constant, J = J and the first condition in (3.10) holds. The second conditionin (3.10) holds also, as d − c = d − c =
0. We choose x , whose column index in thepermuted ΣΣΣ is l = ∈ J . Then we introduce for x , the other variable in block 2 withcolumn index j =
3, a new variable y = x ( d − c ) − v v · x ( d − c ) = x − x . Correspondingly, we append 0 = g = − y + x − x and replace x by y + x in C and A , the equations in block 2.The resulting DAE has the following new equations0 = A = x (cid:48)(cid:48) − (cid:104) c ( x )( x (cid:48) + x (cid:48) ) + x (cid:48) d ( x ) + ( x − x ) (cid:0) a ( x ) + b ( x ) (cid:1) + a ( x ) y (cid:105) = C = x (cid:48)(cid:48) − (cid:104) − c ( x )( x (cid:48) + x (cid:48) ) − x (cid:48) d ( x ) + ( x − x ) (cid:0) a ( x ) − b ( x ) (cid:1) − x (cid:48) c ( x ) − d ( x ) (cid:0) x (cid:48) + x (cid:48) (cid:1) − (cid:0) a ( x ) + b ( x ) (cid:1) y (cid:105) = g = − y + x − x . ΣΣΣ = x x x y x x c i g • B • C • A • D • E • d j J = x x x (cid:48)(cid:48) y (cid:48)(cid:48) x ( ) x ( ) g − B a − a − C (cid:48)(cid:48) a − b − a − b A (cid:48)(cid:48) a + b a D ( ) ∂ D ∂ x ∂ D ∂ x E ( ) ∂ E ∂ x ∂ E ∂ x Now the System Jacobian J is generically nonsingular. The SA reports the correctindex 5, and succeeds at any point where det ( J ) = ( a − a b + b ) sin x (cid:54) = ( ΣΣΣ ) = < = Val ( ΣΣΣ ) .In [8], Pryce fixed the SA’s failure on (3.4), and pointed out that only the intro-duction of x − x as a separate variable is essential to his fix. Example 3.2 verifiesPryce’s argument and shows that the block ES method finds his reformulation in asystematic way.To prove Theorem 3.2, we shall use the following two assumptions.(a) Without loss of generality, we assume the entries (cid:98) v j (cid:54)≡ s positionsof (cid:98) v , that is, (cid:98) v = [ (cid:98) v , . . . , (cid:98) v s , , . . . , ] T . By (3.9), J = ∑ q − m = N m + ∑ q − m = N m + s . onversion Methods, Block Triangularization, and Structural Analysis of DAEs 15 (b) We introduce one more variable y l = x ( d l − c ) l for the chosen l ∈ J , and appendcorrespondingly one more equation 0 = g l = − y l + x ( d l − c ) l .We show first that the signature matrix ΣΣΣ of the resulting DAE can be put in theblock structure as shown in Figure 3.1. Then we construct two ( n + s ) -vectors (cid:101) c and (cid:101) d in (3.14), and prove in Lemma 3.2 that (cid:101) d j − (cid:101) c i > σ i j holds below the block diagonal,while (cid:101) d j − (cid:101) c i ≥ σ i j holds elsewhere. Lastly, we prove Theorem 3.2. Σ , ··· Σ , q − Σ , q −∞ N × s Σ , q + ··· Σ , p ... ... ... ... ... ... ... ... Σ q − , ··· Σ q − , q − Σ q − , q −∞ N q − × s Σ q − , q + ··· Σ q − , p Σ q , ··· Σ q , q − Σ qq , Σ qq , Σ qq , Σ q , q + ··· Σ q , p Σ qq , Σ qq , Σ qq , Σ q + , ··· Σ q + , q − Σ q + , q −∞ N q + × s Σ q + , q + ··· Σ q + , p ... ... ... ... ... ... ... ... Σ p , ··· Σ p , q − Σ p , q −∞ N p × s Σ p , q + ··· Σ p , p f i for i ∈ B < q o f i for i ∈ B q o g r for r ∈ J f i for i ∈ B > q | {z } | {z } | {z } | {z } x j for j ∈ B < q x j for j ∈ B q y j for j ∈ J x j for j ∈ B > q Fig. 3.1: Block structure of
ΣΣΣ of the resulting DAE by the block ES method. Thenotation B < q is short for ∪ q − m = B w , and B > q is short for ∪ pm = q + B m .From the description of the conversion in Theorem 3.2, the substitutions (3.12)only occur in equations f i with i ∈ B q . Hence, in the resulting DAE, variables y j for j ∈ J only appear in f i for i ∈ B q and g r for r ∈ J .Considering the block structure of ΣΣΣ in Figure 3.1, we distinguish between thefour cases for a block submatrix
ΣΣΣ m m : (a) m (cid:54) = q and m (cid:54) = q , (b) m (cid:54) = q and m = q , (c) m = q and m (cid:54) = q , and (d) m = m = q .(a) m (cid:54) = q and m (cid:54) = q . In ΣΣΣ m m , equations f i are of indices i ∈ B < q ∪ B > q . As notedin (3.9), the expression substitutions described in (3.12) only take place in f i (cid:48) with i (cid:48) ∈ M ⊆ B q , so do not happen in such blocks ΣΣΣ m m . Hence, each ΣΣΣ m m remains unchanged in ΣΣΣ : ΣΣΣ m m = ΣΣΣ m m for m (cid:54) = q and m (cid:54) = q .(b) m (cid:54) = q and m = q . In ΣΣΣ m q , we include variables y j for j ∈ J as defined in (3.11).By the same arguments as in (a), the expression substitutions do not happen inthese blocks. That is, y j for j ∈ J do not appear in equations f i for i ∈ B < q ∪ B > q .Hence, we can obtain ΣΣΣ m q by concatenating horizontally ΣΣΣ m q with an N m × s matrix of − ∞ ’s: ΣΣΣ m q = (cid:2) ΣΣΣ m q , −−− ∞∞∞ N m × s (cid:3) for m = , . . . , q − , q + , . . . , p .(c) m = q and m (cid:54) = q . In ΣΣΣ qm , we include equations g r for r ∈ J as defined in(3.13). Also, due to the expression substitutions (3.12) occurring in f i with i ∈ M ⊆ B q , σ ( x j , f i ) and σ (cid:0) x j , f i (cid:1) may not be the same for all i ∈ B q and all j = n . Hence, in contrast to cases (a) and (b), there are no obvious connections between
ΣΣΣ m m and ΣΣΣ m m for m = q and m (cid:54) = q .(d) m = m = q . ΣΣΣ qq contains signature entries for equations f i and g r , where i ∈ B q and r ∈ J , in variables x j and y r , where j ∈ B q and r ∈ J . Similar to theresulting signature matrix ΣΣΣ by the basic ES method [13, § ΣΣΣ qq by the blockES method also has a (sub)block structure; c.f (A.11). We shall use it in the proofof Lemma 3.2 in Appendix A.Let Q = ∑ qm = N m denote the total number of equations (or variables) in the first q blocks. Using a valid ( c ; d ) of ΣΣΣ , we construct an offset pair ( (cid:101) c ; (cid:101) d ) defined as (cid:101) c i = c i if i = Qc if i = Q + Q + sc i − s if i = Q + s + n + s , (cid:101) d j = d j if j = Qc if j = Q + Q + sd j − s if j = Q + s + n + s . (3.14)Then we have the following lemma. Since its proof is rather technical, we present itin Appendix A. Lemma 3.2
Let ( (cid:101) c ; (cid:101) d ) be as constructed in (3.14). In the block structure of ΣΣΣ inFigure 3.1, if a position ( i , j ) in ΣΣΣ is below the block diagonal, then (cid:101) d j − (cid:101) c i > σ i j ;otherwise, (cid:101) d j − (cid:101) c i ≥ σ i j . Using this lemma, we can now prove Theorem 3.2.
Proof
Let T be a transversal of ΣΣΣ . Using Lemma 3.2 and (3.14), we deriveVal ( ΣΣΣ ) = ∑ ( i , j ) ∈ T σ i j ≤ ∑ ( i , j ) ∈ T ( (cid:101) d j − (cid:101) c i ) = n + s ∑ j = (cid:101) d j − n + s ∑ i = (cid:101) c i = (cid:32) Q ∑ j = d j + sc + n + s ∑ j = Q + s + d j − s (cid:33) − (cid:32) Q ∑ j = c i + sc + n + s ∑ i = Q + s + c i − s (cid:33) = n ∑ j = d j − n ∑ i = c i = Val ( ΣΣΣ ) . As in the proof of Theorem 3.1, we prove Val ( ΣΣΣ ) < Val ( ΣΣΣ ) by contradiction.Assume Val ( ΣΣΣ ) = ∑ n + sj (cid:101) d j − ∑ n + si (cid:101) c i = Val ( ΣΣΣ ) , which is ≥
0. With Lemma 3.2, thethree conditions in Lemma 2.3 are satisfied. Then it follows from Lemma 2.3 that(a) ( (cid:101) c ; (cid:101) d ) is a valid offset pair of ΣΣΣ ;(b) the Jacobian pattern S , derived from ΣΣΣ and ( (cid:101) c ; (cid:101) d ) , is in the p × p BTF shown inFigure 3.1;(c) T is the union of HVTs T m of all diagonal blocks ΣΣΣ , . . . , ΣΣΣ pp of ΣΣΣ .We can consider block q of the original DAE as a sub-DAE, with signature matrix ΣΣΣ qq and offset pair ( c q ; d q ) —this follows from Lemma 2.4. The conversion describedin Theorem 3.2 can be regarded as an application of the basic ES method to this sub-DAE, given that the ES conditions [13, (4.11)] hold because of (3.10). By [13, Theo-rem 4.2] for the basic ES method, a conversion results in Val ( ΣΣΣ qq ) < Val ( ΣΣΣ qq ) . Also, onversion Methods, Block Triangularization, and Structural Analysis of DAEs 17 since ΣΣΣ mm = ΣΣΣ mm for m (cid:54) = q , Val ( ΣΣΣ mm ) = Val ( ΣΣΣ mm ) . Then a contradiction Val ( ΣΣΣ ) < Val ( ΣΣΣ ) follows by the same derivations as in (3.8). The assumption Val ( ΣΣΣ ) =
Val ( ΣΣΣ ) is hence false, so Val ( ΣΣΣ ) < Val ( ΣΣΣ ) holds. (cid:117)(cid:116) We demonstrate how to apply the block conversion methods on two DAE problemsoriginated from electrical circuit analysis [3]. They are the transistor amplifier andthe ring modulator. We describe them in § § = f = C ( x (cid:48) − x (cid:48) ) + R − ( x − U e ( t )) = f = − C ( x (cid:48) − x (cid:48) ) − R − U b + x (cid:0) R − + R − (cid:1) − ( α − ) g ( x − x ) = f = C x (cid:48) − g ( x − x ) + R − x = f = C ( x (cid:48) − x (cid:48) ) + R − ( x − U b ) + α g ( x − x ) = f = − C ( x (cid:48) − x (cid:48) ) − R − U b + x (cid:0) R − + R − (cid:1) − ( α − ) g ( x − x ) = f = C x (cid:48) − g ( x − x ) + R − x = f = C ( x (cid:48) − x (cid:48) ) + R − ( x − U b ) + α g ( x − x ) = f = − C ( x (cid:48) − x (cid:48) ) + R − x . ΣΣΣ = x x x x x x x x c i f • f • f • f • f • f • f • f • d j J = x (cid:48) x (cid:48) x (cid:48) x (cid:48) x (cid:48) x (cid:48) x (cid:48) x (cid:48) c i f C − C f − C C f C f C − C f − C C f C f C − C f − C C In the equations, g ( y ) = β (cid:0) e y / U F − (cid:1) and U e ( t ) = . ( π t ) ; α , β , U b , U F , R , R k for k = C k for k = J isidentically singular. The fine BTF reveals that the three 2 × J , J , J are identically singular and have a similar structure. Each block receives the sametreatment when a conversion method is applied. LC method.
One can easily find (cid:98) u = [ , ] T ∈ coker ( J ) , coker ( J ) , coker ( J ) . Weperform on each singular block a conversion, and choose to replace the first equation in each such block. block replace by1 f f = f + f f f = f + f f f = f + f The new equations in the resulting DAE are0 = f = R − ( x − U e ( t )) − R − U b + x (cid:0) R − + R − (cid:1) − ( α − ) g ( x − x ) = f = R − ( x − U b ) + α g ( x − x ) − R − U b + x (cid:0) R − + R − (cid:1) − ( α − ) g ( x − x ) = f = R − ( x − U b ) + α g ( x − x ) + R − x . x x x x x x x x c i f • f • f • f • f • f • f • f • d j x (cid:48) x (cid:48) x (cid:48) x (cid:48) x (cid:48) x (cid:48) x (cid:48) x (cid:48) f (cid:48) R − R − ∂ f ∂ x ∂ f ∂ x f − C C f (cid:48) R − R − + R − ∂ f ∂ x ∂ f ∂ x ∂ f ∂ x f − C C f C f (cid:48) R − R − + R − ∂ f ∂ x f − C C f C The SA still reports index 1, and succeeds with a nonzero constant det ( J ) :det ( J ) = C C C C C (cid:0) R − + R − + R − (cid:1) (cid:0) R − + R − + R − (cid:1) (cid:0) R − + R − (cid:1) (cid:54) = . Now Val ( ΣΣΣ ) = < = Val ( ΣΣΣ ) . ES method.
We can take (cid:98) v = [ , ] T ∈ ker ( J ) , ker ( J ) , ker ( J ) . We show how toperform a conversion on block 1; block 3 and block 5 can be treated in the same way.For block 1, we construct the corresponding v = [ , , T ] T . Using (3.9), we have J = J = (cid:8) j | v j (cid:54)≡ (cid:9) = { , } , s = | J | = , M = { , } , and c = . We choose l = ∈ J , introduce for x a new variable y = x ( d − c ) − v v · x ( d − c ) = x (cid:48) − x (cid:48) , and append correspondingly the equation 0 = h = − y + x (cid:48) − x (cid:48) . Then we replace x (cid:48) by y + x (cid:48) in f , f . onversion Methods, Block Triangularization, and Structural Analysis of DAEs 19 After we complete similar conversions on block 3 and block 5, the resulting DAEhas equations f , f and the following equations:0 = f = − C y + R − ( x − U e ( t )) = f = C y − R − U b + x (cid:0) R − + R − (cid:1) − ( α − ) g ( x − x ) = h = − y + x (cid:48) − x (cid:48) = f = − C y + R − ( x − U b ) + α g ( x − x ) = f = C y − R − U b + x (cid:0) R − + R − (cid:1) − ( α − ) g ( x − x ) = h = − y + x (cid:48) − x (cid:48) = f = − C y + R − ( x − U b ) + α g ( x − x ) = f = C y + R − x = h = − y + x (cid:48) − x (cid:48) . x x y x x y x x x y x c i h • f • f • h • f • f • f • h • f • f • f • d j x (cid:48) x (cid:48) y (cid:48) x (cid:48) x (cid:48) y (cid:48) x (cid:48) x (cid:48) x (cid:48) y (cid:48) x (cid:48) h − f (cid:48) R − C f (cid:48) R − − C ∂ f ∂ x ∂ f ∂ x h − f (cid:48) ∂ f ∂ x C ∂ f ∂ x f (cid:48) R − − C ∂ f ∂ x ∂ f ∂ x f C h − f (cid:48) ∂ f ∂ x C ∂ f ∂ x f (cid:48) R − − C f C Here in J , ∂ f / ∂ x = R − + R − and ∂ f / ∂ x = R − + R − . The SA succeeds witha nonzero constant det ( J ) and Val ( ΣΣΣ ) = < = Val ( ΣΣΣ ) .4.2 Ring modulator.We study the ring modulator problem from [3]. When C s (cid:54) =
0, it is a stiff ODE systemof 15 nonlinear equations. Setting C s = consists of 11 differential and 4 algebraic equations:0 = f = − x (cid:48) + C − (cid:0) x − . x + . x + x − R − x (cid:1) = f = − x (cid:48) + C − (cid:0) x − . x + . x + x − R − x (cid:1) = f = x − q ( U D ) + q ( U D ) = f = − x + q ( U D ) − q ( U D ) = f = x + q ( U D ) − q ( U D ) = f = − x − q ( U D ) + q ( U D ) = f = − x (cid:48) + C − p (cid:0) − R − p x + q ( U D ) + q ( U D ) − q ( U D ) − q ( U D ) (cid:1) = f = − x (cid:48) + − L − h x = f = − x (cid:48) + − L − h x = f = − x (cid:48) + L − s ( . x − x − R g x ) = f = − x (cid:48) + L − s ( − . x + x − R g x ) = f = − x (cid:48) + L − s ( . x − x − R g x ) = f = − x (cid:48) + L − s ( − . x + x − R g x ) = f = − x (cid:48) + L − s ( − x + U in ( t ) − ( R i + R g ) x ) = f = − x (cid:48) + L − s ( − x − ( R c + R g ) x ) . The functions are U D = x − x − x − U in ( t ) q ( U ) = γ ( e δ U − ) U D = − x + x − x − U in ( t ) U in ( t ) = . ( π t ) U D = x + x + x + U in ( t ) U in ( t ) = ( π t ) U D = − x − x + x + U in ( t ) . onversion Methods, Block Triangularization, and Structural Analysis of DAEs 21 We refer the reader to [3] for the nonzero constants C , C p , R , R p , R c , γ , L h , L s , L s , L s , R g , R g , R g , R i , and δ . ΣΣΣ = x x x x x x x x x x x x x x x c i f • f • f • f • f • f • f • f • f • f • f • f • f • f • f • d j Each block of size 1 has a nonsingular sub-Jacobian −
1. Block 8 has an identicallysingular sub-Jacobian J = x x x x f − s − s s − s f − s − s − s s f s − s − s − s f − s s − s − s , where s i = γδ e δ U Di . This is a nonlinear block, since variables x , x , x , x do not occur jointly linearly inequations f , f , f , f . One can also see these variables appear in J . LC method.
We find a constant vector (cid:98) u = [ , − , , − ] T ∈ coker ( J ) , which satis-fies the block LC condition (3.3). Then u = [ T , , − , , − , T ] T . We use (3.2) toderive I = (cid:8) i | u i (cid:54)≡ (cid:9) = { , , , } , c = , and L = L = { , , , } . The row indices in L correspond to the equations f , f , f , f . We can pick any oneof them and replace it by f = u f + u f + u f + u f = f − f + f − f = x + x + x + x . We choose f and replace it by f = f . The resulting DAE has the following ΣΣΣ withVal ( ΣΣΣ ) = < = Val ( ΣΣΣ ) . ΣΣΣ = x x x x x x x x x x x x x x x c i f • f • f • f • f • f • f • f • f • f • f • f • f • f • f • d j Again, each 1 × ∂ f i / ∂ x (cid:48) i = − i =
1, 2, 7,8, 9, 14, 15. The sub-Jacobian of block 4 in the resulting DAE is J = x x x x x (cid:48) x (cid:48) x (cid:48) x (cid:48) f − L − s − f s − s − s − s f − s − s − s s f − s s − s − s f (cid:48) f L − s − f − L − s − f L − s − , whose determinant is det ( J ) = s s s s ( s − + s − + s − + s − )( L − s + L − s ) . TheSA succeeds at any point where det ( J ) (cid:54) =
0, and the DAE is of index 2.
ES method.
Find (cid:98) v = [ − , , − , ] T ∈ ker ( J ) . Then v = [ T , − , , − , , T ] T . Weuse (3.9) to derive J = J = (cid:8) j | v j (cid:54)≡ (cid:9) = { , , , } , s = | J | = , M = J , and c = . We choose column index l = ∈ J in the permuted ΣΣΣ . The variable of this col-umn is x . The other variables in block 8 are x , x , x , so we introduce for them,respectively, y = x − ( v / v ) · x , y = x − ( v / v ) · x , and y = x − ( v / v ) · x . onversion Methods, Block Triangularization, and Structural Analysis of DAEs 23 Then we append the equations corresponding to these variables0 = g = − y + x + x , = g = − y + x − x , and 0 = g = − y + x + x . The equations in block 8 are f , f , f , f . In these equations, we perform the follow-ing substitutions. replace by in x y − x f , f , f x y + x f , f , f x y − x f , f , f The resulting index-2 DAE is of size 18. (We do not display the results of SA here.) Ithas Val ( ΣΣΣ ) = < = Val ( ΣΣΣ ) and det ( J ) = − s s s s ( s − + s − + s − + s − )( L − s + L − s ) . The largest fine block is of size 12, and the other six fine blocks are of size 1.The SA succeeds at any point where det ( J ) (cid:54) = We combined block triangularization with the LC and ES conversion methods for im-proving the Σ -method. When J is identically singular and the DAE has a nontrivialBTF, we can locate each diagonal block whose corresponding sub-Jacobian is identi-cally singular, and perform a conversion on it. We base this strategy on the view thateach diagonal block can be regarded as a sub-DAE, while formulas contributing tothe off diagonal blocks are regarded as driving terms.Compared with the basic conversion methods that work on the whole DAE, theblock methods only work on singular blocks, which are usually smaller than the DAEitself. Hence the block methods require fewer symbolic computations, and can gen-erally find a useful conversion for reducing Val ( ΣΣΣ ) more efficiently. As in the basiccase, a conversion applied on a singular block guarantees (a) a strict decrease in thevalue of the (whole) signature matrix, and (b) the equivalence between the originalDAE and the resulting one. The rationale for choosing a desirable conversion methodis in [13, Table 4.1].We combine M ATLAB ’s Symbolic Math Toolbox [14] with our structural analysissoftware
DAESA [6,10], and have built a prototype code that automates the conversionprocess. We aim to incorporate them in a future version of
DAESA .With our prototype code, we have applied our methods on numerous DAEs onwhich the Σ -method fails. They are either arbitrarily constructed to be SA-failurecases for our investigations, or borrowed from the existing literature. Our conversionmethods succeed in fixing all these solvable DAEs. We believe that our assumptionsand conditions are reasonable for practical problems, and that these methods can helpmake the Σ -method more reliable.We end these two articles with our main conjecture regarding SA’s failure. In allour experiments, when we successfully fix the failure using our conversion methods,the value of a signature matrix always decreases. As Pryce points out in [8], thesolvability of a DAE lies within its inherent nature, not the way it is formulated oranalyzed. Hence we conjecture that a DAE formulation friendly to SA should have a reasonable but never overestimated Val ( ΣΣΣ ) that can be interpreted as the number ofdegrees of freedom (DOF) of the underlying problem. In other words, a DAE shouldnot be formulated to exhibit more DOF than the underlying problem has. However,based on our current knowledge, it appears difficult to show why overestimating DOFcan lead to an identically singular System Jacobian. A Proof of Lemma 3.2.
For
ΣΣΣ = ( σ ij ) in the block structure in Figure 3.1, we write the block sizes in the array N = ( N , N ,..., N q − , N q + s , N q + ,..., N p ) , and also write the block sizes of ΣΣΣ in the array N = ( N , N ,..., N q − , N q , N q + ,..., N p ) . Let blockOf ( i ) denote the block number of a row or column index i in ΣΣΣ . From the construction of N and N , it is notdifficult to show that blockOf ( j ) < q ⇔ ≤ j ≤ q − ∑ w = N w ⇔ blockOf ( j ) < q and (A.1)blockOf ( j + s ) > q ⇔ q ∑ w = N w + s + ≤ j + s ⇔ blockOf ( j ) > q . (A.2)From the construction of ( (cid:101) c ; (cid:101) d ) in (3.14), each variable x j for j = n has the same “variable offset”in ΣΣΣ as x j has in ΣΣΣ . Also, each equation f i for i = n has the same “equation offset” in ΣΣΣ as f i has in ΣΣΣ .Quotation marks are used here because ( (cid:101) c ; (cid:101) d ) is not a valid offset pair of ΣΣΣ ; this vector pair is merely usedfor proving Val ( ΣΣΣ ) < Val ( ΣΣΣ ) in Theorem 3.2.We aim to show that (cid:101) d j − (cid:101) c i (cid:40) > σ ij if blockOf ( j ) < blockOf ( i ) ≥ σ ij if blockOf ( j ) ≥ blockOf ( i ) . (A.3)For the block structure of ΣΣΣ in Figure 3.1, we have shown on page 15 that
ΣΣΣ m m = ΣΣΣ m m if m (cid:54) = q and m (cid:54) = q , and that ΣΣΣ m m = (cid:2) ΣΣΣ m q , −−− ∞∞∞ N m × s (cid:3) if m (cid:54) = q and m = q . Hence, provided m (cid:54) = q , ΣΣΣ m m is below [resp. above] the block diagonal of ΣΣΣ , if
ΣΣΣ m m is below [resp. above] the block diagonal of ΣΣΣ .By (2.4), the inequalities in (A.3) hold for i with blockOf ( i ) (cid:54) = q .What remains to show is the inequalities in (A.3) for i with blockOf ( i ) = q . These inequalities are forthe signature entries in ΣΣΣ qm , the blocks affected by the expression substitutions. We consider three casesfor ΣΣΣ qm : it is (a) below the block diagonal, with m < q , (b) above the block diagonal, with m > q , or (c)the diagonal block ΣΣΣ qq , with m = q .(a) ΣΣΣ qm with m < q . An entry ( i , j ) in this block satisfies blockOf ( j ) < blockOf ( i ) = q . By (A.1),blockOf ( j ) < q and hence j / ∈ B q .Recall from (3.12) that, in each f i with i ∈ M ⊆ B q , we substitute (cid:16) y r + v r v l · x ( d l − c ) l (cid:17) ( c − c i ) for each x ( σ ir ) r with d r − c i = σ ir and r ∈ J \ (cid:8) l (cid:9) ⊂ B q . For a j / ∈ B q ⊃ J \ (cid:8) l (cid:9) , the corresponding derivatives x ( d j − c i ) j are not replaced in the ES conversion, and for r ∈ J \ (cid:8) l (cid:9) (so j , r , l are distinct), σ (cid:18) x j , (cid:16) y r + v r v l · x ( d l − c ) l (cid:17) ( c − c i ) (cid:19) = σ (cid:18) x j , (cid:16) v r v l (cid:17) ( c − c i ) (cid:19) ≤ σ (cid:16) x j , v ( c − c i ) (cid:17) = σ (cid:0) x j , v (cid:1) + ( c − c i ) . (A.4)By (3.10), σ (cid:0) x j , v (cid:1) < d j − c . Using (2.4) and (A.4), we derive σ (cid:0) x j , f i (cid:1) ≤ max (cid:26) σ (cid:0) x j , f i (cid:1) , max r ∈ J \{ l } σ (cid:18) x j , (cid:16) y r + v r v l · x ( d l − c ) l (cid:17) ( c − c i ) (cid:19)(cid:27) ≤ max (cid:8) σ ij , σ (cid:0) x j , v (cid:1) + ( c − c i ) (cid:9) < max (cid:8) d j − c i , ( d j − c ) + ( c − c i ) (cid:9) = d j − c i for i ∈ M ⊆ B q . (A.5)onversion Methods, Block Triangularization, and Structural Analysis of DAEs 25From the ES conversion described in Theorem 3.2, we have σ (cid:0) x j , f i (cid:1) = σ (cid:0) x j , f i (cid:1) < d j − c i for i ∈ B q \ M and (A.6) σ (cid:0) x j , g r (cid:1) ≤ σ (cid:0) x j , v (cid:1) < d j − c for r ∈ J . (A.7)Since blocks ΣΣΣ qm with m < q contain signature entries σ ij for equations f i and g r , where i ∈ B q and r ∈ J , in variables x j with blockOf ( j ) < q , by taking together the inequalities in (A.5)-(A.7), we have σ ij < (cid:40) d j − c i if blockOf ( j ) < q and i ∈ B q d j − c if blockOf ( j ) < q and i ∈ Q + Q + s ;recall Q = ∑ qw = N w . Using (A.1) and the construction of ( (cid:101) c ; (cid:101) d ) in (3.14), we have σ ij < (cid:101) d j − (cid:101) c i for blockOf ( j ) < blockOf ( i ) = q . ( q is the block number of both the original and enlarged diagonal blocks.)(b) ΣΣΣ qm with m > q . An entry ( i , j + s ) in this block satisfies blockOf ( j + s ) > blockOf ( i ) = q . By(A.2), blockOf ( j ) > q and hence j / ∈ B q ⊃ J \ (cid:8) l (cid:9) . By the same arguments as in (a), the correspondingderivatives x ( d j − c i ) j are not replaced in the ES conversion.By (3.10), σ (cid:0) x j , v (cid:1) ≤ d j − c . Then by the same derivations as (A.5–A.7) in (a), we have σ (cid:0) x j , f i (cid:1) ≤ d j − c i for i ∈ M ⊆ B q , (A.8) σ (cid:0) x j , f i (cid:1) = σ (cid:0) x j , f i (cid:1) ≤ d j − c i for i ∈ B q \ M , and (A.9) σ (cid:0) x j , g r (cid:1) ≤ σ (cid:0) x j , v (cid:1) ≤ d j − c for r ∈ J . (A.10)Since blocks ΣΣΣ qm with m > q contain signature entries σ i , j + s for equations f i and g r , where i ∈ B q and r ∈ J , in variables x j with blockOf ( j ) > q , the inequalities (A.8–A.10) yield σ i , j + s ≤ (cid:40) d j − c i if blockOf ( j ) > q and i ∈ B q d j − c if blockOf ( j ) > q and i ∈ Q + Q + s . Using again (A.2) and ( (cid:101) c ; (cid:101) d ) in (3.14), we have σ i , j + s ≤ (cid:101) d j + s − (cid:101) c i for blockOf ( j + s ) > blockOf ( i ) = q ,with j = Q + n . We can rewrite this inequality as σ ij ≤ (cid:101) d j − (cid:101) c i for blockOf ( j ) > blockOf ( i ) = q , with j = Q + + s : n + s .(c) ΣΣΣ qm is ΣΣΣ qq , with m = q . An entry ( i , j ) in ΣΣΣ qq satisfies blockOf ( j ) = blockOf ( i ) = q . We viewblock q of the original DAE as a sub-DAE, with a signature matrix ΣΣΣ qq of size N q and an offset pair ( c q ; d q ) . Given that the ES conditions are satisfied by (3.10), performing the ES conversion as describedin Theorem 3.2 is equivalent to applying the basic ES method to this sub-DAE. After a conversion, theresulting enlarged signature matrix ΣΣΣ qq of size N q + s has the form ΣΣΣ qq = (cid:34) ΣΣΣ qq , ΣΣΣ qq , ΣΣΣ qq , ΣΣΣ qq , ΣΣΣ qq , ΣΣΣ qq , (cid:35) ; (A.11)cf. Figure 3.1 and the details of the basic ES method in [13]. The two block rows of ΣΣΣ qq correspond to f i for i ∈ B q and g j for j ∈ J , respectively. The three block columns of ΣΣΣ qq correspond to x j for j ∈ J , x j for j ∈ B q \ J , and y j for j ∈ J , respectively. If we apply the same arguments in the proof of [13, Lemma 4.4]for the basic ES method, then we have (cid:101) d j − (cid:101) c i ≥ σ ij for all entries in ΣΣΣ qq . (cid:117)(cid:116) Acknowledgements
The authors acknowledge with thanks the financial support for this research: GT issupported in part by the McMaster Centre for Software Certification through the Ontario Research Fund,Canada, NSN is supported in part by the Natural Sciences and Engineering Research Council of Canada,and JDP is supported in part by the Leverhulme Trust, the UK.6 Guangning Tan et al.
References
1. Campbell, S.L., Griepentrog, E.: Solvability of general differential-algebraic equations. SIAM Journalon Scientific Computing (2), 257–270 (1995)2. Duff, I., Erisman, A., Reid, J.: Direct Methods for Sparse Matrices. Oxford Science Publications.Clarendon Press, Oxford (1986)3. Mazzia, F., Iavernaro, F.: Test set for initial value problem solvers. Tech. Rep. 40, Department ofMathematics, University of Bari, Italy (2003). http://pitagora.dm.uniba.it/~testset/
4. McKenzie, R.: Structural analysis based dummy derivative selection for differential-algebraic equa-tions. Tech. rep., Cardiff University, UK (2015). Submitted to BIT Numerical Mathematics, Oct2015.5. McKenzie, R.: Reducing the index of differential-algebraic equations by exploiting underlying struc-tures. PhD Thesis, School of Mathematics, Cardiff University, Senghennydd Road, Cardiff CF244AG, UK (2016)6. Nedialkov, N.S., Pryce, J.D., Tan, G.: Algorithm 948: DAESA—a Matlab tool for structural analysisof differential-algebraic equations: Software. ACM Trans. Math. Softw. (2), 12:1–12:14 (2015)7. Nedialkov, N.S., Tan, G., Pryce, J.D.: Exploiting fine block triangularization and quasilinearity indifferential-algebraic equation systems (2016). McMaster University, Cardiff University. In prepara-tion8. Pryce, J.D.: Solving high-index DAEs by Taylor Series. Numerical Algorithms , 195–211 (1998)9. Pryce, J.D.: A simple structural analysis method for DAEs. BIT Numerical Mathematics (2), 364–394 (2001)10. Pryce, J.D., Nedialkov, N.S., Tan, G.: DAESA—a Matlab tool for structural analysis of differential-algebraic equations: Theory. ACM Trans. Math. Softw. (2), 9:1–9:20 (2015)11. Pryce, J.D., Nedialkov, N.S., Tan, G.: Fine block triangular structure of DAEs and its uses (2016).Cardiff University, McMaster University. In preparation12. Tan, G., Nedialkov, N., Pryce, J.: Symbolic-numeric methods for improving structural analysis ofdifferential-algebraic equation systems. Tech. rep., Department of Computing and Software, Mc-Master University, 1280 Main St. W., Hamilton, ON, L8S4L8, Canada (2015). CAS-15-07-NN, 84pages13. Tan, G., Nedialkov, N.S., Pryce, J.D.: Conversion methods for improving structural analysis ofdifferential-algebraic equation systems. BIT Numerical Mathematics (2016). Submitted14. The MathWorks, Inc.: Matlab Symbolic Math Toolbox (2016).(2), 9:1–9:20 (2015)11. Pryce, J.D., Nedialkov, N.S., Tan, G.: Fine block triangular structure of DAEs and its uses (2016).Cardiff University, McMaster University. In preparation12. Tan, G., Nedialkov, N., Pryce, J.: Symbolic-numeric methods for improving structural analysis ofdifferential-algebraic equation systems. Tech. rep., Department of Computing and Software, Mc-Master University, 1280 Main St. W., Hamilton, ON, L8S4L8, Canada (2015). CAS-15-07-NN, 84pages13. Tan, G., Nedialkov, N.S., Pryce, J.D.: Conversion methods for improving structural analysis ofdifferential-algebraic equation systems. BIT Numerical Mathematics (2016). Submitted14. The MathWorks, Inc.: Matlab Symbolic Math Toolbox (2016).