Structure preserving model order reduction of a class of second-order descriptor systems via balanced truncation

JID:APNUM AID:3724 /FLA

[m3G; v1.261; Prn:18/12/2019; 9:15] P.1 (1-14)

Applied Numerical Mathematics ••• (••••) •••–•••

Contents lists available at ScienceDirect

Applied Numerical Mathematics www.elsevier.com/locate/apnum

Structure preserving model order reduction of a class of second-order descriptor systems via balanced truncation M. Monir Uddin Department of Mathematics and Physics, North South University, Dhaka-1229, Bangladesh

a r t i c l e

i n f o

Article history: Received 27 August 2019 Received in revised form 16 November 2019 Accepted 10 December 2019 Available online xxxx Keywords: Second-order index-3 descriptor systems Sparsity Structure preserving model reduction Balanced truncation Low-rank alternating direction implicit iterations

a b s t r a c t Large sparse second-order index-3 descriptor system arises in various disciplines of science and engineering including constraint mechanics, mechatronics (where mechanical and electrical elements are coupled) and circuit designs. Simulation, controller design and design optimization are some applications of such models. These tasks become challenging when the dimension of the system is high. This paper discusses a method to obtain a reduced second-order model from a large sparse second-order index-3 system using the Balanced Truncation. For this purpose, the low-rank alternating direction implicit iteration is modified to solve the Lyapunov equations of the index-3 structure system efficiently in an implicit way. Numerical resultants are discussed to show the reflectivity and efficiency of the techniques. © 2019 Published by Elsevier B.V. on behalf of IMACS.

1. Introduction In classical mechanics or multibody dynamics, a linearized holonomically constraint1 equation of motion [11] can often be represented as

M 1 ξ¨ (t ) + D 1 ξ˙ (t ) + K 1 ξ(t ) + G 1T η(t ) = H 1 u (t ), G 1 η(t ) = 0, y (t ) = L 1 ξ(t ),

(1a) (1b) (1c)

where ξ ∈ Rnξ is a vector with position coordinates related to the degree of freedom of the individual mass and η ∈ Rnη is a vector with nη (nη < nξ ) unknown parameters, the so-called Lagrange multipliers [15]. Moreover, the matrices M 1 , K 1 , D 1 ∈ Rnξ ×nξ are possibly large and sparse and known as mass, stiffness, damping respectively, and G 1 ∈ Rnη ×nξ is known as constraint matrix. Furthermore, H 1 ∈ Rnξ ×m is the input matrix corresponding to the input vector u (t ) ∈ Rm , and L 1 ∈ Rq×nξ is the output matrix associated to the output vector y (t ) ∈ Rq . In the constrained equations of motion the constraints i.e., the equations in (1b) define a manifold of free motion and the constraint forces are responsible for the constraint to be satisfied. It can be shown that the constraint forces are orthogonal to the manifold. Therefore, in Equation (1a), G 1T η(t ) appears with nη unknown parameters η .

E-mail address: [email protected]. When constraints are imposed in the position label (as seen in Equation (1b) then they are called holonomic. Otherwise, the constraints are called non-holonomic. 1

https://doi.org/10.1016/j.apnum.2019.12.010 0168-9274/© 2019 Published by Elsevier B.V. on behalf of IMACS.

JID:APNUM AID:3724 /FLA

2

[m3G; v1.261; Prn:18/12/2019; 9:15] P.2 (1-14)

M. Monir Uddin / Applied Numerical Mathematics ••• (••••) •••–•••

System in (1) is called second-order index-3 descriptor system due to the analogy of first-order conversion [26]. Such structure system arises in a large variety of applications, in particular the modeling of constrained vibrational systems, multibody systems, or electrical networks (see e.g. [11,20]). In many applications, models with the above structure can become very large and complex. For instance, in vibrational analysis, simulations of the models are often generated by the finite element method. If the grid resolution is very fine, as many geometrical details need to be resolved, then the size of the system gets very large. When this happens, it is very expensive to simulate, control and optimize the system. Therefore, we want to reduce the complexity of the model by applying model order reduction (MOR), i.e. we seek for an approximation to the original model that well-approximates the behavior of the original model but is much faster to evaluate. For systems governed by ordinary differential equations (ODEs), there already exists a well-developed theory. The two most frequently applied modern MOR methods are balanced truncation (BT) [17] and rational interpolation of the transfer function by the iterative rational Krylov algorithm (IRKA) [12]. Both approaches have been extended to general descriptor systems [22,13]. In that case, one has to apply certain projections to the model in order to decouple the differential and the algebraic equations. However, the computation and application of the projectors can be very demanding in terms of computational complexity, memory requirements and robustness issues. An easier approach exists, where the projector has to be constructed explicitly onto the differential part of the system from the given problem data and thus eliminating the algebraic equations [14,13]. Then, by exploiting certain technical properties of the projector, one can directly apply the methods for ODE systems to the projected system. This procedure has however been studied for structured first order index-2 systems in both the BT and IRKA based model reduction method. In case of structured systems of index-3 the BT and IRKA have been discussed in [26] and [1] respectively. These days, structure preserving MOR (SPMOR) of second-order systems received a lot of attention. Note that in some applications, preserving the structure is essential if the simulation of a system is required, using a software designed for a particular system. The SPMOR has been discussed in several research articles. See, e.g., [19] and the references therein. Recently, the idea of SPMOR is extended to the second-order index-1 descriptor systems [5,9]. Although, authors in [1] have investigated the SPMOR of second-order index-3 system via IRKA, until now, and to the best of the author’s knowledge, there is no such investigation for the BT technique which is addressed in this paper. This paper discusses how to perform the SPMOR via balanced truncation of the second-order index-3 system (1) without forming an ODE system explicitly. One should know that although the BT based SPMOR generates reduce-order balanced systems, in general, it does not guarantee the stability of the system and does not have a global error bound, see, e.g., [19,6] for details. Another drawback of this method is to solve the two Lyapunov equations for determining the system Gramian factors for the derivation of the balancing and truncating transformations. For a large sparse LTI system low-rank alternating direction implicit (LR-ADI) [4,18] iteration is one of the efficient methods to compute these Gramian factors. This paper also discusses how to solve the projected Lyapunov equations using the LR-ADI iteration. The Lyapunov equations can also be solved by other iterative methods e.g., Krylov based projection methods [21] and Sign Function methods [3] for computing the low-rank Gramian factors. This paper is particularly interested on the LR-ADI method since it computes the real Gramian factors [7] and the method has residual-factor based stopping criterion [6] to terminate the iteration. Note that the LR-ADI iteration for solving such projected Lyapunov equations of the underlying system was already discussed in [26]. In this work, the LR-ADI method is updated by improving the efficiency in solving the linear system at each iteration step. Moreover instead of heuristic shifts, this paper uses adaptive shifts for a better convergence of the LR-ADI iterations. These two approaches significantly reduce the computational time of the method to solve the Lyapunov equations. Finally the theoretical results are investigated by numerical experiments and the results are compared (if possible) with the existing algorithms. The rest of this article is organized as follows. The fundamental idea of SPMOR via balanced truncation and the LR-ADI for standard second-order system are repeated briefly in Section 2. Section 3 presents the ideas of the SPMOR of secondorder index-3 descriptor system. Section 4 consists the LR-ADI algorithm for the implicit solution of the projected Lyapunov equations and also the techniques for selecting adaptive shift parameters for the model considered in this paper. The subsequent section illustrates numerical results. Section 6 provides the conclusive remarks. 2. Balanced truncation for standard second-order systems This section briefly discusses the BT method for standard second-order LTI systems from previous literature. The ideas of this section will be used for index-3 system in the next sections. Consider a second-order LTI system

M x¨ (t ) + D x˙ (t ) + K x(t ) = H u (t ),

y (t )= Lx(t ),

(2)

where M , D and K are nonsingular matrices, and x(t ) ∈ Rn is a vector with the states. Unlike the index-3 system, here the matrix M is invertible. Therefore, the system is called standard system. The input-output relation of this dynamical system in the Laplace domain can be represented as

(s) = L (s2 M + sD + K )−1 H

(3)

where s ∈ C which is known as transfer function matrix of the system. To perform the balanced truncation of this system one should solve the following two Lyapunov equations [10,19]

JID:APNUM AID:3724 /FLA

[m3G; v1.261; Prn:18/12/2019; 9:15] P.3 (1-14)

M. Monir Uddin / Applied Numerical Mathematics ••• (••••) •••–•••

A P E T + E P AT = −B B T T

T

and

(4a)

T

A Q E + E Q A = −C C , where

 E=

M 0





0 0 , A= M −K

3

(4b)



M ,B= −D



0 H







and C = L

0 .

The solutions P ∈ R2n×2n and Q ∈ R2n×2n of the Lyapunov equations (4a) and (4b) are symmetric positive definite (SPD) matrices and respectively known as the controllability and observability Gramians of the system. The Lyapunov equations can be solved efficiently by using the low-rank alternating direction implicit [4] (LR-ADI) iteration. Solving the Lyapunov equations (4a) and (4b), the method can compute the approximate low-rank Gramian factors



R=

Rp Rv





and

L=

Lp Lv



(5)

respectively, where P ≈ R R T and Q ≈ LL T . In both cases the computed R ∈ R2n×k and L ∈ R2n×k (k  n) are the thin rectangular matrices and the n × k matrices R v , R p , L v and L p are known as the controllability-velocity, controllability-position, observability-velocity and observability-position Gramian factors respectively. Once we have R p and L p , the balancing and truncating transformation can be formed using the singular value decomposition (SVD)



T R Tp M L p = U pp  pp V pp = U pp ,1

U pp ,2



   pp ,1

 pp ,2

T V pp ,1



T V pp ,2

,

and defining −1

T l := L p U pp ,1  pp2,1 ,

−1

T r := R p V pp ,1  pp2,1 .

(6)

Here U pp ,1 and V pp ,1 are composed of the leading k columns of U pp and V pp respectively and  pp ,1 is the first k × k block of the matrix  pp . Now, the reduced order model is

ˆ x¨ˆ (t ) + Dˆ x˙ˆ (t ) + Kˆ xˆ (t ) = Hˆ u (t ), M

yˆ (t ) = Lˆ xˆ (t ),

(7)

where xˆ (t ) ∈ Rk (k  n) is formed by constructing the reduced matrices as

ˆ = T T M T r , Dˆ = T T D T r , Kˆ = T T K T r , Hˆ = T T H and Lˆ = LT r . M l l l l

(8)

The goal is to minimize  y − yˆ H∞ sufficiently which can be measured alternatively by

ˆ H∞ ,  − 

(9)

ˆ + s Dˆ + Kˆ )−1 Hˆ is the transfer function matrix of the reduced system and .H∞ denotes the H∞ -norm ˆ s) = Lˆ (s2 M where ( of the system. Definition 1. The H∞ -norm of the stable system (2) is defined by

H∞ = sup σmax ((jω)),

(10)

ω∈R

where

σmax denotes the maximum singular value of (jω).

The balanced truncation of the first-order system has a global error bound [2] which can be defined priori by a user. In the structure preserving model reduction there is no such an error bound which will be reflected in the numerical results. However, one may monitor the ratio of the entries in  pp and can determine the dimension of the reduced model once [6]

σ j, pp ≤  (truncation-tolerance) σ1, pp

(11)

is since it uses  fulfilled. The above procedure of constructing the reduced model is called position-position  (pp)  balancing R p , L p for constructing the balancing and truncating transformations. Similarly when R p , L v , R v , L p and { R v , L v } are used for constructing the balancing and truncating transformations, the procedure of model reduction are called positionvelocity (pv), velocity-position (vp) and velocity-velocity (vv) balancing respectively. The fundamental drawback of the balancing based model reduction is that it computes two Gramian factors by solving two respective Lyapunov equations i.e., (4). Among several approaches, recently the LR-ADI method is considered as one of the efficient and most commonly used approach for the large-scale dynamical systems. Readers are referred to [4] and

JID:APNUM AID:3724 /FLA

[m3G; v1.261; Prn:18/12/2019; 9:15] P.4 (1-14)

M. Monir Uddin / Applied Numerical Mathematics ••• (••••) •••–•••

4

Algorithm 1: LR-ADI iteration for solving APE T + EPA T = −BB T . J

Input : E , A, B , {μi }i =1 and imax .

Output: R = Z i , such that P ≈ RR T . 1 W 0 = B. 2 Z 0 = [ ], i = 1. 3 while  W iT−1 W i −1  ≥ τ (tolerance) or i ≤ i max do 4 5 6 7

Compute V i = (A + μi E )−1 W i −1 . if Im (μi )= 0 then  √ −2μi V i . Z i = Z i −1 W i = W i −1 − 2μi E V i .

8

else Re (μi ) γ = −2 Re (μi ), δ = Im (μi ) ,

Z i = 2γ (Re ( V i ) + δ Im ( V i ))



Z¯ i = 2γ (δ 2 + 1) Im ( V i )

9 10 11





Z i +1 = Z i −1 Z i Z¯ i , W i +1 = W i −1 + 2γ E (Re ( V i ) + δ Im ( V i )). i =i+1

12 13 14

i =i+1

15

the references therein for details on the LR-ADI approach. The recent updates of this prominent method are obtained in [25]. This paper follows [25, Algorithm 5] to solve the Lyapunov equations in (4). For a convenience this algorithm is again presented here (i.e., Algorithm 1) as LR-ADI iteration. Note that in [25] the name of this algorithm was G-LRCF-ADI. This algorithm can solve Lyapunov equations (4a) and (4b) by replacing inputs (E , A, B ) by ( E , A , B ) for the controllability Gramian factor (R) and ( E T , A T , C T ) for the observability Gramian factor (L) respectively. The algorithm is provided a set of J shift parameters {μi }i =1 . If the number of iterations imax is greater than J , then the shift-parameters are used cyclic way. 3. Balanced truncation for second-order index-3 descriptor systems This section turns to the second-order index-3 descriptor system (1). Convert the index-3 system into an index-0 system form

 M 1 T ξ¨ (t ) +  D 1 T ξ˙ (t ) +  K 1 T ξ(t ) =  H 1 u (t ),

(12a)

T

y (t ) = L 1  ξ(t ),

(12b)

where

 := I nξ − G 1T (G 1 M 1−1 G 1T )−1 G 1 M 1−1 ,

(13)

is called projector onto the null-space of G 1 . A complete analogy of such formulation is available in [26],  and there the author also shows that the projector  fulfills the following properties.  M 1 = M 1  T , Null () = Range G 1T , Range () =

Null G 1 M 1−1



and

G 1 ξ = 0 iff  T ξ = ξ.

(14)

The dynamical system (12) still has redundant equations due to the singularity of . The redundant part can be avoided by applying the thin singular value decomposition (SVD):

  T   1 0  V1  = U V T = U 1 U 2 = U 1 1 V 1T = l rT , T 0

1

0

(15)

V2

1

where l = U 1 12 , r = V 1 12 and U 1 , V 1 ∈ Rnξ ×nm consist of the corresponding leading nm columns of U , V ∈ Rnξ ×nξ . Moreover, l and r satisfy

lT r = I nm .

(16)

Inserting the decomposition of  from (15) into (12) and considering ξ˜ (t ) = l ξ(t ), the resulting dynamical system leads to T

¨˜ t ) + T D ξ( ˙˜ T T ˜ rT M 1 r ξ( 1 r t ) + r K 1 r ξ (t ) = r H 1 u (t ), r y (t ) = L 1 r ξ˜ (t ).

(17a) (17b)

JID:APNUM AID:3724 /FLA

[m3G; v1.261; Prn:18/12/2019; 9:15] P.5 (1-14)

M. Monir Uddin / Applied Numerical Mathematics ••• (••••) •••–•••

5

This system is now a standard second-order system as described by equation (2). This system can be seen as the system (12) with the redundant equations being removed through the r projection. Note that the coefficient matrices of (17) are dense when compared to (1). Therefore, for the large-scale index-3 system explicit computation of (17) is forbidden. The dynamical systems (2), (12) and (17) are equivalent in the sense that they are the different realizations of the same transfer function. Moreover, their finite spectra are the same which has been proven in the sequel. Equivalent finite spectra. The quadratic matrix polynomial associated with the index-3 DAE system (1) is defined by



Q (λ) = λ2



M1 0











0 D1 +λ 0 0

M

 D









0 K + 1 0 G1



G 1T , 0

 K

(18)



where λ ∈ C . Although Q (λ) is regular, due to the singularity of M, it contains some infinite eigenvalues as well. If the ˜ where n˜ = nξ + nμ , then Q (λ) has r finite and 2n˜ − r infinite eigenvalues [23]. The quadratic degree of det ( Q (λ)) is r < 2n, matrix polynomials corresponding to the systems (12) and (17) are

Q˜ (λ) = λ2  M 1  T + λ D 1  T +  K 1  T

(19)

Q¯ (λ) = λ2 rT M 1 r + λ rT D 1 r + rT K 1 r ,

(20)

and

respectively. We know that M 1 r ∈ R is nonsingular and Q¯ is regular. Therefore, all of the eigenvalues of the  polynomial Q¯ (λ) are finite [23]. The degree of det Q¯ (λ) is 2(nξ − nμ ). Hence the number of finite eigenvalues of Q¯ (λ) is exactly 2(nξ − nμ ) and the number of infinite eigenvalues is 2n˜ − 2(nξ − nμ ) = 2nξ + 2nμ − 2nξ + 2nμ = 4nμ infinite eigenvalues. By applying the appropriate projectors onto the index-3 DAE system (1), we can preserve all the finite eigenvalues of the system. The following theorem proves that both the original and the projected systems have the same finite eigenvalues. nξ ×nξ

rT

Theorem 3.1. Let us consider the matrix polynomials Q (λ) and Q¯ (λ) respectively defined in (18) and (20). An eigenvalue λ1 is a



finite eigenvalue of Q (λ) with corresponding eigenvector v 1T eigenvector v˜ 1 where v˜ 1 = lT v 1 with l defined in (16).

v 2T

T

if and only if λ1 is an eigenvalue of Q¯ (λ) with corresponding



Proof. Suppose λ1 is a finite eigenvalue of Q (λ) corresponding to the eigenvector v 1T problem of the matrix polynomial (18) is





λ21

M1 0



0 + λ1 0



D1 0





0 K1 + 0 G1

G 1T 0

 

v1 v2



v 2T

T

. Then the quadratic eigenvalue

 

=

0 . 0

(21)

The last equation of (21) gives G 1 v 1 = 0 which implies v 1 is in the null-space of G 1 . By using (14), we obtain  T v 1 = v 1 . Plugging  T v 1 = v 1 into the first equation of (21) and multiplying the resulting equation from the left by  and taking into account that G 1T = 0 leads to

(λ21  M 1 T + λ1  D 1 T +  K 1 T ) v 1 = 0,

(22)

which is the eigenvalue problem for the matrix polynomial Q˜ (λ), see equation (19). Applying the decompositions of  as defined above to (15) and using v˜ 1 = lT v 1 we obtain

l (λ21 rT M 1 r + λ1 rT D 1 r + rT K 1 r ) v˜ 1 = 0. Multiplying by r from the left and using (16) yields

(λ21 rT M 1 r + λ1 rT D 1 r + rT K 1 r ) v˜ 1 = 0,

(23)

which is the eigenvalue problem of the matrix polynomial (19) with λ1 as an eigenvalue. Conversely, we want to demon-



T

strate that if v˜ 1 is an eigenvector of Q¯ (λ) to the corresponding eigenvalue λ1 i.e. equation (23) holds, then v 1T v 2T is an eigenvector of Q (λ) with the same eigenvalue. Again, plugging v˜ 1 = lT v 1 into (23) and multiplying the resulting equation by l from the left we obtain

(λ21  M 1 T + λ1  D 1 T +  K 1 T ) v 1 = 0. Since the projector  satisfies  T v 1 = v 1 , (24) gives

(24)

JID:APNUM AID:3724 /FLA

[m3G; v1.261; Prn:18/12/2019; 9:15] P.6 (1-14)

M. Monir Uddin / Applied Numerical Mathematics ••• (••••) •••–•••

6

(λ21 M 1 v 1 + λ1 D 1 v 1 + K 1 v 1 ) = 0,





which means that λ21 M 1 v 1 + λ1 D 1 v 1 + K 1 v 1 is in the null-space of . We know that Null () = Range G 1T . Therefore, there exists a vector v 2 such that

λ21 M 1 v 1 + λ1 D 1 v 1 + K 1 v 1 = −G 1T v 2 , which implies

λ21 M 1 v 1 + λ1 D 1 v 1 + K 1 v 1 + G 1T v 2 = 0.

(25)

Again if  v 1 = v 1 yields T

G1 v1 = 0

(26)

Equations (25) and (26) yield (21).

2

Once the index-3 system (1) is converted into the index-0 system (17) then the balanced truncation can be applied to the converted system. However, the converted system as described by equation (17), is dense and therefore, considering the computational complexity forming the index-0 system is very costly. In the following an implementation of a model reduction technique which does not require explicit forming of the index-0 system is proposed. Model reduction of the projected system. The balancing idea introduced in Section 2 can be applied to the projected system (17). However and as already mentioned, this is infeasible for a large-scale system. Therefore, the BT technique can be applied to the equivalent system (12) instead. For this purpose, recalling the discussion in Section 2, the projected Lyapunov equations:

˜ P˜ E˜ T + E˜ P˜ A˜ T = − B˜ B˜ T , A

(27a)

˜ T Q˜ E˜ + E˜ T Q˜ A˜ = −C˜ T C˜ , A

(27b)

where

E˜ =



I nξ 0

0  M 1 T



˜= A

,

B˜ =



0

 

0 − K 1 T



C˜ = L 1  T

and

H1

I nξ − D 1 T

 , 

0 ,

are to be solved. The LR-ADI iteration can compute the low-rank controllability Gramian factor R˜ ( R˜ R˜ T ≈ P˜ ) and observability Gramian factor L˜ ( L˜ L˜ T ≈ Q˜ ) by solving the Lyapunov equations (27a) and (27b) respectively. Due to the block structure of P˜ and Q˜ , the low-rank Gramian factors can be partitioned as

R˜ =





R˜ p , R˜ v

L˜ =





L˜ p , L˜ v

where R˜ p , L˜ p and R˜ v , L˜ v denote the low-rank position, velocity controllability and observability Gramian factors respectively. Let us consider the Gramian factors R˜ p and L˜ p to compute the thin SVD

  ˜  ˜L Tp M 1 R˜ p = U˜ pp ,1 U˜ pp ,1  pp ,1 0

0

˜ pp ,2 

 ˜ T V

pp ,1 T V˜ pp ,2



,

(28)

and construct the balancing and truncating transformations −1

˜ 2 , T˜ l = L˜ p U˜ pp ,1  pp ,1

−1

˜ 2 . T˜ r = R˜ p U˜ pp ,1  pp ,1

(29)

Now applying T˜ l and T˜ r to the system (12) we can construct the reduced model

¨

˙

ˆ 1 ξ( ˆ t ) + Dˆ 1 ξ( ˆ t ) + Kˆ 1 ξˆ (t ) = Hˆ 1 u (t ), M yˆ (t ) = Lˆ 1 ξˆ (t ).

(30)

Due to the properties of the projector  T T˜ l = T˜ l and  T T˜ r = T˜ r , the coefficient matrices in (30) can be obtained as

ˆ 1 = T˜ T M 1 T˜ r , M l

ˆ 1 = T˜ T D 1 T˜ r , D l

Kˆ 1 = T˜ lT M 1 T˜ r ,

ˆ 1 = T˜ T H 1 , H l

Lˆ 1 = L 1 T˜ r ,

(31)

JID:APNUM AID:3724 /FLA

[m3G; v1.261; Prn:18/12/2019; 9:15] P.7 (1-14)

M. Monir Uddin / Applied Numerical Mathematics ••• (••••) •••–•••

7

Algorithm 2: Balanced truncation for second-order index-3 system. Input : M 1 , D 1 , K 1 , H 1 , L 1 . ˆ 1 , Dˆ 1 , Kˆ 1 , Hˆ 1 , Lˆ 1 . Output: M ˜ p and L˜ p . 1 Solve the Lyapunov equations (27a) to compute R 2 Compute the balancing and truncating transformations as in (29).

ˆ 1 , Dˆ 1 , Kˆ 1 , Hˆ 1 , Lˆ 1 as in (31). 3 Construct M

which prevents from constructing the projected system (12). The process of obtaining a reduced model by using a pair of low-rank controllability and observability position Gramian factors, i.e. ( R˜ p , L˜ p ) is summarized in Algorithm 2. The constructed reduced model via ( R˜ p , L˜ p ) is called position-position (pp) balancing. Likewise, the reduced models are called velocity-velocity (vv), velocity-position (vp), and position-velocity (pv) balancing if we use the pairs ( R˜ v , L˜ v ), ( R˜ v , L˜ p ) and ( R˜ p , L˜ v ) respectively. In Algorithm 2 the most challenging task is to solve the projected Lyapunov equations (27) for computing R˜ p and L˜ p , which can be solved by applying Algorithm 1. Now the remaining question is whether one can avoid the projector  as given by equation (13) for the solution of the Lyapunov equations. The answer to this question is given in the following section. 4. Solution of the projected Lyapunov equations The previous section pointed out that in order to implement the balancing based MOR for the second-order index-3 descriptor system (1), two projected Lyapunov equations as defined in (27) have to be solved for computing the low-rank ˜ This section discusses how to apply the LR-ADI iteration i.e. Algorithm 1 to solve such projected Gramian factors R˜ and L. Lyapunov equations efficiently by first focusing on the solution of the controllability Lyapunov equation (27a). Following subsections provide a step by step procedure to update Algorithm 1 for solving the underlying Lyapunov equation. 4.1. Initial residual factor Step 1 in Algorithm 1 computes the initial Lyapunov residual factor W 0 . For solving the Lyapunov equation (27a) consider ˜ 0 which can be partitioned as the initial residual factor W

˜ 0 = B˜ = W



0 H1



 =

˜ W 0

(1 )

˜ W 0

(2 )



(32)

,

˜ (1) = 0 and W ˜ (2) =  H 1 . The matrix  of the linear system where W 0 0



M1 G1

G 1T 0



    H1 =  0

(33)

gives  = M 1−1  H 1 . Now multiplying  by M 1 ,  H 1 can be obtained by avoiding the explicit computation of  while computing the initial Lyapunov residual factor. 4.2. Efficient solution of the linear systems Solving a linear system at each iteration is the most expensive task in the LR-ADI iteration. Applying Algorithm 1 to equations (27a) for the controllability Gramian, at each i-th iteration step one needs to solve the linear system

˜ i −1 , ( A˜ + μi E˜ ) V i = W



˜ i −1 is the ADI residual factor computed from the (i − 1)-st iteration. Now partitioning W ˜ i −1 as W ˜ i −1 = where W equation (34) can be written as





μi I n ξ

I nξ

− K 1 T

− D 1 T + μi  M 1 T

(1 )

Vi

(2 )

Vi



 =

˜ W i −1 (1 )

˜ W i −1 (2 )

(34) 

˜ (1) W i −1 ˜ (2) W i −1

,

 ,

(35)

which gives the linear systems

μi V i(1) + V i(2) = W˜ i(−1)1 , T

(1 )

K1 V i

T

T

(2 )

+ ( D 1  + μi  M 1  ) V i

˜ (2 ) . =W i −1

(36) (37)

JID:APNUM AID:3724 /FLA

[m3G; v1.261; Prn:18/12/2019; 9:15] P.8 (1-14)

M. Monir Uddin / Applied Numerical Mathematics ••• (••••) •••–•••

8

(1)

(1)

(2)

Now,  T V i = V i ,  T V i fore, equation (37) yields (1 )

( K 1 V i (1)

i.e. K 1 V i

(2 )

+ D1 V i

(2)

= Vi

˜ (2) =  W ˜ (2) . Thereand it can be shown that residual factor in every step satisfies W i −1 i −1 (2 )

+ μi M 1 V i

˜ ( 2 ) ) = 0, −W i −1

(38)



T

˜ (2) is in the null-space of . Since Null () = Range G , there exists a , such that + D 1 V i(2) + μi M 1 V i(2) − W 1 i −1 (1 )

(2 )

˜ (2) = −G 1T , −W i −1

(2 )

˜ + G 1T  = W . i −1

(2 )

+ μi M 1 V i

(2 )

+ μi M 1 V i

+ D1 V i

K1 V i

which implies (1 )

+ D1 V i

K1 V i

(2 )

(39)

Equation (36) gives (1 )

1

=

Vi

μi

˜ (1 ) − V (2 ) ) (W i −1 i

(40)

Plugging this into (39) yields

K1

1

μi

˜ (1) − V (2) ) + D 1 V (2) + μi M 1 V (2) + G 1T  = W ˜ (2 ) , (W i −1 i i i i −1

which implies

˜ (2 ) − K 1 W ˜ (1 ) . (μ2i M 1 + μi D 1 − K 1 ) V i(2) + G 1T  = μi W i −1 i −1 (2)

= V i , then V i

(2 )

= 0.

With  T V i

G1 V i

(2)

(2)

is in the null-space of G 1 , i.e.

(42)

Combining (41) and (42) yields



(41)

μ2i M 1 + μi D 1 − K 1 G 1T G1



(2 )



0





μi W˜ i(−2)1 − K 1 W˜ i(−1)1 . =

Vi

(43)

0



(2)

Therefore, instead of solving (35) (i.e. a projected linear system), Equation (43) can be solved for V i



(40). The vector

˜ W i −1 (1)



(1)

and then V i

from

(for the MIMO case it is a matrix) is updated after each iteration computed from the Lyapunov ˜ (2) W i −1 residual factor of the previous step. In each iteration the Lyapunov residual factor can be computed by (see, Step 7 in Algorithm 1)

˜i=W ˜ i −1 − 2 Re (μi ) E˜ V i , W which can be partitioned as



˜ W i

(1 )

˜ W i

(2 )



 =

(1 )

˜ W i −1 (2 )

 =

˜ W i −1





I − 2 Re (μi ) nξ 0

0  M 1 T



˜ W i −1 − 2 Re (μi ) V i (1 )

(1 )

˜ W − 2 Re (μi ) M 1 T V i i −1 (2 )



(2 )

(1 )

Vi

(2 )

Vi

 .

. (2)

Exploiting the properties of , i.e.  M 1 = M 1  T and  T V i

˜ W i

˜ (1) − 2 Re (μi ) V (1) , =W i −1 i ˜ (2 ) = W ˜ (2) − 2 Re (μi ) M 1 V (2) . W

(2)

= V i , the above equation results in

(1 )

i −1

i

(44)

i

If the two consecutive shift parameters are complex conjugates of each other, that is when {μi , μi +1 := μi }, then for Step 13 in Algorithm 1)

˜ i +1 = W ˜ i −1 − 4 Re (μi ) E˜ (Re ( V i ) + δ Im ( V i )) , W where δ =

Re (μi ) Im (μi )

gives

JID:APNUM AID:3724 /FLA

[m3G; v1.261; Prn:18/12/2019; 9:15] P.9 (1-14)

M. Monir Uddin / Applied Numerical Mathematics ••• (••••) •••–•••

9

Algorithm 3: LR-ADI for the second-order index-3 systems. J

Input : M 1 , D 1 , K 1 , G 1 , H 1 , {μi }i =1 and imax .

Output: R˜ = Z i , such that P˜ ≈ R˜ R˜ T . 1 Set Z 0 = [ ].

˜ 2 i = 1 and W 0 = 0. ˜ (2) = M 1 . 3 Solve (33) for  and compute W 0 ˜ (1) T W ˜ (2) T W ˜ (1) + W ˜ (2)  ≥ τ (tolerance) and i ≤ i max do 4 while  W 0 0 0 0 (1)

5

(2)

Solve (43) for V i



(1)

and then compute V i

(1) T Vi

(2) T Vi

T

6

Compute V i =

7 8

if Im (μi )= 0 then  √ −2μi Re ( V i ) , Z i = Z i −1

˜ W i

(1)

9 10

˜ =W − 2μi V i , i −1 (1)

(2)

˜ =W − 2μi M 1 V i . i −1 (2)

(2)

else

11

γ = −2 Re (μi ), δ =

12

Z i +1 =

13

V i +1 = Re ( V i



2γ (Re ( V i ) + δ Im ( V i ))

(1)

˜ (1) W i +1

=

i=i+1

15

Re (μi ) , Im (μi )



Z i −1

(1)

14

16

˜ W i

(1)

from (40).

.

˜ (1) W i −1

(1)



2γ

(2)



 (δ 2 + 1) Im ( V i ) ,

(2)

(2)

) + δ Im ( V i ), V i +1 = Re ( V i ) + δ Im ( V i ), (1) ˜ (2) + 2γ M 1 V (2) . ˜ (2) = W + 2γ V , W i +1

i +1

i −1

i +1

i =i+1



˜ W i +1 (1 )



 =

˜ W i +1 (2 )

˜ W i −1 (1 )

˜ W i −1 (2 )

 =





I − 4 Re (μi ) nξ 0

0  M 1 T



(1 ) (2 )

(1 )

(2 )



(1 )

˜ W − 4 Re (μi ) M 1 T (Re ( V i ) + δ Im ( V i )) i −1 (2 )



Re ( V i ) + δ Im ( V i )

˜ W − 4 Re (μi )(Re ( V i ) + δ Im ( V i )) i −1 (1 )

(1 )

Re ( V i ) + δ Im ( V i )

(2 )

(2 )

.

Again, following the properties of  the above relation results in

˜ ˜ W =W − 4 Re (μi )(Re ( V i ) + δ Im ( V i )), i +1 i −1 (1 )

(1 )

(1 )

(1 )

(45)

˜ ˜ W =W − 4 Re (μi ) M 1 (Re ( V i ) + δ Im ( V i )). i +1 i −1 (2 )

(2 )

(2 )

(2 )

Therefore, the Lyapunov residual can be computed without computing the projector  explicitly. The described procedure to compute the low-rank controllability Gramian factor by solving the controllability Lyapunov equation (27a) is stated in Algorithm 3. 4.3. Solution of observability Lyapunov equation Algorithm 3 can solve the projected observability Lyapunov equation (27a) considering a few changes. First, in the input ˜ (1) =  L T , W ˜ (2) = 0. Replacing H 1 by L T , solve the linear matrices replace H 1 by L 1 . The initial residual factors are W 1 0 0 (1)

system (33) for  to compute W 0



= M 1 . In Step 5, solve the linear system

 T

μ2i M 1T − μi D 1T + K 1T G 1 G1

0

(2 )

Vi







μi W˜ i(−2)1 − W˜ i(−1)1 =



0

(2)

to compute V i , then compute (1 )

Vi



=

1

μi

for V i = V (1) i

T

˜ (1) +  K 1T V (2) ) (W i −1 i (2) T

Vi

T

. Note that  K 1T can be computed implicitly as  H 1 discussed in subsection 4.1. Applying these

changes in Algorithm 3, compute L˜ = Z i such that Q˜ ≈ L˜ L˜ T . 4.4. ADI shift parameter selection Algorithm 3 depends on certain shift parameters that are crucial for fast convergence of the method. Recently two types of ADI shift parameters such as Penzl’s heuristic [18] and adaptive [8] shift parameters are often used for solving large-scale Lyapunov equations. The Penzl’s approach was already implemented in [26] for the system (1). This paper investigates the

JID:APNUM AID:3724 /FLA

[m3G; v1.261; Prn:18/12/2019; 9:15] P.10 (1-14)

M. Monir Uddin / Applied Numerical Mathematics ••• (••••) •••–•••

10

Fig. 1. A damped spring-mass system (DSMS) with holonomic constraint ([16]).

second type of shift computation strategy, i.e. the adaptive approach which is rather simple and more efficient. In this case the shift parameters are associated with the eigenvalues of the quadratic matrix polynomial (19). From the deliberation in Section 3 it is already known that the matrix polynomial (19) includes all of the finite eigenvalues of the index-3 system in (1). For the initialization of the shifts, the matrix polynomial (19) is projected onto the range of a random thin rectangular matrix Bˇ ∈ Rnξ ×k , where k  nξ , yielding

U T Q˜ (λ)U = λ2 U T  M 1  T U + λU T  D 1  T U + U T  K 1  T U , (46)

 where U = Range Bˇ . Once all the initial shift parameters are set, the shifts are updated considering U with all the V i (in Algorithm 3) generated by the previous shift values. Now from the eigenvalues of (46), select some desired number of shifts by solving the so called ADI min-max problem as discussed in the Penzol’s heuristic procedure [18]. Note that for any orthogonal bases U the projector  satisfies  T U = U . As an immediate result of this fact, the projected matrix polynomial in (46) can be replaced by

Pˇ (λ) = λ2 U T M 1 U + λU T D 1 U + U T K 1 U .

(47)

5. Numerical results This section discusses some numerical experiments to show the accuracy and efficiency of the model reduction techniques discussed in the above sections with two model examples. The first example is a damped spring-mass system (DSMS) with holonomic constraint as shown in Fig. 1 which is taken from [16]. Following the discussion in [16] the i-th mass of weight mi is connected to the (i − 1)-st mass by a spring and a damper with constants ki and δi , respectively. Moreover, the first mass is connected to the last one by a rigid bar and influenced by the control u (t ). This numerical experiments takes

m1 = m2 = · · · = mξ = 100, k 1 = · · · = k ξ − 1 = κ 2 = · · · = κξ − 1 = 3 , d1 = · · · = dξ −1 = δ2 = · · · = δξ −1 = 1,

κ1 = κ ξ = 6 , δ1 = δξ = 2.

Therefore the resultant coefficient matrices are M 1 = 100 I nξ ,

⎡

−9

⎢ 3 ⎢ ⎢ K1 = ⎢ ⎢ ⎢ ⎣

3 −9 3

⎤ ⎥ ⎥ ⎥ ⎥, . −9 ⎥ ⎥ .. .. . . 3 ⎦ 3 −9 3

..

⎡

−3

⎢ 1 ⎢ ⎢ D1 = ⎢ ⎢ ⎢ ⎣

1 −3 1

⎤ ⎥ ⎥ ⎥ ⎥ . −3 ⎥ ⎥ .. .. . . 1 ⎦ 1 −3 1

..

and the constraint matrix G = [1, 0, · · · , 0, −1]. The Input matrix H 1 = e 1 and the output matrix L 1 = [e 1 , e 2 , enξ −1 ] T , where e i denotes the i −th column of the identity matrix I nξ lead the system to one input and three outputs. Considering nξ = 2000 masses a 2001 dimensional second-order index-3 system is formed. The second example considered here is a triple chain oscillator model (TCOM) from [24] with the index-3 setup described in [26]. As shown in Fig. 2 TCOM contains three chains and each of the chains is coupled to a fixed mounting by an additional damper on one end and fixed rigidly to a large mass coupling all three of them. The large mass is bound to a fixed mount by a single spring element. Each of the chains consists of nξ equal masses and of equal stiffness in the spring elements. Therefore, the model parameters are the masses m1 , m2 , m3 and the corresponding stiffnesses k1 , k2 , k3 of the three oscillator chains, the coupling mass m0 with its spring stiffness k0 , the viscosities ϑ of the additional wall-mount-dampers and the  length nξ of the single oscillator chains. The resulting system then is of order nξ = 3n1 + 1. The mass matrix M 1 = diag m1 I nξ , m2 I nx i , m3 I n1 , m0 ; the stiffness matrix K 1 and damping matrix D 1 consist of a leading block diagonal matrix (consisting of the three stiffness matrices for the three oscillator chains) and coupling terms in the

JID:APNUM AID:3724 /FLA

[m3G; v1.261; Prn:18/12/2019; 9:15] P.11 (1-14)

M. Monir Uddin / Applied Numerical Mathematics ••• (••••) •••–•••

11

Fig. 2. A triple chain oscillator model (TCOM) ([24]).

Table 1 The performances of the heuristic and adaptive shifts in the LR-ADI iteration in terms of convergence rate and computational cost. Models

DSMS TCOM

Tolerance (τ )

Heuristic shifts

τ

Iterations

10−8 10−8

Adaptive shifts

R˜

L˜

Time (sec.) ( R˜ + L˜ + μ)

109 166

122 170

6 37

Iterations R˜

L˜

Times (sec.) ˜ ( R˜ + L)

86 153

120 240

3 4

Table 2 Number of optimal shift parameters used in heuristic approach and each cycle of adaptive approach. Models

Number of shift parameters Heuristic shifts

DSMS TCOM

Adaptive shifts

R˜

L˜

R˜

L˜

25 100

25 100

10 10

20 50

last row and column at positions n1 , 2n1 , 3n1 and in the diagonal element. For the numerical experiment, the values of the variables are considered as follows:

m1 = 1, m2 = 2, m3 = 3, m0 = 10, k1 = 10, k2 = 20, k3 = 1, k0 = 50,

ϑ = 5. To have the index-3 structure, consider 5000 algebraic constraints in which each constraint is obtained by connecting two masses using a rigid bar as considered for DSMS. The constraints are chosen randomly. If nξ = 2000 then G 1 becomes a 5000 × 6001 matrix and the dimension of the second-order index-3 system becomes 11001. Input matrix H 1 ∈ R6001×1 consists of all elements with one and the output matrix L 1 = H 1T . Hence, the system is single-input-single-output (SISO). The experiments are carried out with MATLAB® R2015a (8.5.0.197613) on a board with 4×Intel®Core™i5-4460s CPU with a 2.90 GHz clock speed and 16 GB RAM. The low-rank Gramian factors are computed with Algorithm 3 for both models considering the tolerance τ = 10−8 . To implement this algorithm, the performances of both the heuristic and adaptive shifts are investigated. The results are shown in Table 1. For both models, the adaptive shift approach requires less computational time than the heuristic approach. Note that the adaptive shift selection approach is discussed in Subsection 4.4 and the heuristic shifts are selected by following [26, Section 4]. The number of appropriate optimal shift parameters selection is an important issue for better convergence of the LR-ADI. In these experiments the number of selected shift parameters for DSMS and TCOM in both approaches are shown in Table 2. In the case of heuristic approach DSMS model considers 25 optimal shifts from 40 large and 30 small magnitude Ritz-values. On the other hand for TCOM 100 optimal shifts are selected from 120 large and 180 small magnitude Ritz-values. For the initialization of the adaptive shifts in both models considering Bˇ = Rnξ ×100 , 10 optimal shift parameters are selected following the approach discussed in Subsection 4.4. The next of shifts are updated after steps every 10 iterations. It is already mentioned that this paper updates the LR-ADI iteration presented in [26, Algorithm 3] for solving the projected Lyapunov equations arising from second-order index 3 systems. The proposed algorithm i.e., Algorithm 3 in this paper is compared with [26, Algorithm 3] in Fig. 3 in terms of computational time. The figure shows that [26, Algorithm 3] is more expensive than the proposed algorithm in this paper. The computed Gramian factors are used in Algorithm 2 to obtain the reduced order models for both the systems considering the truncation tolerances  = 10−5 and 10−3 for DSMS and TCOM respectively. Several reduced models of varying

JID:APNUM AID:3724 /FLA

[m3G; v1.261; Prn:18/12/2019; 9:15] P.12 (1-14)

M. Monir Uddin / Applied Numerical Mathematics ••• (••••) •••–•••

12

Fig. 3. Comparison of Algorithm 3 and [26, Algorithm 3].

Table 3 Dimensions of original and various reduced systems with different balancing levels. Models

DSMS TCOM

nξ + nη



2001 11001

10−5

10−3

k pp

pv

vp

vv

66 44

59 44

66 38

55 41

Fig. 4. Comparison of different dimensional second-order reduced and original models in the frequency domain for the DSMS.

sizes are computed with different balancing levels and summarized in Table 3. Figs. 4 and 5 compare the reduced models with the original models of the DSMS and TCOM respectively. The frequency responses of the full and the different reduced models and their absolute and relative errors for the DSMS are shown in Fig. 4 over the frequency interval ω = [10−2 , 100 ]. Fig. 4a shows that the frequency responses of all reduced models are matching with the original model with good accuracy. The absolute and relative deviations between full and different reduced-order models are shown in Figs. 4b and 4c, respectively. On the other hand, Fig. 5 depicts the frequency responses, absolute and relative errors of the full and the reduced models of the TCOM application over the frequency interval ω = [10−3 , 100 ]. This figure also shows (in Fig. 5a) that the frequency responses of the reduced models with various balancing levels are matching correctly with the original model. Figs. 5b and 5c show a good approximation between the original and reduced models using the absolute and relative errors, respectively. Note that for both model the absolute and relative errors are respectively measured by

ˆ j ω)H∞ ( j ω) − (

and

ˆ j ω)H∞ ( j ω) − ( ( j ω)H∞

where ω = [ωmin , ωmax ] is a frequency interval. Although, in all balancing levels of DSMS the absolute errors are bounded by the prescribed truncation tolerance, the TCOM does not obey such error bound in all frequencies. This reflects that structure preserving model reduction does not have global error bound.

JID:APNUM AID:3724 /FLA

[m3G; v1.261; Prn:18/12/2019; 9:15] P.13 (1-14)

M. Monir Uddin / Applied Numerical Mathematics ••• (••••) •••–•••

13

Fig. 5. Comparison of different dimensional second-order reduced and original models in the frequency domain for the TCOM.

6. Conclusions This paper has discussed a balancing based technique to find a reduced second-order system from a large-scale sparse second-order index-3 system. In particular, the model considered here was obtained by linearizing the differential equation of motion with holonomic constrains around an equilibrium point in structural dynamics. The paper has shown that projecting the index-3 system onto the hidden manifold, then the algebraic constraints can be removed and the system can be converted into the standard form. However, model reduction can be applied without forming the projected system explicitly. The paper has also discussed an efficient way to solve the projected Lyapunov equations arising from the secondorder index-3 system using the LR-ADI iteration. Moreover the computation of optimal ADI shift parameters for a better convergence of the LR-ADI iteration has been discussed. The efficiency of the proposed techniques are shown by numerical experiments applying to two structured dynamical models. Appendix A. Supplementary material Supplementary material related to this article can be found online at https://doi.org/10.1016/j.apnum.2019.12.010. References [1] M.I. Ahmad, P. Benner, Interpolatory model reduction techniques for linear second-order descriptor systems, in: 2014 European Control Conference, ECC, IEEE, 2014, pp. 1075–1079. [2] A. Antoulas, Approximation of Large-Scale Dynamical Systems, SIAM, Philadelphia, PA, ISBN 978-0-89871-529-3, 2005. [3] U. Baur, Control-Oriented Model Reduction for Parabolic Systems, Ph.D. Thesis, Technische Universität Berlin, Germany, 2008. [4] P. Benner, J.-R. Li, T. Penzl, Numerical solution of large-scale Lyapunov equations, Riccati equations, and linear-quadratic optimal control problems, Numer. Linear Algebra Appl. 15 (9) (2008) 755–777, https://doi.org/10.1002/nla.622. [5] P. Benner, J. Saak, M.M. Uddin, Second order to second order balancing for index-1 vibrational systems, in: 7th International Conference on Electrical & Computer Engineering, ICECE, 2012, IEEE, 2012, pp. 933–936. [6] P. Benner, P. Kürschner, J. Saak, An improved numerical method for balanced truncation for symmetric second order systems, Math. Comput. Model. Dyn. Syst. 19 (6) (2013) 593–615, https://doi.org/10.1080/13873954.2013.794363. [7] P. Benner, P. Kürschner, J. Saak, Efficient handling of complex shift parameters in the low-rank Cholesky factor ADI method, Numer. Algorithms 62 (2) (2013) 225–251, https://doi.org/10.1007/s11075-012-9569-7. [8] P. Benner, P. Kürschner, J. Saak, Self-generating and efficient shift parameters in ADI methods for large Lyapunov and Sylvester equations, Electron. Trans. Numer. Anal. 43 (2014) 142–162, https://doi.org/10.17617/2.2071065. [9] P. Benner, J. Saak, M.M. Uddin, Balancing based model reduction for structured index-2 unstable descriptor systems with application to flow control, Numer. Algebra Control Optim. 6 (1) (2016) 1–20, https://doi.org/10.3934/naco.2016.6.1. [10] Y. Chahlaoui, D. Lemonnier, A. Vandendorpe, P. Van Dooren, Second-order balanced truncation, Linear Algebra Appl. 445 (2–3) (2006) 373–384. [11] E. Eich-Soellner, C. Führer, Numerical Methods in Multibody Dynamics, European Consortium for Mathematics in Industry, B. G. Teubner GmbH, Stuttgart, ISBN 3-519-02601-5, 1998. [12] S. Gugercin, A.C. Antoulas, C.A. Beattie, H2 model reduction for large-scale dynamical systems, SIAM J. Matrix Anal. Appl. 30 (2) (2008) 609–638, https://doi.org/10.1137/060666123. [13] S. Gugercin, T. Stykel, S. Wyatt, Model reduction of descriptor systems by interpolatory projection methods, SIAM J. Sci. Comput. 35 (5) (2013) B1010–B1033, https://doi.org/10.1137/130906635. [14] M. Heinkenschloss, D.C. Sorensen, K. Sun, Balanced truncation model reduction for a class of descriptor systems with application to the Oseen equations, SIAM J. Sci. Comput. 30 (2) (2008) 1038–1063, https://doi.org/10.1137/070681910. [15] N. Liu, A.E. Jeffers, A geometrically exact isogeometric Kirchhoff plate: feature-preserving automatic meshing and C 1 rational triangular Bézier spline discretizations, Int. J. Numer. Methods Eng. 115 (3) (2018) 395–409, https://doi.org/10.1002/nme.5809.

JID:APNUM AID:3724 /FLA

14

[m3G; v1.261; Prn:18/12/2019; 9:15] P.14 (1-14)

M. Monir Uddin / Applied Numerical Mathematics ••• (••••) •••–•••

[16] V. Mehrmann, T. Stykel, Balanced truncation model reduction for large-scale systems in descriptor form, in: Dimension Reduction of Large-Scale Systems, Springer-Verlag, ISBN 978-3-540-24545-2, 2005, pp. 357–361, Ch. 20. [17] B.C. Moore, Principal component analysis in linear systems: controllability, observability, and model reduction, IEEE Trans. Autom. Control AC–26 (1) (1981) 17–32, https://doi.org/10.1109/TAC.1981.1102568. [18] T. Penzl, A cyclic low rank Smith method for large sparse Lyapunov equations, SIAM J. Sci. Comput. 21 (4) (2000) 1401–1418, https://doi.org/10.1137/ S1064827598347666. [19] T. Reis, T. Stykel, Balanced truncation model reduction of second-order systems, Math. Comput. Model. Dyn. Syst. 14 (5) (2008) 391–406, https:// doi.org/10.1080/13873950701844170. [20] R. Riaza, Differential-Algebraic Systems. Analytical Aspects and Circuit Applications, World Scientific Publishing Co. Pte. Ltd., Singapore, ISBN 978-981-279-180-1, 2008. [21] V. Simoncini, A new iterative method for solving large-scale Lyapunov matrix equations, SIAM J. Sci. Comput. 29 (3) (2007) 1268–1288, https:// doi.org/10.1137/06066120X. [22] T. Stykel, Gramian-based model reduction for descriptor systems, Math. Control Signals Syst. 16 (4) (2004) 297–319, https://doi.org/10.1007/s00498004-0141-4. [23] F. Tisseur, K. Meerbergen, The quadratic eigenvalue problem, SIAM Rev. 43 (2) (2001) 235–286, https://doi.org/10.1137/S0036144500381988. ´ Bounds on the trace of a solution to the Lyapunov equation with a general stable matrix, Syst. Control Lett. 56 (7–8) (2007) [24] N. Truhar, K. Veselic, 493–503, https://doi.org/10.1016/j.sysconle.2007.02.003. [25] M.M. Uddin, Computational methods for model reduction of large-scale sparse structured descriptor systems, Ph.D. Thesis, Otto-von-GuerickeUniversität, Magdeburg, Germany, 2015. [26] M.M. Uddin, Gramian-based model-order reduction of constrained structural dynamic systems, IET Control Theory Appl. 12 (17) (2018) 2337–2346, https://doi.org/10.1049/iet-cta.2018.5580.