Model-order reduction of coupled DAE systems via ε-embedding technique and Krylov subspace method

Available online at www.sciencedirect.com

Journal of the Franklin Institute 349 (2012) 3027–3045 www.elsevier.com/locate/jfranklin

Model-order reduction of coupled DAE systems via e-embedding technique and Krylov subspace method Yao-Lin Jiangn, Chun-Yue Chen, Hai-Bao Chen Department of Mathematical Sciences, School of Science, Xi’an Jiaotong University, Xi’an, Shaanxi 710049, China Received 29 January 2012; received in revised form 10 September 2012; accepted 11 September 2012 Available online 23 September 2012

Abstract For a class of coupled systems which include differential-algebraic equation (DAE) subsystems, we discuss model-order reduction (MOR) for them via e-embedding technique and Krylov subspace method. First we change the coupled system into a closed-loop form. Then, both one-sided and twosided Arnoldi algorithms are used to reduce the closed-loop system. To preserve the interconnected structure, each subsystems can be reduced by the two same algorithms. Next, we show some facts on the error, stability, and passivity for the resulted systems. Finally, we demonstrate the effectiveness of our approach in two examples. & 2012 The Franklin Institute. Published by Elsevier Ltd. All rights reserved.

1. Introduction Coupled systems that consist of ordinary differential equations, differential-algebraic equations and partial differential equations can be naturally obtained by modelling and simulation of complex physical and technical processes. Such systems arise in many practical applications such as very large system integrated (VLSI) chip designs, elastic multibody systems, and micro-electro-mechanical systems (MEMS), see [1–4]. By using model-order reduction techniques, a large-scale system is changed into a system of lower order that has nearly the same response characteristics. Both in the area of control theory and in numerical analysis, a wealth of MOR techniques have been developed over the decades. Pade approximation technique, Balanced truncation, Krylov subspace method, and proper orthogonal decomposition have been extensively used to reduce linear n

Corresponding author. E-mail address: [email protected] (Y.-L. Jiang).

0016-0032/$32.00 & 2012 The Franklin Institute. Published by Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.jfranklin.2012.09.004

Y.-L. Jiang et al. / Journal of the Franklin Institute 349 (2012) 3027–3045

3028

ODE systems and coupled systems, see [5–8]. In recent years, the MOR techniques for DAE systems have been received increased attention, see [9,10]. Moreover, there are also some methods that are directly used to solve DAE systems such as waveform relaxation, see [11–13]. In this paper, we consider a system of k coupled linear time-invariant generalized state space subsystems in the following form: Ei

dxi ðtÞ ¼ Ai xi ðtÞ þ Bi ui ðtÞ, dt

yi ðtÞ ¼ Ci xi ðtÞ,

ð1:1Þ

with the initial condition xi ð0Þ ¼ xi0 , where Ei ,Ai 2 Rni ni , Bi 2 Rni pi , Ci 2 Rmi ni , xi ðtÞ 2 Rni is the internal descriptor vector, ui ðtÞ 2 Rpi is the internal input variable, and yi ðtÞ 2 Rmi is the internal output variable. These subsystems are coupled through the input–output relations ( ui ðtÞ ¼ Ki1 y1 ðtÞ þ    þ Kik yk ðtÞ þ Gi uðtÞ, ð1:2Þ yðtÞ ¼ R1 y1 ðtÞ þ    þ Rk yk ðtÞ, i ¼ 1,2, . . . ,k, where Kil 2 Rpi ml ðl ¼ 1,2, . . . ,kÞ, Gi 2 Rpi p , Ri 2 Rmmi , uðtÞ 2 Rp is the external input variable, and yðtÞ 2 Rm is the external output variable. We suppose that coupled system (1.1)–(1.2) composed of k1 DAE subsystems and k2 ODE subsystems, where k1 þ k2 ¼ k. Without loss of generality, assume that the first k1 subsystems are DAE systems. According to [14], the matrix pencil (E,A) is called regular if and only if there exists some l 2 C such that detðlEAÞa0. When the pencil ðEi ,Ai Þ is regular, the matrices Ei and Ai ^ i ,In^ i2 g, can be changed into the block form Pi Ei Qi ¼ diagfIn^ i1 , N^ i g and Pi Ai Qi ¼ diagfM ni ni where n^ i1 þ n^ i2 ¼ ni , Pi ,Qi 2 R are both nonsingular, In^ i1 and In^ i2 are the identity matrices of orders n^ i1 and n^ i2 , respectively, N^ i 2 Rn^ i2 n^ i2 is a nilpotent matrix with index of nilpotency ni , which corresponds to the eigenvalue at infinity of the matrix pencil ðEi ,Ai Þ, ^ i 2 Rn^ i1 n^ i1 is the Jordan block corresponding to the finite eigenvalues of ðEi ,Ai Þ. and M This block decomposition is called the Weierstrass–Kronecker canonical form of ðEi ,Ai Þ [15]. In this paper, we discuss the case when the matrices Ei and Ai have the form Ei ¼ diagfIni1 ,Ni g,

Ai ¼ diagfMi ,Ini2 g,

ni1 þ ni2 ¼ ni , i ¼ 1,2, . . . ,k1 :

In this case, the vector xi(t) can be partitioned as xi ðtÞ ¼ ½xTi1 ðtÞ xTi2 ðtÞT , where xi1 ðtÞ 2 Rni1 and xi2 ðtÞ 2 Rni2 . Furthermore, we have xi ð0Þ ¼ ½xTi1 ð0Þ xTi2 ð0ÞT . The paper is organized as follows. In Section 2, we report two methods for reducedorder modelling of coupled systems. In Section 3, model-order reduction of the subsystems is discussed in detail. Also the structure-preserving MOR method is applied to the subsystems. The error, stability, and passivity for our MOR approach in Sections 2 and 3 are considered in Section 4. In Section 5, several experiments are used to verify the effectiveness of our approach. Finally, we show some conclusions in Section 6.

2. Model-order reduction of the closed-loop system In this section, we first consider all subsystems together in the closed-loop form. Then, our MOR approach which is firstly introduced in [16] is applied to the closed-loop system. Moreover, a structure-preserving method is also used to reduce the closed-loop system.

Y.-L. Jiang et al. / Journal of the Franklin Institute 349 (2012) 3027–3045

3029

2.1. e-embedding of the closed-loop system Let n ¼ n1 þ n2 þ    þ nk , p0 ¼ p1 þ p2 þ    þ pk , and m0 ¼ m1 þ m2 þ    þ mk . Introduce some notations E ¼ diagfE1 ,E2 , . . . ,Ek g 2 Rnn , A ¼ diagfA1 ,A2 , . . . ,Ak g 2 Rnn , B ¼ diagfB1 ,B2 , . . . ,Bk g 2 Rnp0 , C ¼ diagfC1 ,C2 , . . . ,Ck g 2 Rm0 n , 2 T 3T 2 3 2 3 R1 K11 K12    K1k G1 6 T7 6 K21 K22    K2k 7 6 G2 7 R 6 2 7 6 7 6 7 7 2 Rmm0 : K ¼6 7 2 Rp0 m0 , G ¼ 6 7 2 Rp0 p , R ¼ 6 6 7 ^ ^ & ^ 5 4 ^ 4 ^ 5 4 5 RTk Kk1 Kk2    Kkk Gk Then, we can rewrite coupled system (1.1)–(1.2) in the closed-loop form dxðtÞ ¼ AxðtÞ þ BuðtÞ, yðtÞ ¼ CxðtÞ, ð2:1Þ dt where E ¼ E 2 Rnn , A ¼ A þ BKC 2 Rnn , B ¼ BG 2 Rnp , C ¼ RC 2 Rmn , and xðtÞ ¼ ½xT1 ðtÞ xT2 ðtÞ    xTk ðtÞT 2 Rn is the descriptor vector. Its initial condition is xð0Þ ¼ ½xT10 xT20    xTk0 T . Let 3 2 0 0 0   0 0   0 In11 7 6 0 In21 0   0 0   0 7 6 0 7 6 6           7 7 6 7 6 0 0 0    0 Ink1 1 0   0 7 6 0 7 6 6 0 0 0 0  0 0 0 Ink1 þ1 0 ^ 7 7 6 7 6 ^ ^ ^ ^ & ^ ^ ^ & & 0 7 2 Rnn : 6 P¼6 7 6 0 0 0 0  0 0 0  0 Ink 7 7 6 7 6 6 0 In12 0 0   0 0   0 7 7 6 6 0 0 0 In22 0  0 0   0 7 7 6 7 6 6           7 5 4 0 0 0 0   0 Ink1 2 0  0 E

^ ^ ¼ PxðtÞ and insert into Eq. (2.1). We change the descriptor vector x(t) into xðtÞ as xðtÞ Then, we get an equivalent DAE system as follows: ^ d xðtÞ ^ xðtÞ ^ ^ þ BuðtÞ, E^ ¼A dt

^ ¼ C^ xðtÞ, ^ yðtÞ

ð2:2Þ

^ ¼ yðtÞ and where yðtÞ E^ ¼ PEP1 ¼ diagfIn11 , . . . ,Ink1 1 ,Ek1 þ1 , . . . ,Ek ,N1 , . . . ,Nk1 g 2 Rnn , ^ ¼ PAP1 2 Rnn , A

B^ ¼ PB 2 Rnp ,

C^ ¼ CP1 2 Rmn :

Its initial condition is of the form ^ xð0Þ ¼ ½xT11 ð0Þ    xTk1 1 ð0Þ xTk1 þ1 ð0Þ    xTk ð0Þ xT12 ð0Þ    xTk1 2 ð0ÞT :

3030

Y.-L. Jiang et al. / Journal of the Franklin Institute 349 (2012) 3027–3045

We set n^ 1 ¼ n11 þ n21 þ    þ nk1 1 , n^ 2 ¼ n12 þ n22 þ    þ nk1 2 , and N ¼ diagf N1 , . . . ,Nk1 g. Applying an e-embedding method to DAE system (2.2) results in a system, given by ^ d xðtÞ ^ xðtÞ ^ ^ þ BuðtÞ, ^ ¼ C^ xðtÞ, ^ ¼A yðtÞ ð2:3Þ E^ e dt where E^ e ¼ diagfIn^ 1 ,Ek1 þ1 , . . . ,Ek ,eIn^ 2 þ Ng and the parameter e such that e40, e51. It has the same initial condition as system (2.2). It is not difficult to see that E^ e is nonsingular, when eojlN j holds for each nonzero eigenvalue lN of the matrix N. In this case, Eq. (2.3) is an ODE system. The transfer functions of systems (2.2) and (2.3) are in the following form: ^ 1 B, ^ E^ AÞ ^ ^ ðsÞ ¼ Cðs H

^ 1 B, ^ E^ e AÞ ^ ^ e ðsÞ ¼ Cðs H

^ ðsÞH ^ e ðsÞJ2 as follows: respectively. Then, we can bound JH ^ 1 J2 JðsE^ e AÞ ^ 1 J2 ^ 2 JBJ ^ 2 JðsE^ AÞ ^ ðsÞH ^ e ðsÞJ2 rejsjJCJ JH ^ 1 J2 ^ 2 JBJ ^ 2 JðsE^ AÞ ejsjJCJ 2 r : ^ 1 J2 1ejsjJðsE^ AÞ Since the equality PT ¼ P1 holds, it is easy to see that T T ^ T ðsE^ AÞ ^ 1 ¼ PðsEAÞT ðsEAÞ1 P1 : C^ C^ ¼ PCT CP1 , B^ B^ ¼ BT B, ðsE^ AÞ ^ 1 J2 ¼ JðsEAÞ1 J2 . Therefore, ^ 2 ¼ JBJ2 , and JðsE^ AÞ ^ 2 ¼ JCJ2 , JBJ This implies that JCJ we have 1 2 ^ ðsÞH ^ e ðsÞJ2 r ejsjJCJ2 JBJ2 JðsEAÞ J2 : JH 1ejsjJðsEAÞ1 J2

ð2:4Þ

Thus, system (2.3) is a good approximation of closed-loop system (2.1) and DAE system ^ ðsÞ ¼ H ^ e ðsÞ. (2.2), if the parameter e is very small. Clearly, when e ¼ 0, we have H 2.2. Model-order reduction of the closed-loop system In this section, one-sided and two-sided Arnoldi algorithms are used to reduce embedding system (2.3), respectively. Furthermore, we also discuss the structurepreserving methods in details. Throughout the following paper, assume that the matrix ^ is nonsingular. A 1. One-sided Arnoldi algorithm for the closed-loop system. Consider the Krylov subspace ^ 1 E^ e , A ^ 1 BÞ, ^ the block Arnoldi algorithm is used to construct a column-orthonormal Kr ðA nq ^ into xðtÞ ~ as matrix Ve 2 R , where q ¼ rp and q5n. We change the descriptor vector xðtÞ ^xðtÞ  Ve xðtÞ. ~ Substituting this approximation into system (2.3), we get a reduced-order system as follows: ~ d xðtÞ ~ e xðtÞ ~ ¼ C~ e xðtÞ, ~ þ B~ e uðtÞ, yðtÞ ~ ¼A ð2:5Þ E~ e dt ^ e 2 Rqq , B~ e ¼ V T B^ 2 Rqp , C~ e ¼ CV ~ e ¼ V T AV ^ e 2 Rmq , where E~ e ¼ VeT E^ e Ve 2 Rqq , A e e q m ~ 2 R is the output variable. Its initial condition ~ 2 R is the descriptor vector, and yðtÞ xðtÞ T ~ e Þ1 B~ e . ~ e ðsÞ ¼ C~ e ðsE~ e A ~ ^ is xð0Þ ¼ Ve xð0Þ. The transfer function has the form H

Y.-L. Jiang et al. / Journal of the Franklin Institute 349 (2012) 3027–3045

3031

The matrix Ve can be expanded as Ve ¼ V0 þ V1 e þ   , where Vl 2 Rnq , l ¼ 0,1, . . .. Since Ve is the column-orthonormal matrix, then V0 ¼ Ve je ¼ 0 is also a columnorthonormal matrix. Specially, if we choose the parameter e ¼ 0 in system (2.5), we get ^ 1 E^ e , A ^ 1 BÞ ^ ¼ colspan a reduced-order system of original system (2.2). Note that Kr ðA 1

1

^ E^ , A ^ BÞ ^ ¼ colspanfV0 g. Thus, this special reduced-order fVe g, clearly we have Kr ðA system is equivalent to the system that is constructed by the direct one-sided Arnoldi MOR algorithm. In general, the coefficient matrix E~ e in Eq. (2.5) does not have the same special block structure of the matrix E^ e in Eq. (2.3). Next based on Ve , we give a structure-preserving MOR algorithm [18]. The matrix Ve can be partitioned as T T T Ve2    Ve,kþk T , Ve ¼ ½Ve1 1

ð2:6Þ

where Vei 2 Rni1 q ði ¼ 1,2, . . . ,k1 Þ; Vei 2 Rni q ði ¼ k1 þ 1,k1 þ 2, . . . ,kÞ; Vei 2 Rnik,2 q ði ¼ k þ 1,k þ 2, . . . ,k þ k1 Þ. Let V^ e ¼ diagfVe1 ,Ve2 , . . . ,Ve,kþk1 g 2 Rnðkþk1 Þq . Applying the Gram–Schmidt (GS) scheme [19] to the matrix V^ e derives a column-orthonormal matrix, given by V~ e ¼ diagfV~ e1 , V~ e2 , . . . , V~ e,kþk1 g 2 Rnq~ , where V~ ei 2 Rni1 q~ i ði ¼ 1,2, . . . ,k1 Þ; V~ ei 2 Rni q~ i ði ¼ k1 þ 1,k1 þ 2, . . . ,kÞ; V~ ei 2 Rnik,2 q~ i ði ¼ k þ 1,k þ 2, . . . ,k þ k1 Þ. Here Pkþk1 ~ q~ i rq, qZq ~ ~ ~ i ¼ q, and q5n. i¼1 q ^  V~ e x~ q~ ðtÞ, and inserting into Eq. (2.3) results in a Taking the approximation xðtÞ reduced-order system in the form d x~ q~ ðtÞ ~ eq~ x~ q~ ðtÞ þ B~ eq~ uðtÞ, y~ ðtÞ ¼ C~ eq~ x~ q~ ðtÞ, ¼A E~ eq~ q~ dt T ^ with the initial condition x~ q~ ð0Þ ¼ V~ e xð0Þ, where x~ q~ ðtÞ 2 Rq~ and

ð2:7Þ

T T T E~ eq~ ¼ V~ e E^ e V~ e ¼ diagfIq~ 0 , V~ e,k1 þ1 Ek1 þ1 V~ e,k1 þ1 , . . . , V~ e,k Ek V~ e,k , T T ~ q~ eIq~ þ V~ N1 V~ e,kþ1 , . . . ,eIq~ þ V~ Nk1 V~ e,kþk1 g 2 Rq , kþ1

e,kþ1

~ q~ q ~ eq~ ¼ V~ T A ^ ~ A , e Ve 2 R

kþk1

e,kþk1

T ~ B~ eq~ ¼ V~ e B^ 2 Rqp ,

C~ eq~ ¼ C^ V~ e 2 Rmq~ ,

q~ 0 ¼

k1 X

q~ i :

i¼1

~ eq~ Þ1 B~ eq~ . Clearly, the matrix E~ eq~ ~ eq~ ðsÞ ¼ C~ eq~ ðsE~ eq~ A The transfer function is given by H preserves the block diagonal structure of the coefficient matrix E^ e in system (2.3). 2. Two-sided Arnoldi algorithm for the closed-loopT system.TWe apply the block Arnoldi ^ E^ T , A ^ C^ T Þ, and obtain a columnalgorithm to the r^ th block Krylov subspace Kr^ ðA e orthonormal matrix We 2 Rnq , where q ¼ r^ m ¼ rp. Then from the projection matrices We and Ve , a reduced-order system can be constructed as d x~ We ðtÞ ~ W x~ W ðtÞ þ B~ W uðtÞ, y~ ðtÞ ¼ C~ W x~ W ðtÞ, ¼A ð2:8Þ E~ We We e e e e e dt ^ e , B~ W ¼ W T B, ~ W ¼ W T AV ^ e . It has the same initial ^ and C~ W ¼ CV where E~ We ¼ WeT E^ e Ve , A e e e e e ~ condition as system (2.5). However, the matrix E We is dense, in general. Similarly, the matrix W01¼ We1 je ¼ 0 ð2 Rnq Þ is also a column-orthonormal matrix. ^ ^ ^ ^ ¼ colspanfW0 g. Therefore, when the parameter It satisfies the relation Kr ðA E , A BÞ e ¼ 0 in system (2.5), it is a reduced-order system of original system (2.2). And this special reduced-order system is equivalent to what is generated by the direct two-sided Arnoldi MOR algorithm.

Y.-L. Jiang et al. / Journal of the Franklin Institute 349 (2012) 3027–3045

3032

T T T Let We ¼ ½We1 We2    We,kþk T be the partitioning of We corresponding to the block 1 ^ e ¼ diagfWe1 ,We2 , sizes of the matrix Ve in Eq. (2.6). Then, we set W . . . ,We,kþk1 g 2 Rnðkþk1 Þq . A column-orthonormal can be computed by using the GS scheme, given by

~ e ¼ diagfW ~ e1 , W ~ e2 , . . . , W ~ e,kþk1 g: W ~ e has the same block sizes as the matrix V~ e . Then, substituting Suppose that the matrix W ^  V~ e x~ q~ ðtÞ into system (2.3), and premultiplying the state equation the approximation xðtÞ ~ by W e , we derive a reduced-order system d x~ q~ ðtÞ ~ ~ x~ q~ ðtÞ þ B~ ~ uðtÞ, ¼A E~ W~ e q~ W e q~ W e q~ dt

y~ q~ ðtÞ ¼ C~ W~ e q~ x~ q~ ðtÞ,

ð2:9Þ

where ~ ~ T E^ e V~ e ¼ diagfW ~ T V~ e,1 , . . . , W ~ T V~ e,k1 , W ~ T E~ W~ e q~ ¼ W e e,1 e,k1 e,k1 þ1 Ek1 þ1 V e,k1 þ1 , ~ T Ek V~ e,k ,eW ~ T V~ e,kþ1 þ W ~ T N1 V~ e,kþ1 , . . . , . . . ,W e,k e,kþ1 e,kþ1 T T ~ q~ ~ ~ ~ ~ V e,kþk1 þ W Nk1 V e,kþk1 g 2 Rq , eW e,kþk1

~ q~ q ~ ~ ¼W ^ ~ ~ TA A , e Ve 2 R W e q~

e,kþk1

~ ~ T B^ 2 Rqp B~ W~ e q~ ¼ W , e

C~ W~ e q~ ¼ C^ V~ e 2 Rmq~ :

Its initial condition is the same as that of system (2.7). And the matrix E~ W~ e q~ also preserves the block diagonal structure of the matrix E^ e in Eq. (2.3). However, its first k1 block diagonal matrices are not identity matrices. The transfer functions can be written as ~ ~ Þ1 B~ ~ . ~ ~ ðsÞ ¼ C~ ~ ðsE~ ~ A H W e q~ W e q~ W e q~ W e q~ W e q~ Replacing ni1 , ni2, and nj by q~ i , q~ kþi and q~ j , respectively in the matrix P, we get a new 1 T matrix P~ such that P~ ¼ P~ , where i ¼ 1,2, . . . ,k1 and j ¼ k1 þ 1,k1 þ 2, . . . ,k. Substitut1 ing the transformation x~ P~ ðtÞ ¼ P~ x~ q~ ðtÞ into system (2.7), an equivalent system can be obtained, given by d x~ P~ ðtÞ 1 1 ~ ~ ~ ~ ðtÞ þ P~ 1 B~ eq~ uðtÞ, ¼ P~ A P~ E~ eq~ P~ eq~ P x P dt 1 ^ with the initial condition x~ P~ ð0Þ ¼ P~ VeT xð0Þ, where

y~ P~ ðtÞ ¼ C~ eq~ P~ x~ P~ ðtÞ,

ð2:10Þ

1 T T P~ E~ eq~ P~ ¼ diagfIq~ 1 ,eIq~ kþ1 þ V~ e,kþ1 N1 V~ e,kþ1 , . . . ,Iq~ k1 ,eIq~ kþk1 þ V~ e,kþk1 Nk1 V~ e,kþk1 , T T ~ q~ : V~ e,k1 þ1 Ek1 þ1 V~ e,k1 þ1 , . . . , V~ e,k Ek V~ e,k g 2 Rq

~T ~ Tq2 ~ Tq,kþk We can partition x~ q~ ðtÞ as x~ q~ ðtÞ ¼ ½x~ Tq1 ðtÞT , where x~ qi ~ ðtÞ ¼ V ei xi1 ðtÞ 2 ~ ðtÞ x ~ ðtÞ    x ~ 1 q~ i ~T ~T x~ qi Rq~ i ði ¼ 1,2, . . . ,k1 Þ; x~ qi ~ ðtÞ ¼ V ei xi ðtÞ 2 R ði ¼ k1 þ 1,k1 þ 2, . . . ,kÞ; ~ ðtÞ ¼ V ei xi2 ðtÞ 2 q~ i R ði ¼ k þ 1,k þ 2, . . . ,k þ k1 Þ. Then, x~ P~ ðtÞ can be written as follows: T ~ Tq,kþ1 ~ Tq,kþk ~ Tq,k ðtÞ    x~ Tqk ðtÞ x~ Tq,k x~ P~ ðtÞ ¼ ½x~ Tq1 ~ ðtÞ x ~ ~ 1 ðtÞ x ~ ~ 1 þ1 ðtÞ    x ~ ðtÞ : 1

It should be pointed out that the sort order of x~ P~ ðtÞ is the same as that of x(t) in Eq. (2.1). 1 Similarly, by using the transformation x~ P~ ðtÞ ¼ P~ x~ q~ ðtÞ, reduced-order system (2.9) can be changed into d x~ P~ ðtÞ 1 1 ~ ~ 1 x~ ~ ðtÞ þ P~ 1 B~ ~ uðtÞ, ¼ P~ A P~ E~ W~ e q~ P~ ~ e q~ P W P W e q~ dt

y~ P~ ðtÞ ¼ C~ W~ e q~ P~ x~ P~ ðtÞ,

ð2:11Þ

Y.-L. Jiang et al. / Journal of the Franklin Institute 349 (2012) 3027–3045

3033

where 1 ~ T V~ e1 ,eW ~ T V~ e,kþ1 P~ E~ W~ e q~ P~ ¼ diagfW e1 e,kþ1 ~ T N1 V~ e,kþ1 , . . . , W ~ T V~ ek1 ,eW ~ T þW

~

e,kþ1 ek1 e,kþk1 V e,kþk1 T T ~ ~ ~ ~ ~ ~ T þW e,kþk1 Nk1 V e,kþk1 , W e,k1 þ1 Ek1 þ1 V e,k1 þ1 , . . . , W e,k Ek V e,k g

~ q~ 2 Rq :

It has the same initial condition as system (2.10). ^ 1 E^ e , A ^ 1 BÞ ^ ¼ colspanfVe gDcolspanfV~ e g, both According to [17], since Kr ðA reduced-order systems (2.5) and (2.7) match the first r moments of original system (2.3). ^ T E^ T , A ^ T C^ T Þ ¼ colspanfWe gDcolspanfW ~ e g holds. Therefore, The relation Kr^ ðA e systems (2.8) and (2.9) that are constructed by two-sided Arnoldi MOR algorithm preserves the first r þ r^ moments. In addition, because systems (2.10) and (2.11) are equivalent to systems (2.7) and (2.9), respectively, thus they match the first r and r þ r^ moments of system (2.3). 3. Model-order reduction of the subsystems In this section, we directly apply our approach to the reduction of the subsystems. To preserve the special block structure of the subsystems, a structure-preserving MOR algorithm is also used. 3.1. e-embedding of the subsystems From the assumption mentioned in Section 1, coupled system (1.1)–(1.2) composed of k1 DAE subsystems and k2 ODE subsystems. Furthermore, the first k1 subsystems are the DAE systems which can be described as " # " # " # Mi 0 Bi1 Ini1 0 dxi ðtÞ ¼ ð3:1Þ x ðtÞ þ ui ðtÞ, yi ðtÞ ¼ ½Ci1 Ci2 xi ðtÞ, 0 Ini2 i 0 Ni Bi2 dt where Bi1 2 Rni1 pi , Bi2 2 Rni2 pi , Ci1 2 Rmi ni1 , Ci2 2 Rmi ni2 , and i ¼ 1,2, . . . ,k1 . The last k2 subsystems are in the following form: dxj ðtÞ ¼ Aj xj ðtÞ þ Bj uj ðtÞ, yj ðtÞ ¼ Cj xj ðtÞ, j ¼ k1 þ 1,k1 þ 2, . . . ,k: ð3:2Þ dt The coupled relations between these subsystems are still of form (1.2). By using the e-embedding, DAE subsystem (3.1) is changed into an ODE system Ej

dxi ðtÞ ¼ Ai xi ðtÞ þ Bi ui ðtÞ, yi ðtÞ ¼ Ci xi ðtÞ, ð3:3Þ dt where Eei ¼ diagfIni1 ,eIni2 þ Ni g is nonsingular and the system is related to the parameter e. Here e40 and e51. Its initial condition is xi ð0Þ ¼ xi0 which can be rewritten as xi0 ¼ ½xTi1,0 xTi2,0 T , where xi1,0 2 Rni1 and xi2,0 2 Rni2 . The transfer functions of systems (3.1) and (3.3) can be described in the form " #" # Bi1 0 ðsI ni1 Mi Þ1 Hi ðsÞ ¼ ½Ci1 Ci2  Bi2 0 ðsN i Ini2 Þ1 Eei

¼ Ci1 ðsI ni1 Mi Þ1 Bi1 þ Ci2 ðsN i Ini2 Þ1 Bi2 ,

3034

Y.-L. Jiang et al. / Journal of the Franklin Institute 349 (2012) 3027–3045

"

Hei ðsÞ ¼ ½Ci1

ðsI ni1 Mi Þ1 Ci2  0

ðsðeIni2

0 þ Ni ÞIni2 Þ1

#"

Bi1 Bi2

#

¼ Ci1 ðsI ni1 Mi Þ1 Bi1 þ Ci2 ðsðeIni2 þ Ni ÞIni2 Þ1 Bi2 , respectively. Similar to derivation of Eq. (2.4), JHi ðsÞHei ðsÞJ2 can be bounded as follows: JHi ðsÞHei ðsÞJ2 r

ejsjJCi2 J2 JBi2 J2 JðsN i Ini2 Þ1 J22 : 1ejsjJðsN i Ini2 Þ1 J2

Let SDC be a compact set such that for any s 2 S, the matrix sN i Ini2 is nonsingular. Thus, when e-0, we have Hei ðsÞ-Hi ðsÞ. In this case, system (3.3) is a good approximation of subsystem (3.1). 3.2. Model-order reduction of the subsystems In the following, we use one-sided and two-sided Arnoldi algorithms to reduce the each subsystems, respectively. 1 1. One-sided Arnoldi algorithm for the subsystems. For Kri ðA1 i Eei ,Ai Bi Þ, we can get a ni qi column-orthonormal matrix Vei 2 R by the block Arnoldi algorithm. Here qi ¼ ri pi and qi 5ni . Then, taking the approximation x~ i ðtÞ  VeiT xi ðtÞ derives a reduced-order system d x~ i ðtÞ ¼ A~ ei x~ i ðtÞ þ B~ ei ui ðtÞ, y~ i ðtÞ ¼ C~ ei x~ i ðtÞ, E~ ei ð3:4Þ dt ~ ei ¼ V T Ai Vei , B~ ei ¼ V T Bi with the initial condition x~ i ð0Þ ¼ VeiT xi0 , where E~ ei ¼ VeiT Eei Vei , A ei ei 1 ~ ~ ~ ~ ei ðsÞ ¼ C ei ðsE~ ei A ei Þ B~ ei . and C ei ¼ Ci Vei . Its transfer function is H Let V0i ¼ Vei je ¼ 0 2 Rni qi . Thus, V0i ¼ Vei is also a column-orthonormal matrix and 1 satisfies the relation Kri ðA1 i Ei ,Ai Bi Þ ¼ colspanfV0i g. In this case, when the parameter e is chosen as zero in system (3.4), it is regarded as a reduced-order system of original system (3.3). Furthermore, it is also equivalent to the reduced-order system resulting from the direct one-sided Arnoldi algorithm to Eq. (3.3). 1 A column-orthonormal matrix Vj 2 Rnj qj can be generated based on Krj ðA1 j Ej ,Aj Bj Þ, where qj ¼ rj pj ðqj 5nj Þ and j ¼ k1 þ 1,k1 þ 2, . . . ,k. Therefore, a reduced-order system can be obtained from the matrix Vj, given by d x~ j ðtÞ ~ j x~ j ðtÞ þ B~ j uj ðtÞ, y~ j ðtÞ ¼ C~ j x~ j ðtÞ, ¼A E~ j ð3:5Þ dt ~ j ¼ V T Aj Vj , B~ j ¼ V T Bj , C~ j ¼ Cj Vj , and the initial condition is where E~ j ¼ VjT Ej Vj , A j j T x~ j ð0Þ ¼ Vj xj0 . Subsystems (3.4) and (3.5) are coupled through the relations ( ul ðtÞ ¼ Kl1 y~ 1 ðtÞ þ    þ Klk y~ k ðtÞ þ Gl uðtÞ, ð3:6Þ ~ ¼ R1 y~ 1 ðtÞ þ    þ Rk y~ k ðtÞ, l ¼ 1,2, . . . ,k: yðtÞ The transfer functions of systems (3.2) and (3.5) have the form ~ j Þ1 B~ j : ~ j ðsÞ ¼ C~ j ðsE~ j A Hj ðsÞ ¼ Cj ðsE j Aj Þ1 Bj , H ~ ei do not preserve the block structure of the matrices The coefficient matrices E~ ei and A Eei and Ai in Eq. (3.3). Instead of using the matrix Vei directly for the projection, the new MOR algorithm employs a modified version of this matrix which trivially leads to the T T T structure preservation. We partition the column-orthonormal Vei as Vei ¼ ½Vei,1 Vei,2  ,

Y.-L. Jiang et al. / Journal of the Franklin Institute 349 (2012) 3027–3045

3035

where its block sizes correspond to those of the matrices Eei and Ai. We use the GS scheme on the columns of V^ ei ¼ diagfVei,1 ,Vei,2 g to generate a new column-orthonormal matrix V~ ei ¼ diagfV~ ei,1 , V~ ei,2 g 2 Rni q~ i ðq~ i 5ni Þ, where V~ ei,1 2 Rni q~ i1 and V~ ei,2 2 Rni q~ i2 , in which 1 q~ i ¼ q~ i1 þ q~ i2 , q~ il rqi ðl ¼ 1,2Þ and q~ i Zqi . It is easy to see that Kri ðA1 i Eei ,Ai Bi ÞD colspanfV~ ei g. Based on the matrix V~ ei , a new reduced-order system can be computed as follows: d x~ q~ i ðtÞ ~ eq~ x~ q~ ðtÞ þ B~ eq~ ui ðtÞ, ¼A E~ eq~ i i i i dt

y~ q~ i ðtÞ ¼ C~ eq~ i x~ q~ i ðtÞ,

ð3:7Þ

T

where the initial condition is x~ q~ i ð0Þ ¼ V~ ei xi0 , and " T E~ eq~ i ¼ V~ ei Eei V~ ei

¼ "

~ eq~ ¼ V~ T Ai V~ ei A ei i

¼

Iq~ i1 0

2

#

0

, T eIq~ i2 þ V~ ei,2 Ni V~ ei,2

T V~ ei,1 Mi V~ ei,1

0

0

Iq~ i2

#

T B~ eq~ i ¼ V~ ei Bi

¼4

T V~ ei,1 Bi1 T V~ ei,2 Bi2

"

,

C~ eq~ i

ðCi1 V~ ei,1 ÞT ~ ¼ Ci V ei ¼ ðCi2 V~ ei,2 ÞT

3 5,

#T :

~ eq~ have the same diagonal block structure as the matrices Eei and The matrices E~ eq~ i and A i Ai in Eq. (3.3). From Eq. (3.6), we have 8 < ul ðtÞ ¼ Kl1 y~ q~ 1 ðtÞ þ    þ Klk1 y~ q~ k ðtÞ þ Kl,k1 þ1 y~ k1 þ1 ðtÞ þ    þ Klk y~ k ðtÞ þ Gl uðtÞ, 1

~ ¼ R1 y~ q~ 1 ðtÞ þ    þ Rk1 y~ q~ k ðtÞ þ Rk1 þ1 y~ k1 þ1 ðtÞ þ    þ Rk y~ k ðtÞ, : yðtÞ 1 Pk1

l ¼ 1,2, . . . ,k: ð3:8Þ

Pk

~ Let q~ ¼ i ¼ 1 q~ i þ j ¼ k1 þ1 qj . Coupled system (3.5), (3.7), and (3.8) is of order q. 2. Two-sided Arnoldi algorithm for the subsystems. We construct a column-orthonormal matrix Wei 2 Rni qi , the columns of which are a basis of the Krylov subspace T T T Kr^ i ðAT i Ei ,Ai Ci Þ, where qi ¼ r^ i mi ¼ ri pi . Then, a reduced-order system is obtained, given by d x~ Wei ðtÞ ~ W x~ W ðtÞ þ B~ W ui ðtÞ, ¼A E~ Wei ei ei ei dt

y~ Wei ðtÞ ¼ C~ Wei x~ Wei ðtÞ,

ð3:9Þ

with the same initial condition as Eq. (3.4), where E~ Wei ¼ WeiT Eei Vei , A~ Wei ¼ WeiT Ai Vei , B~ Wei ¼ WeiT Bi and C~ Wei ¼ Ci Vei . Similarly, if we choose the parameter e ¼ 0, system (3.9) is a reduced-order system of original system (3.1). Furthermore, it is equivalent to the system constructed by the direct two-sided Arnoldi MOR algorithm. Next, we discuss the structure-preserving in details. The matrix Wei can be partitioned as T T T Wei ¼ ½Wei,1 Wei,2  with the same block sizes of the matrix Vei . Suppose that the two ~ ei,1 2 Rni q~ i1 and W ~ ei,2 2 Rni q~ i2 can be obtained by column-orthonormal matrices W running the GS scheme on the columns of the matrices Wei,1 and Wei,2 , respectively. Let ~ ei ¼ diagfW ~ ei,1 , W ~ ei,2 g. We have Kr^ i ðAT E T ,AT C T Þ ¼ colspanfWei gDcolspan W i i i i ~ ei g. fW ~ ei and V~ ei , we get a reduced-order system as follows: According to the matrices W d x~ q~ i ðtÞ ~ eIi x~ q~ ðtÞ þ B~ eIi ui ðtÞ, ¼A E~ eIi i dt

y~ q~ i ðtÞ ¼ C~ eIi x~ q~ i ðtÞ,

ð3:10Þ

3036

Y.-L. Jiang et al. / Journal of the Franklin Institute 349 (2012) 3027–3045

where the initial condition is the same as Eq. (3.7), and 2 T 3 ~ ~ W 0 ei,1 V ei,1 T 5, ~ Eei V~ ei ¼ 4 E~ eIi ¼ W ei ~ T Ni V~ ei,2 ~ T V~ ei,2 þ W 0 eW ei,2 ei,2 2 T 3 ~ W Bi1 ~ T Bi ¼ 4 ei,1 5, B~ eIi ¼ W ei T ~ W ei,2 Bi2 2 T 3 " #T ~ ~ W 0 ðCi1 V~ ei,1 ÞT ei,1 Mi V ei,1 T 5, C~ eIi ¼ Ci V~ ei ¼ ~ Ai V~ ei ¼ 4 A~ eIi ¼ W : ei ~ T V~ ei,2 ðCi2 V~ ei,2 ÞT 0 W ei,2 ~ Obviously, the matrices E~ eIi and Coupled system (3.5), (3.8), and (3.10) is also of order q. ~ A eIi preserve the diagonal block structure of the matrices Eei and Ai in Eq. (3.3). However, ~ T V~ ei,2 are dense in general. ~ T V~ ei,1 and W the matrices W ei,1 ei,2 1 It can be seen from [17] that since Kri ðA1 i Eei ,Ai Bi Þ ¼ colspanfVei g, reduce-order systems (3.4) and (3.7) have the same first ri moments of original system (3.3). Moreover, system (3.5) preserves the first rj moments of subsystem (3.2), because of 1 Krj ðA1 j Ej ,Aj Bj Þ ¼ colspanfVj g. In addition, note that T T T ~ Kr^ i ðAT i Ei ,Ai Ci Þ ¼ colspanfWei gDcolspanfW ei g:

Therefore, both systems (3.9) and (3.10) which are resulting from the two-sided Arnoldi algorithm match the first ri þ r^ i moments of Eq. (3.3). 4. Error estimate, stability, and passivity In this section, we discuss the error between coupled system (1.1)–(1.2) and its reducedorder systems resulting from the MOR methods given in Sections 2 and 3. Furthermore, the stability and passivity for our MOR approach are also considered. 4.1. Error estimate In the following, we first estimate the error between closed-loop system (2.1) and ^ 1 E^ e , A ^ 1 BÞ, ^ we reduced-order system (2.5). Applying the block Arnoldi algorithm to Kr ðA get the matrix Ve ¼ ½V0e V1e    Vr1,e  and an upper block Hessenberg matrix H ¼ ðHil Þrr 2 Rqq , where Vei 2 Rnp ði ¼ 0,1, . . . ,r1Þ and Hil 2 Rpp ðl ¼ 0,1, . . . ,r1Þ. According to [20], the matrices Ve and H satisfy the following relations: 2 3 ^ 1 E^ e Ve , H ¼ VeT A

^ 1 E^ e Ve ¼ Ve H þ Vre Hr,r1 40    0 Ip 5: A |fflfflfflffl{zfflfflfflffl}

ð4:1Þ

r1

~ 1 V T  ^ eA Let Oðe,sÞ ¼ jsjJCJ2 JBJ2 JðsEAÞ1 J2 =ð1ejsjJðsEAÞ1 J2 Þ and LðeÞ ¼ ðAV e e ^ re Hr,r1 ½0    0 Ip . From [16], the error JHðsÞH ~ e ðsÞJ2 can be estimated as In ÞAV ~ e Þ1 J2 JLðeÞJ2 Þ: ~ e ðsÞJ2 rOðe,sÞðeJðsEAÞ1 J2 þ JW J2 JðsE~ e A JHðsÞH

Y.-L. Jiang et al. / Journal of the Franklin Institute 349 (2012) 3027–3045

3037

Since the matrix V~ e satisfies colspanfVe gDcolspanfV~ e g, there exists a matrix C 2 T ~ Rqq such that Ve ¼ V~ e C. Clearly, we have C ¼ V~ e Ve and CT C ¼ Iq . Substituting Ve ¼ V~ e C into Eq. (4.1) yields equivalent relations T

^ H ¼ CT V~ e A

1

E^ e V~ e C,

^ 1 E^ e V~ e C ¼ V~ e CH þ Vre Hr,r1 ½ 0    0 Ip : A

ð4:2Þ

^ 1 E^ e , A ^ 1 BÞDcolspanf ^ In addition, the relation Kr ðA V~ e g is true. There exists a matrix 1 ~ ^ B^ ¼ V~ e F0 . Then, the matrix B~ eq~ in Eq. (2.7) satisfies the such that A F0 2 Rqp ~ 1 B~ eq~ ¼ A ^ 1 B. ^ In the following, we estimate the error between the transfer expression V~ e A eq~

functions of systems (2.1) and (2.7). ~ eq~ ðsÞ be the transfer functions of closed-loop system (2.1) and Theorem 4.1. Let H(s) and H reduced-order system (2.7), respectively. If the parameter e satisfies the relation ~ eq~ ðsÞJ2 can be bounded as ejsjJðsEAÞ1 J2 o1, then the error JHðsÞH ~ eq~ Þ1 J2 ðJL1 ðeÞJ2 þ JL2 ðeÞJ2 ÞÞ, ~ eq~ ðsÞJ2 rOðe,sÞðeJðsEAÞ1 J2 þ JðsE~ eq~ A JHðsÞH ~ 1 E~ T E^ e V~ e ÞðIq~ CCT Þ and ^ V~ e A where L1 ðeÞ ¼ ðA eq~ eq~ ~ 1 V~ T In ÞAV ^ V~ e A ^ re Hr,r1 ½0    0 Ip CT : L2 ðeÞ ¼ ðA e eq~ ~ eq~ ðsÞJ2 of the form Proof. We compute the error JHðsÞH ~ eq~ ðsÞJ2 ¼ JH ^ ðsÞH ~ eq~ ðsÞJ2 rJH ^ ðsÞH ^ e ðsÞJ2 þ JH ^ e ðsÞH ~ eq~ ðsÞJ2 : JHðsÞH

ð4:3Þ

^ e ðsÞH ~ eq~ ðsÞJ2 . By using the ^ ðsÞH ^ e ðsÞJ2 and JH This implies that we need to estimate JH 1 1 ~ B~ eq~ ¼ A ^ B, ^ we have equality V~ e A eq~

~ 1 E~ eq~ E^ e V~ e ÞðsE~ eq~ A ^ 1 ðA ^ V~ e A ~ eq~ Þ1 B~ eq~ J2 : ^ E^ e AÞ ~ eq~ ðsÞJ2 ¼ JsCðs ^ e ðsÞH JH eq~

ð4:4Þ

It is not difficult to see that from Eq. (4.2), we get ^ V~ e A ~ 1 E~ eq~ E^ e V~ e ¼ L1 ðeÞ þ ðA ~ 1 E~ eq~ E^ e V~ e ÞCCT ^ V~ e A A eq~ eq~ ^ 1 E^ e V~ e CCT Þ, ^ V~ e ðA ~ 1 E~ eq~ CCT CCT V~ T A ¼ L1 ðeÞDðeÞ þ A e eq~ ^ re Hr,r1 ½0    0 Ip CT . Using Eq. (4.2) again, we can compute the matrix where DðeÞ ¼ AV ^ 1 E^ e V~ e CCT as ~ 1 E~ eq~ CCT CCT V~ T A A eq~

e

~ 1 E~ eq~ CCT CCT V~ T A ^ 1 E^ e V~ e CCT ¼ A ^ 1 E^ e V~ e CCT Þ ~ 1 ðE~ eq~ CCT A ~ eq~ CCT V~ T A A e e eq~ eq~ ^ re Hr,r1 ½0    0 Ip CT : ~ 1 V~ T AV ¼A e eq~ ~ eq~ ðsÞJ2 can be written as ^ e ðsÞH From expression (4.4), JH ^ 1 ðL1 ðeÞ þ L2 ðeÞÞðsE~ eq~ A ~ eq~ ÞB~ eq~ J2 : ^ E^ e AÞ ~ eq~ ðsÞJ2 ¼ JsCðs ^ e ðsÞH JH ^ 1 J2 ¼ JðsEAÞ1 J2 , according to ^ 2 ¼ JCJ2 , JBJ ^ 2 ¼ JBJ2 , and JðsE^ AÞ Notice that JCJ Eqs. (2.4) and (4.3), we can get the conclusion of Theorem 4.1. This completes the proof. &

3038

Y.-L. Jiang et al. / Journal of the Franklin Institute 349 (2012) 3027–3045

Similarly, for system (2.9) resulting from the structure-preserving MOR algorithm, if ejsjJðsEAÞ1 J2 o1, the error between the transfer functions of Eqs. (2.1) and (2.9) is ~ ~ ðsÞJ2 rOðe,sÞðeJðsEAÞ1 J2 JHðsÞH W e q~

~ ~ Þ1 J2 ðJL^ 1 ðeÞJ2 þ JL^ 2 ðeÞJ2 ÞÞ, þ JðsE~ W~ e q~ A W e q~ ~ 1 ^ V~ e A ~ ~ E^ e V~ e ÞðIq~ CCT Þ and ^ 1 ðeÞ ¼ ðA where L ~ e q~ E W ~ W eq ~ 1 ^ V~ e A ^ re Hr,r1 ½0    0 Ip CT : ~ T In ÞAV L^ 2 ðeÞ ¼ ðA ~ e q~ W e W In the following, the error estimation of the reduced-order systems which are constructed by the subsystem MOR algorithms mentioned in Section 3, can be derived. However, it is not very suitable to measure the error directly, rather it is useful as a guide to understand the convergence properties. A column-orthonormal matrix Vei and a qi  qi block upper Hessenberg matrix Hei can be 1 generated from the block Krylov subspace Kri ðA1 i Eei ,Ai Bi Þ ði ¼ 1,2, . . . ,k1 Þ. They have the form Vei ¼ ½V0,ei V1,ei    Vri 1,ei  and Hei ¼ ðHjl,ei Þri ri . In addition, based on the Krylov 1 subspace Krj ðA1 j Ej ,Aj Bj Þðj ¼ k1 þ 1,k1 þ 2, . . . ,kÞ, we obtain the matrix Vj ¼ ½V0j V1j    Vrj 1,j  and an upper block Hessenberg matrix Hj ¼ ðHil,j Þrj rj . Let Oi ðe,sÞ ¼ ~ ei ðsÞ and jsjJCi J2 JBi J2 JðsE i Ai Þ1 J2 =ð1ejsjJðsE i Ai Þ1 J2 Þ and Hl ðsÞ ðl ¼ 1,2, . . . ,kÞ, H ~ H j ðsÞ denote the transfer functions of systems (1.1), (3.4), and (3.5), respectively. If the relation ~ ei ðsÞJ2 and JHj ðsÞH ~ j ðsÞJ2 can be estimated as ejsjJðsE i Ai Þ1 J2 o1 holds, JHi ðsÞH ( ~ ei Þ1 J2 Þ, ~ ei ðsÞJ2 rOi ðe,sÞðeJðsE i Ai Þ1 J2 þ JL~ i ðeÞJ2 JðsE~ ei A JHi ðsÞH ð4:5Þ ~ j Þ1 J2 , ~ j J2 JðsE j Aj Þ1 J2 JðsE~ j A ~ j ðsÞJ2 rjsjJCj J2 JBj J2 JL JHj ðsÞH ~ i ðeÞ ¼ ðAi Vei A~ ei V T In ÞAi Vr ,ei Hr ,r 1,ei ½0    0 Ip  and where L i i i i i ei ~ 1 V T In ÞAj Vj Hr ,r 1,j ½0    0 Ip : L~ j ¼ ðAj Vj A j j j j j j Following, we establish the fact on the error between the transfer functions of coupled systems (1.1)–(1.2) and (3.4)–(3.6). The transfer function of (1.1)–(1.2) can be written as HðsÞ ¼ RðIm0 H0 ðsÞKÞ1 H0 ðsÞG, where H0 ðsÞ ¼ diagfH1 ðsÞ, . . . ,Hk ðsÞg. Similarly, for coupled system (3.4)–(3.6), its transfer function is in the form ~ eI ðsÞ ¼ RðIm0 H ~ eI0 ðsÞG, ~ eI0 ðsÞKÞ1 H H ~ eI0 ðsÞ ¼ diagfH ~ e1 ðsÞ, . . . , H ~ ek1 ðsÞ, H ~ k1 þ1 ðsÞ, . . . , H ~ k ðsÞg. where H Theorem 4.2. Let a ¼ JðIm0 H0 ðsÞKÞ1 J2 , P1 ðs,e,iÞ and P2 ðs,jÞ be the right parts of Eq. (4.5). Assume that ejsjJðsE i Ai Þ1 J2 o1, aJKJ2 P1 ðs,e,iÞo1 and aJKJ2 P2 ðs,jÞo1 for ~ eI ðsÞJ2 can be bounded as i ¼ 1, . . . ,k1 , j ¼ k1 þ 1, . . . ,k. Then, the error JHðsÞH aJRJ2 JGJ2 ð1 þ aJH0 ðsÞJ2 JKJ2 Þ imax fP1 ðs,e,iÞ,P2 ðs,jÞg ¼ 1,...,k , ~ eI ðsÞJ2 r JHðsÞH

1 j ¼ k1 þ1,...,k

1aJKJ2 i max fP1 ðs,e,iÞ,P2 ðs,jÞg: ¼ 1,...,k , 1 j ¼ k1 þ1,...,k

~ eI ðsÞJ2 can be computed as Proof. The error JHðsÞH ~ eI ðsÞJ2 ¼ JRððIm0 H0 ðsÞKÞ1 ðIm0 H ~ e0 ðsÞKÞ1 ÞH0 ðsÞG JHðsÞH ~ e0 ðsÞKÞ1 ðH0 ðsÞH ~ e0 ðsÞÞGJ2 þRðIm0 H

Y.-L. Jiang et al. / Journal of the Franklin Institute 349 (2012) 3027–3045

3039

~ e0 ðsÞKÞ1 J2 JH0 ðsÞH ~ e0 ðsÞJ2 rJRJ2 JGJ2 JðIm0 H ð1 þ aJKJ2 JH0 ðsÞJ2 Þ: It is easy to see that ~ e0 ðsÞKÞ1 ¼ ðIm0 ðIm0 H ~ e0 ðsÞKÞ1 ðH0 ðsÞH ~ e0 ðsÞÞKÞðIm0 H0 ðsÞKÞ1 : ðIm0 H ~ e0 ðsÞKÞ1 J2 is rewritten in the form Thus, JðIm0 H a ~ e0 ðsÞKÞ1 J2 r JðIm0 H : ~ e0 ðsÞJ2 JKJ2 1aJH0 ðsÞH ~ e0 ðsÞJ2 can be estimated as follows: In addition, JH0 ðsÞH ~ eI ðsÞJ2 ¼ JdiagfH1 ðsÞH ~ e1 ðsÞ, . . . ,Hk1 ðsÞH ~ ek1 ðsÞ, JHðsÞH ~ ~ Hk1 þ1 ðsÞH k1 þ1 ðsÞ, . . . ,Hk ðsÞH k ðsÞgJ2 : Note that JdiagfA1 ,A2 gJ2 ¼ maxfJA1 J2 ,JA2 J2 g. Therefore, we have ~ eI ðsÞJ2 ¼ max fJHi ðsÞH ~ ei ðsÞJ2 ,JHj ðsÞH ~ j ðsÞJ2 g: JHðsÞH i ¼ 1,...,k , 1 j ¼ k1 þ1,...,k

From Eq. (4.5), we can obtain the result of Theorem 4.2. This completes the proof.

&

Next, consider the error estimation of the reduced-order systems which are generated by the structure-preserving MOR algorithms in Section 3. Since the expression colspanfVei gDcolspanfV~ ei g holds, there exists a matrix Ci 2 Rq~ i qi such that ~ eq~ ðsÞ and H ~ eIi ðsÞ be the transfer functions of systems (3.7) and (3.10), Vei ¼ V~ ei Ci . Let H i ~ eq~ ðsÞJ2 and respectively. Suppose that ejsjJðsE i Ai Þ1 J2 o1. Then, JHi ðsÞH i ~ JHi ðsÞH eIi ðsÞJ2 have the form ( ^ i2 ðeÞJ2 ÞÞ, ~ eq~ ðsÞJ2 rOi ðe,sÞðeJðsE i Ai Þ1 J2 þ JðsE~ eq~ A~ eq~ Þ1 J2 ðJL^ i1 ðeÞJ2 þ JL JHi ðsÞH i

i

i

~ i2 ðeÞJ2 ÞÞ, ~ eIi ðsÞJ2 rOi ðe,sÞðeJðsE i Ai Þ1 J2 þ JðsE~ eIi A~ eIi Þ1 J2 ðJL~ i1 ðeÞJ2 þ JL JHi ðsÞH

ð4:6Þ ~ 1 E~ eq~ Eei V~ ei ÞðIq~ Ci CT Þ, ^ i1 ðeÞ ¼ ðAi V~ ei A L eq~ i T i i i 1

where Ci CTi Þ, L^ i2 ðeÞ ¼ ðAi V~ ei A~ eq~ i V~ ei Ini ÞAi Vri ,ei Hri ,ri 1,ei ½0

~ 1 E~ eIi Eei V~ ei ÞðIq~  ~ i1 ðeÞ ¼ ðAi V~ ei A L eIi i T    0 Ipi Ci , and

~ 1 W ~ i2 ðeÞ ¼ ðAi V~ ei A ~ T Ini ÞAi Vri ,ei Hri ,ri 1,ei ½0    0 Ipi CT : L i eIi ei Applying the Laplace transformation to coupled system (3.5), (3.7), and (3.8) yields its ~ eII0 ðsÞG, where ~ eII0 ðsÞKÞ1 H ~ eII ðsÞ ¼ RðIm0 H transfer function H ~ eq~ ðsÞ, . . . , H ~ eq~ ðsÞ, H ~ k1 þ1 ðsÞ, . . . , H ~ k ðsÞg: ~ eII0 ðsÞ ¼ diagfH H 1

k1

~ eIII ðsÞ be the transfer function of coupled system (3.5), (3.8), and (3.10). Obviously, Let H ~ ~ eIII0 ðsÞ ¼ RðIm0 H ~ eIII0 ðsÞKÞ1 H ~ eIII ðsÞ ¼ ~ eIII0 ðsÞG, where H H eIII ðsÞ has the form H ~ ~ ~ ~ ~ diagfH eI1 ðsÞ, . . . , H eIk1 ðsÞ, H k1 þ1 ðsÞ, . . . , H k ðsÞg. Some results on JHðsÞH eII ðsÞJ2 and ~ eIII ðsÞJ2 can be stated in the following theorem. JHðsÞH ~ eII ðsÞJ2 and JHðsÞH ~ eIII ðsÞJ2 . For the error JHðsÞH ^ 1 ðs,e,iÞ and P ~ 1 ðs,e,iÞ be the right parts of (4.6). If the relations Theorem 4.3. Let P ejsjJðsE i Ai Þ1 J2 o1,

^ 1 ðs,e,iÞ, P ~ 1 ðs,e,iÞ,P2 ðs,jÞgo1 aJKJ2 i max fP ¼ 1,...,k , 1 j ¼ k1 þ1,...,k

3040

Y.-L. Jiang et al. / Journal of the Franklin Institute 349 (2012) 3027–3045

~ eII ðsÞJ2 and JHðsÞH ~ eIII ðsÞJ2 have the form hold, the error JHðsÞH ^ 1 ðs,e,iÞ,P2 ðs,jÞg aJRJ2 JGJ2 ð1 þ aJH0 ðsÞJ2 JKJ2 Þ imax fP ¼ 1,...,k , 1 j ¼ k1 þ1,...,k

~ eII ðsÞJ2 r JHðsÞH

^ 1 ðs,e,iÞ,P2 ðs,jÞg, 1aJKJ2 i max fP ¼ 1,...,k , 1 j ¼ k1 þ1,...,k

~ 1 ðs,e,iÞ,P2 ðs,jÞg aJRJ2 JGJ2 ð1 þ aJH0 ðsÞJ2 JKJ2 Þ imax fP ¼ 1,...,k , ~ eIII ðsÞJ2 r JHðsÞH

1 j ¼ k1 þ1,...,k

~ 1 ðs,e,iÞ,P2 ðs,jÞg: fP 1aJKJ2 i max ¼ 1,...,k , 1 j ¼ k1 þ1,...,k

The proof of Theorem 4.3 follows the same way as the proof of Theorem 4.2. Here, we omit the details. 4.2. Stability and passivity In this section, for the systems resulting from the Krylov subspace MOR methods in Sections 2 and 3, we consider their stability and passivity. Given a matrix A 2 Rnn , if all of its eigenvalues have strictly negative real part, then it is said to be a Hurwitz matrix. It should be noted that in general, the ODE systems which are constructed by the e-embedding method may not preserve the stability of the original DAE systems. We ^ and B^ of system (2.3) as follows: partition the coefficient matrices E^ e , A, " # " # ^ ^ ^ ^ ¼ A 11 A 12 , B^ ¼ B 1 , ^E e ¼ diagfE^ 1 ,eIn^ þ Ng, A 2 ^ ^ B^ 2 A 21 A 22 ^ 11 2 Rðnn^ 2 Þðnn^ 2 Þ , and B^ 1 2 Rðnn^ 2 Þp . Then, the state variable where E^ 1 2 Rðnn^ 2 Þðnn^ 2 Þ , A T T T ^ In the ^ ¼ ½x^ 1 ðtÞ x^ 2 ðtÞ is the partitioning of xðtÞ ^ corresponding to the block sizes of A. xðtÞ ^ ^ ^ ^ case that the matrix N ¼ 0, and that A 0 þ B 0 K1 and A 22 þ B 2 K2 are Hurwitz. Here, ^ 21 Þ, and ^ 0 ¼ E^ 1 ðA ^ 11 ðA ^ 12 þ B^ 1 K2 ÞðA ^ 22 þ B^ 2 K2 Þ1 A K1 2 Rpðnn^ 2 Þ , K2 2 Rpn^ 2 , A 1

1 ^ 12 þ B^ 1 K2 ÞðA ^ 22 þ B^ 2 K2 Þ1 B^ 2 Þ. Therefore, according to [21] system (2.3) B^ 0 ¼ E^ 1 ðB^ 1 ðA with the special input variable uðtÞ ¼ K1 x^ 1 ðtÞ þ K2 x^ 2 ðtÞ is asymptotically stable for each e 2 ð0,e0 Þ, where e0 is a positive value. Similarly, for system (3.3) with Eei ¼ diagfIni1 ,eIni2 g and uðtÞ ¼ ½K^ i1 K^ i2 xi ðtÞ, where K^ i1 2 Rpi ni1 and K^ i2 2 Rpi ni2 , if the two matrices Ini2 þ Bi2 K^ i2 and Mi þ Bi1 ðIni1 K^ i2 ðIni2 þ Bi2 K^ i2 Þ1 Bi2 ÞK^ i1 are Hurwitz, system (3.3) is asymptotically stable for e 2 ð0,ei0 Þ. Here, ei0 is a positive value. In the following, the stability of systems (2.8) and (3.9) is considered, respectively. First, we give a fact on the stability of system (2.8). 1 ^ is Hurwitz and there exists a matrix We such that Theorem 4.4. If E^ e A

WeT E^ e Ve 40,

^ e þ V TA ^ T We o0, WeT AV e

then, both systems (2.3) and (2.8) are asymptotically stable.

ð4:7Þ

Y.-L. Jiang et al. / Journal of the Franklin Institute 349 (2012) 3027–3045

3041

^ is Hurwitz, ^ in Eq. (2.3) such that E^ 1 A Proof. Note that the coefficient matrices E^ and A this implies that system (2.3) is asymptotically stable. From Eq. (4.7), we have ^ e þ V TA ^ T We o0. Let le 2 sðE~ W , A ~W þA ~ T ¼ W T AV ~ W Þ and E~ We ¼ WeT E^ e Ve 40 and A e e e We e e xle 2 Cq be a right eigenvector corresponding to the eigenvalue le . Thus, xHle ~ W Þxl ¼ 0 holds. Furthermore, we get xH ðl e E~ T A ~ T Þxl þ xH ðle E~ W  A ~WÞ ðle E~ We A e e e e e We We le le T ~W þA ~ o0, it is easy to see that ðle þ l e Þ xle ¼ 0. Since E~ We 40 and A e W e ~W þA ~ T Þxl o0 and xH E~ W xl 40, which means that Reðle Þo0. TherexHle E~ We xle ¼ xHle ðA e e e e We le fore, system (2.5) is also asymptotically stable. This completes the proof. &

Moreover, we have the similar conclusion on the stability of system (3.9). Assume that Eei1 Ai is Hurwitz and Wei is selected satisfying WeiT Eei Vei 40,

WeiT Ai Vei þ VeiT ATi Wei o0:

ð4:8Þ

Then, system (3.9) preserves the asymptotical stability of system (3.3). In the following, we present some results on the passivity of systems (2.8) and (3.9). From [22], a matrix-valued function H: C-Cpp is positive real if and only if the following conditions hold: (1) H(s) is analytic in Cþ ¼ fs 2 C : ReðsÞ40g, (2) HðsÞ ¼ HðsÞ for all s 2 C, and (3) HðsÞ þ H H ðsÞZ0 for all s 2 Cþ . In addition, when the transfer function of a linear time-invariant system is positive real, the system is passive. 1 ^ is Hurwitz and system (2.3) is passive. If we Theorem 4.5. Suppose that the matrix E^ e A T choose a matrix We satisfying (4.7) and WeT B^ ¼ VeT C^ , then, system (2.8) preserves the asymptotical stability and passivity of system (2.3).

The proof of Theorem 4.5 can be found in [23]. Furthermore, let Eei1 Ai be Hurwitz and T T ^ system (3.3) be passive. If there is a matrix Wei satisfying Eq. (4.8) and We,i B i ¼ VeiT C^ i , then, reduced-order system (3.9) preserves the asymptotical stability and passivity of original system (3.3).

5. Numerical experiments In this section, we provide two numerical examples to illustrate the effectiveness of our MOR approach for coupled DAE systems. Example 5.1. Consider a heated beam whose temperature is steered by a controller which is derived from [22]. The controller can be described as the descriptor system E1

dx1 ðtÞ ¼ A1 x1 ðtÞ þ b1 u1 ðtÞ, dt

y1 ðtÞ ¼ c1 x1 ðtÞ,

where 2 E1 ¼ diagfI200 ,0800 g,

A1 ¼ I1000 ,

b1 ¼ ½bT11 bT12 T ,

and

3

c1 ¼ 41 0    0 15, |fflfflfflffl{zfflfflfflffl} 998

Y.-L. Jiang et al. / Journal of the Franklin Institute 349 (2012) 3027–3045

3042

in which 2

3T

b11 ¼ 40:1    0:1 5 |fflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflffl}

2 and

200

3T

b12 ¼ 40:1    0:1 5 : |fflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflffl} 800

By a spatial discretization of the heated beam, we obtain E2 dxdt2 ðtÞ ¼ A2 x2 ðtÞ þ b2 u2 ðtÞ, y2 ðtÞ ¼ c2 x2 ðtÞ, where E2 ¼ I4 and 2 3 2 3 2 3T 25 25 0 0 5 0 6 25 50 25 7 607 607 0 6 7 6 7 6 7 A2 ¼ 6 7 , b2 ¼ 6 7 , c 2 ¼ 6 7 : 4 0 405 405 25 50 25 5 0 0 25 25 0 1 The interconnection of the controller and the beam is expressed by the relations u1 ðtÞ ¼ uðtÞy2 ðtÞ, u2 ðtÞ ¼ y1 ðtÞ, and yðtÞ ¼ y2 ðtÞ. The initial conditions of the above coupled system are 2 3T x1 ð0Þ ¼ 40 0    0 5 |fflfflfflfflfflffl{zfflfflfflfflfflffl}

and

x2 ð0Þ ¼ ½0 0 0 1T :

1000

Its external input variable is given by uðtÞ ¼ e10t . The positive parameter e is chosen as e ¼ 1012 . Since the order of the first subsystem is 1000, it should be approximated by our MOR approach on the interval [0, 1]. Moreover, its reduced-order is set to 10. For comparison, the direct Arnoldi algorithm and the POD method [10] are used to 10−2

10−4

10−6

10−8

10−10

10−12

10−14

Embedding MOR Direct Arnoldi MOR POD MOR

10−16 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

Fig. 1. Relative error for Example 5.1.

0.8

0.9

1

Y.-L. Jiang et al. / Journal of the Franklin Institute 349 (2012) 3027–3045

3043

approximate the first subsystem. The reduction order for the direct Arnoldi algorithm is R1 10. Since the first two eigenvalues of the correlation matrix M1 ¼ 0 x1 ðtÞxT1 ðtÞ dt are much larger than the others, then the reduction order is chosen as 2. Fig. 1 shows the relative error between the original coupled and its reduced systems. It is easy to see that our approach is better than the direct Arnoldi algorithm and the POD method. According to Table 1, the simulation time of our approach and the direct Arnoldi algorithm are less than that of the POD method. Example 5.2. A coupled system with two subsystems is discussed. The first subsystem has the form E1 dxdt1 ðtÞ ¼ A1 x1 ðtÞ þ b1 u1 ðtÞ, y1 ðtÞ ¼ c1 x1 ðtÞ, where 2 E1 ¼ C1 diagfI200 ,01000 g,

2

3T

b 1 ¼ 41 0    0 5 , |fflfflfflffl{zfflfflfflffl} 1199

3

c1 ¼ 40    0 15, |fflfflfflffl{zfflfflfflffl}

"

and

A11 A1 ¼ A21

# A12 , A22

1199

Table 1 Simulation time comparison of the three methods. Model

Full (s)

Embedding (s)

Arnoldi (s)

POD (s)

Example 5.1 Example 5.2

10.0624 34.0302

7.6355 20.1366

1.6881 1.9729

1.0157 1.2789

10−1

Embedding MOR Direct Arnoldi MOR POD MOR

10−2

10−3

10−4

10−5

10−6

10−7 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

Fig. 2. Relative error for Example 5.2.

0.8

0.9

1

Y.-L. Jiang et al. / Journal of the Franklin Institute 349 (2012) 3027–3045

3044

in which A11 is a zero matrix of order 600  600, A22 ¼ R1 I600 and 2 3 2 3 1 1 1 6 1 1 7 6 7 1 & 6 7 6 7 A12 ¼ 6 7, A21 ¼ 6 7: 4 5 4 & & & 1 5 1

1

1

The second subsystem can be written as follows: " #     0 0 1 1 1 dx2 ðtÞ ¼ x2 ðtÞ þ u2 ðtÞ, 0 C2 R2 dt 1 0 0

y2 ðtÞ ¼ ½0 R2 x2 ðtÞ:

The above two subsystems are coupled through the input–output relations u1 ðtÞ ¼ uðtÞ, u2 ðtÞ ¼ y1 ðtÞ, and yðtÞ ¼ y2 ðtÞ. The initial conditions are 2 3T x1 ð0Þ ¼ 40    0 5 |fflfflfflffl{zfflfflfflffl}

and

x2 ð0Þ ¼ ½0 1T :

1200

The external input is uðtÞ ¼ 20e10t . The positive parameter e is also chosen as e ¼ 1012 . Here, we give the parameter values R1 ¼ 10, R2 ¼ 100, and C1 ¼ C2 ¼ 109 . On the interval [0, 1], the first subsystem is reduced by our approach with a reduced model of order 10. Similarly, we use the direct Arnoldi algorithm to reduce the order to 10. Furthermore, for the POD method [10], we set the reduction order to 2, because the first two eigenvalues of its correlation matrix are much larger. In Fig. 2, the relative error between the original coupled and its reduced systems are shown. Clearly, the relative error resulting from our MOR approach is smaller than those resulting from the direct Arnoldi algorithm and the POD method. Table 1 shows that the POD method has much more simulation time. 6. Conclusions In this paper, we have discussed the reduction of linear coupled systems including DAE subsystems via e-embedding Krylov subspace method. The investigated technique transfers DAE subsystems into ODE systems. Then, both one-sided and two-sided Arnoldi algorithms are applied to reduction of the closed-loop system and the subsystems, respectively. In addition, we conclude some results on the error, stability, and passivity for the resulting systems. Acknowledgments This work was supported by the Natural Science Foundation of China (NSFC) under Grant 11071192 and the International Science and Technology Cooperation Program of China under Grant 2010DFA14700.

References [1] T. Bechtold, Model order reduction of electro-thermal MEMS, Ph.D. Dissertation, Albert-Ludwigs University, 2005.

Y.-L. Jiang et al. / Journal of the Franklin Institute 349 (2012) 3027–3045

3045

[2] C.A. Felippa, K.C. Park, C. Farhat, Partitioned analysis of coupled mechanical systems, Computer Methods in Applied Mechanics and Engineering 190 (2001) 3247–3270. [3] C. Tischendorf, Coupled systems of differential algebraic and partial differential equations in circuit and device simulation: modeling and numerical analysis, Habilitation Dissertation, Humboldt University, 2003. [4] B. Simeon, DAE’s and PDE’s in elastic multibody systems, Numerical Algorithms 19 (1998) 235–246. [5] G. Parmar, S. Mukherjee, R. Prasad, System reduction using eigen spectrum analysis and Pade approximation technique, International Journal of Computational and Applied Mathematics 84 (2007) 1871–1880. [6] B. Moore, Principal component analysis in linear systems: controllability, observability, and model reduction, IEEE Transactions on Automatic Control 26 (1981) 17–32. [7] Y.Q. Lin, L. Bao, Y.M. Wei, Order reduction of bilinear MIMO dynamical systems using new block Krylov subspaces, Journal of Computational and Applied Mathematics 58 (2009) 1093–1102. [8] Y.H. Cao, J. Zhu, Z.D. Luo, I.M. Navon, Reduced-order modeling of the upper tropical pacific ocean model using proper orthogonal decomposition, Journal of Computational and Applied Mathematics 52 (2006) 1373–1386. [9] C.L. Sun, J. Hahn, Reduction of stable differential-algebraic equation systems via projections and system identification, Journal of Process Control 15 (2005) 639–650. [10] F. Ebert, A note on POD model reduction methods for DAEs, Mathematical and Computer Modelling of Dynamical Systems 16 (2010) 115–131. [11] Y.L. Jiang, R.M.M. Chen, O. Wing, Convergence analysis of waveform relaxation for nonlinear differentialalgebraic equations of index one, IEEE Transactions on Circuits and Systems Part I: Fundamental Theory and Applications 47 (2000) 1639–1645. [12] Y.L. Jiang, A general approach to waveform relaxation solutions of nonlinear differential-algebraic equations: the continuous-time and discrete-time cases, IEEE Circuits and Systems: Regular Papers 51 (2004) 1770–1780. [13] Y.L. Jiang, R.M.M. Chen, Computing periodic solutions of linear differential-algebraic equations by waveform relaxation, Mathematical and Computer 74 (2005) 781–804. [14] R. Riaza, Differential-Algebraic Systems: Analytical Aspects and Circuit Applications, World Scientific, Hackensack, London, 2008. [15] F.R. Gantmacher, The Theory of Matrices, Chelsea, New York, 1959. [16] C.Y. Chen, Y.L. Jiang, H.B. Chen, An e model-order reduction approach for differential-algebraic equation systems, Mathematical and Computer Modelling of Dynamical Systems 18 (2012) 223–241. [17] B. Salimbahrami, B. Lohmann, Order reduction of large scale second-order systems using Krylov subspace methods, Linear Algebra Applications 415 (2006) 385–405. [18] R.W. Freund, SPRIM: Structure-Preserving Reduced-Order Interconnect Macromodeling, in: Proceedings of the International Conference on Computer Aided Design (ICCAD), IEEE Computer Society, Los Alamitos, California, 2004, pp. 80–87. [19] L. Giraud, J. Langou, M. Rozloznik, The loss of orthogonality in the Gram–Schmidt orthogonalization process, Computers and Mathematics with Applications 50 (2005) 1069–1075. [20] A.C. Antoulas, Approximation of Large-Scale Dynamical Systems, SIAM, Philadelphia, PA, 2005. [21] A.I. Klimushchev, N.N. Krasovskii, Uniform asymptotic stability of systems of differential equations with small parameter in the derivative terms, Journal of Applied Mathematics and Mechanics 25 (1961) 1011–1025. [22] W.H.A. Schilders, H.A. van der Vorst, J. Rommes, Model Order Reduction: Theory, Research Aspects and Applications, Springer-Verlag, Berlin, 2008. [23] R.W. Freund, Krylov-subspace methods for reduced-order modeling in circuit simulation, Journal of Computational and Applied Mathematics 123 (2000) 395–421.