Abstract

Some new and simple Gramian-based model order reduction algorithms are presented on second-order linear dynamical systems, namely, SVD methods. Compared to existing Gramian-based algorithms, that is, balanced truncation methods, they are competitive and more favorable for large-scale systems. Numerical examples show the validity of the algorithms. Error bounds on error systems are discussed. Some observations are given on structures of Gramians of second order linear systems.

1. Introduction

Models of linear dynamical systems are often given in second-order form𝑀̈𝑞(𝑡)+𝐺̇𝑞(𝑡)+𝐾𝑞(𝑡)=𝐵0𝑦𝑢(𝑡)(𝑡)=𝐶0𝑞(𝑡)+𝐷0̇𝑞(𝑡),(1)where 𝑀,𝐺,𝐾𝑛×𝑛, 𝐵0𝑛×𝑚, 𝐶0, and 𝐷0𝑝×𝑛 are given matrices, 𝑢(𝑡)𝑚 is the given vector of inputs, 𝑦(𝑡)𝑝 is the unknown vector of outputs, 𝑞(𝑡)𝑛 is the unknown vector of internal variables, 𝑛 is the dimension of the system, and 𝑀 is assumed to be invertible.

Models of mechanical systems in particular are usually of second-order form (1). For such a system, 𝑀=𝑀𝑇, 𝐺 and 𝐾=𝐾𝑇 are respectively the mass, damping, and stiffness matrices, 𝑢(𝑡)𝑚 is the input, 𝐵0𝑢(𝑡)=𝑓(𝑡)𝑛 is the vector of external forces, and 𝑞(𝑡)𝑛 is the vector of internal variables.

Many applications lead to very large models where 𝑛 is very large, while 𝑚,𝑝𝑛. Due to limitations on time and computer storage, it is often difficult or impossible to directly simulate or control these large-scale systems. Therefore, it is often desirable to reduce the system to a smaller-order model:𝐵𝑀̈𝑞(𝑡)+𝐺̇𝑞(𝑡)+𝐾𝑞(𝑡)=0𝑢𝐶(𝑡)̂𝑦(𝑡)=0𝑞𝐷(𝑡)+0̇𝑞(𝑡),(2) such that (1)𝑦̂𝑦/𝑢<tol, that is, outputs corresponding to the same inputs are close. (2)The reduced system preserves important system properties such as stability and passivity. (3)The procedure is computationally efficient.

By letting 𝑥=𝑞̇𝑞, second-order system (1) is then transformed to a corresponding 2𝑛-dimensional first-order system Σ̇𝑥(𝑡)=𝐴𝑥(𝑡)+𝐵𝑢(𝑡)𝑦(𝑡)=𝐶𝑥(𝑡)+𝐷𝑢(𝑡),(3)

where 𝐴=0𝐼𝑀1𝐾𝑀1𝐺=0𝐼𝐾𝑀𝐺𝑀0𝑀𝐵=1𝐵0=0𝐵𝑀𝐶𝐶=0𝐷0𝐷=0.(4)

Since algorithms and theory are well developed for first-order model reduction, it is natural to use these techniques to develop algorithms for second-order systems.

Closely related to first-order system are the two Lyapunov equations𝐴𝒫+𝒫𝐴+𝐵𝐵𝐴=0,𝒬+𝒬𝐴+𝐶𝐶=0.(5)Under assumption that first-order system is stable, it is well known that the reachability and observability Gramians 𝒫 and 𝒬 are the unique solutions to Lyapunov equation (5), where both 𝒫 and 𝒬 are symmetric positive definite. This Gramian 𝒫 has also the following variational interpretation proposed in [1]. Let 𝐽(𝑢,𝑡1,𝑡2)=𝑡2𝑡1𝑢(𝑡)𝑢(𝑡)𝑑𝑡=𝑡2𝑡1𝑢(𝑡)2𝑑𝑡. The optimal value to optimization problem min𝑢𝐽(𝑢,,0)𝑠.𝑡.̇𝑥=𝐴𝑥+𝐵𝑢,𝑥(0)=𝑥0(6)

is 𝑥0𝒫1𝑥0, that is, the minimal energy required to steer the state of the system from 𝑡= to state 𝑥0 at time 𝑡=0 is 𝑥0𝒫1𝑥0. Similarly, the optimal value to its dual system is 𝑥0𝒬1𝑥0.

Partition 𝒫 and 𝒬 into four equal blocks: 𝒫𝒫=11𝒫12𝒫21𝒫22𝒬,𝒬=11𝒬12𝒬21𝒬22.(7)

In [24], it is shown that the optimal value to optimization problemmiṅ𝑞0𝑛min𝑢𝐽(𝑢,,0)𝑠.𝑡.𝑀̈𝑞(𝑡)+𝐺̇𝑞(𝑡)+𝐾𝑞(𝑡)=𝐵𝑢(𝑡),𝑞(0)=𝑞0(8) is 𝑞0𝒫111𝑞0, that is, the minimal energy required to reach the given position 𝑞0 over all past inputs and initial velocities. And the optimal value to optimization problemmin𝑞0𝑛min𝑢𝐽(𝑢,,0)s.t.𝑀̈𝑞(𝑡)+𝐺̇𝑞(𝑡)+𝐾𝑞(𝑡)=𝐵𝑢(𝑡),̇𝑞(0)=̇𝑞0(9) is ̇𝑞0𝒫122̇𝑞0, that is, the minimal energy required to reach the given velocity ̇𝑞0 over all past inputs and initial positions. In [3, 4], 𝒫11 and 𝒫22 are defined to be position Gramian and velocity Gramian, respectively. By duality, 𝒬11 and 𝒬22 are position Gramian and velocity Gramian, respectively corresponding to 𝒬.

It is well known that 𝒫=0𝑒𝐴𝜏𝐵𝐵𝑇𝑒𝐴𝑇𝜏𝑑𝜏,𝒬=0𝑒𝐴𝑇𝜏𝐶𝑇𝐶𝑒𝐴𝜏𝑑𝜏.(10)

Suppose 𝑥(𝑡) is the solution of first-order system for impulse input 𝑢(𝑡)=𝛿(𝑡). Then𝒫=0𝑥(𝑡)𝑥(𝑡)𝑑𝑡.(11) For second-order system and its corresponding first-order system, 𝒫=0𝑥(𝑡)𝑥=(𝑡)𝑑𝑡0𝑞𝑞(𝑡)̇𝑞(𝑡)(𝑡)̇𝑞=(𝑡)𝑑𝑡0𝑞(𝑡)𝑞(𝑡)𝑞(𝑡)̇𝑞(𝑡)̇𝑞(𝑡)𝑞(𝑡)̇𝑞(𝑡)̇𝑞(𝑡)𝑑𝑡.(12)

Therefore,𝒫11=0𝑞(𝑡)𝑞𝒫(𝑡)𝑑𝑡,22=0̇𝑞(𝑡)̇𝑞(𝑡)𝑑𝑡,(13)where 𝒫11 and 𝒫22 are position and velocity Gramians, respectively. Similar results can be applied to 𝒬, 𝒬11, and 𝒬22 through dual systems.

For Gramian based structure preserving mode order reduction of second-order systems, in [3, 5], some balanced truncation methods were derived. In this paper, we present two new algorithms based on SVD rather than balanced truncation. The algorithms are derived from second-order Gramians and system 2-norm. Compared to existing techniques, they are simple and more suitable for large-scale settings. In Section 4, we apply these methods to three numerical examples, the computational results show the validity of our algorithms. Error bounds and structures of Gramians on second-order systems are also discussed.

2. Balanced Truncation Method

For first-order systemΣ̇𝑥(𝑡)=𝐴𝑥(𝑡)+𝐵𝑢(𝑡),𝑦(𝑡)=𝐶𝑥(𝑡).(14) balanced truncation method is to transform the system to another coordinate system such that both 𝒫 and 𝒬 are equal and diagonal: 𝜎𝒫=𝒬=diag1,𝜎2,,𝜎𝑁,(15)

where 𝜎1𝜎2𝜎𝑁. It then truncates the model by keeping the states corresponding to 𝑘 largest eigenvalues of the product 𝒫𝒬, that is, 𝑘 largest Hankel singular values of the system.

The main problem is to find balancing projection matrices. A numerically stable way to get the balancing truncated system is given in [6]: suppose the Cholesky factorization of 𝒫 and 𝒬 are is 𝒫=𝑈𝐶𝑈𝑇𝐶,𝒬=𝐿𝑂𝐿𝑇𝑂,(16)

where 𝑈𝐶 is upper triangular and 𝐿𝑂 is lower triangular matrices. Take SVD of 𝑈𝑇𝐶𝐿𝑂, 𝑈𝑇𝐶𝐿𝑂=𝑈Σ𝑉𝑇.(17)

Let Σ𝑘 be the first 𝑘 by 𝑘 principle submatrix of Σ, 𝑈𝑘, and 𝑉𝑘 consist of the first 𝑘 columns of 𝑈 and 𝑉, respectively. Then the balanced projection matrices are 𝑊=𝐿𝑂𝑉𝑘Σ𝑘1/2,𝑉=𝑈𝐶𝑈𝑘Σ𝑘1/2.(18)

The reduced system is obtained in: 𝐴=𝑊𝑇𝐴𝑉,𝐵=𝑊𝑇𝐵,𝐶=𝐶𝑉.(19)

Balanced truncation method has a global error bound [1]: Σerr=ΣΣred𝜎2𝑘+1++𝜎𝑛.(20)

For second-order system 𝑀̈𝑞(𝑡)+𝐺̇𝑞(𝑡)+𝐾𝑞(𝑡)=𝐵0𝑦𝑢(𝑡)(𝑡)=𝐶0𝑞(𝑡)+𝐷0̇𝑞(𝑡),(21)

there are two balanced truncation methods. One was given in [5] which is to balance both 𝒫11 and 𝒬11 to form the projection matrices 𝑊1𝑟 and 𝑉1𝑟 and then get the reduced system by keeping 𝑟 largest eigenvalues of 𝒫11𝒬11. It is equivalent to performing projections to its corresponding first-order system as following:𝐴=𝑊𝑇=𝑊𝐴𝑉𝑇1𝑟𝑊𝑇1𝑟0𝐼𝑀1𝐾𝑀1𝐺𝑉1𝑟𝑉1𝑟=0𝐼𝑊11𝑟𝑀1𝐾𝑉1𝑟𝑊11𝑟𝑀1𝐺𝑉1𝑟,𝑊𝐵=𝑊𝐵=𝑇1𝑟𝑊𝑇1𝑟0𝑀1𝐵0=0𝑊𝑇1𝑟𝑀1𝐵0𝐶𝐶=𝐶𝑉=0𝐷0𝑉1𝑟𝑉1𝑟=𝐶0𝑉1𝑟𝐷0𝑉1𝑟.(22)The reduction rules are then𝑀=𝐼𝑟,𝐾=𝑊𝑇1𝑟𝑀1𝐾𝑉1𝑟,𝐺=𝑊𝑇1𝑟𝑀1𝐺𝑉1𝑟,𝐵0=𝑊𝑇1𝑟𝑀1𝐵0,𝐶0=𝐶0𝑉1𝑟,𝐷0=𝐷0𝑉1𝑟.(23)Similar results can be applied to 𝒫22 and 𝒬22. This gives the following algorithm.

Algorithm 1 (Balanced truncation method based on 𝒫11 and 𝒬11 (or 𝒫22 and 𝒬22, resp.) [5]). (1) Compute 𝒫11 and 𝒬11.(2) Compute the balanced truncation matrices 𝑊1𝑟,𝑉1𝑟𝑛×𝑟 such that 𝑊𝑇1𝑟𝑉1𝑟=𝐼𝑟 and 𝑊𝑇1𝑟𝒫11𝑊1𝑟=𝑉𝑇1𝑟𝒬11𝑉1𝑟=diag(𝜎1,𝜎2,,𝜎𝑟), where 𝜎21,𝜎22,,𝜎2𝑟 are the 𝑟 largest eigenvalues of 𝒫11𝒬11. (3) Perform projection to (𝑀,𝐺,𝐾,𝐵0,𝐶0,𝐷0) and get (𝐵𝑀,𝐺,𝐾,0,𝐶0,𝐷0): 𝑀=𝐼𝑟,𝐾=𝑊𝑇1𝑟𝑀1𝐾𝑉1𝑟,𝐺=𝑊𝑇1𝑟𝑀1𝐺𝑉1𝑟,𝐵0=𝑊𝑇1𝑟𝑀1𝐵0,𝐶0=𝐶0𝑉1𝑟,𝐷0=𝐷0𝑉1𝑟.(24)

In [5], it also gave two other similar methods. One is to balance 𝒫11 and 𝒬22 to get the projection matrices, the other is to balance 𝒫22 and 𝒬11 to get the projection matrices.

Another second-order balanced truncation method is given in [3] which is to balance both 𝒫11 and 𝒬11 to get the projection matrices 𝑊1𝑟 and 𝑉1𝑟 by keeping 𝑟 largest eigenvalues of 𝒫11𝒬11, balance both 𝒫22 and 𝒬22 to get the projection matrices 𝑊2𝑟 and 𝑉2𝑟, then use 𝑊=𝑊1𝑟𝑊2𝑟 and 𝑉=𝑉1𝑟𝑉2𝑟 as projection matrices for the corresponding first-order system: 𝐴=𝑊𝑇=𝑊𝐴𝑉𝑇1𝑟𝑊𝑇2𝑟0𝐼𝑀1𝐾𝑀1𝐺𝑉1𝑟𝑉2𝑟=0𝑊𝑇1𝑟𝑉2𝑟𝑊𝑇2𝑟𝑀1𝐾𝑉1𝑟𝑊𝑇2𝑟𝑀1𝐺𝑉2𝑟,𝐵=𝑊𝑇𝑊𝐵=𝑇1𝑟𝑊𝑇2𝑟0𝑀1𝐵0=0𝑊𝑇2𝑟𝑀1𝐵0𝐶𝐶=𝐶𝑉=0𝐷0𝑉1𝑟𝑉2𝑟=𝐶0𝑉1𝑟𝐷0𝑉2𝑟.(25)

In order to let reduced system have the companion form of second-order system, it then takes the following transformation by letting 𝐻1=𝑊𝑇1𝑟𝑉2𝑟,𝐼𝐻𝐴=1𝐴𝐼𝐻=𝐼𝐻10𝑊𝑇1𝑟𝑉2𝑟𝑊𝑇2𝑟𝑀1𝐾𝑉1𝑟𝑊𝑇2𝑟𝑀1𝐺𝑉2𝑟𝐼𝐻=0𝐼𝐻1𝑊𝑇2𝑟𝑀1𝐾𝑉1𝑟𝐻1𝑊𝑇2𝑟𝑀1𝐺𝑉2𝑟𝐻,𝐼𝐻𝐵=1𝐼𝐻𝐵=10𝑊𝑇2𝑟𝑀1𝐵0=0𝐻1𝑊𝑇2𝑟𝑀1𝐵0,𝐶𝐼𝐻=𝐶𝐶=0𝑉1𝑟𝐷0𝑉2𝑟𝐼𝐻=𝐶0𝑉1𝑟𝐷0𝑉2𝑟𝐻.(26)The corresponding algorithm is as follows:

Algorithm 2 (Balanced truncation method based on 𝒫11𝒫22 and 𝒬11𝒬22 [3]). (1) Compute 𝒫11, 𝒫22, 𝒬11, and 𝒬22.(2) Compute the balanced truncation matrices for 𝒫11 and 𝒬11 to get 𝑊1𝑟,𝑉1𝑟𝑛×𝑟, and balanced transformation matrices for 𝒫22 and 𝒬22 to get 𝑊2𝑟,𝑉2𝑟𝑛×𝑟.(3) Perform projection to get reduced system (𝐵𝑀,𝐺,𝐾,0,𝐶0,𝐷0):𝑀=𝐼𝑟,𝐾=𝐻1𝑊𝑇2𝑟𝑀1𝐾𝑉1𝑟,𝐺=𝐻1𝑊𝑇2𝑟𝑀1𝐺𝑉2𝑟𝐵𝐻,0=𝐻1𝑊𝑇2𝑟𝑀1𝐵0,𝐶0=𝐶0𝑉1𝑟,𝐷0=𝐷0𝑉2𝑟𝐻,(27) where 𝐻1=𝑊𝑇1𝑟𝑉2𝑟.

3. Some New Algorithms Based on SVD

In this section, we propose some new algorithms based on SVD method directly.

First, 2 norm is an important quantity for bounding linear systems.

Lemma 3 (see [7]). Given a first-order system 𝐴Σ=𝐵𝐶,(28)suppose 𝒫 and 𝒬 are reachability and observability Gramians. Then, 2 norm of the system can be expressed as Σ22=trace𝐶𝒫𝐶𝑇𝐵=trace𝑇.𝒬𝐵(29)

2 norm for second-order system (1) with 𝐷0=0 can then be stated as in the following proposition.

Proposition 4. Given a second-order system (1) with 𝐷0=0, suppose 𝐴Σ=𝐵𝐶(30)is the transformed first-order system, where 𝐴=0𝐼𝑀1𝐾𝑀1𝐺,0𝑀𝐵=1𝐵0,𝐶𝐶=00.(31)Assume reachability and observability Gramians 𝒫 and 𝒬 are divided into four equal blocks, and 𝒫=𝒫11𝒫12𝒫𝑇12𝒫22, 𝒬=𝒬11𝒬12𝒬𝑇12𝒬22. Then the following holds: Σ22𝐵=trace𝑇0𝑀𝑇𝒬22𝑀1𝐵0𝐶=trace0𝒫11𝐶𝑇0.(32)

Proof. From Lemma 3, Σ22𝐵=trace𝑇0𝑀𝒬𝐵=trace1𝐵0𝑇𝒬11𝒬12𝒬𝑇12𝒬220𝑀1𝐵0𝐵=trace𝑇0𝑀𝑇𝒬22𝑀1𝐵0,(33)Σ22𝐶=trace00𝒫11𝒫12𝒫𝑇12𝒫22𝐶𝑇00𝐶=trace0𝒫11𝐶𝑇0.(34)

From (8), we know the positions which are difficult to reach are spanned by the eigenvectors corresponding to small eigenvalues of 𝒫11. We would like to eliminate those positions in reduced systems. Proposition 4 shows that 𝒫11 can independently decide the system 2 norm. Therefore, we may get the reduced system by keeping the positions corresponding to 𝑟 largest eigenvalues of 𝒫11. This motivates the following procedure for model order reduction. Suppose SVD of 𝒫11 is 𝒫11=𝑊1𝑆𝑉𝑇1, where 𝑆=diag(𝜎1,𝜎2,,𝜎𝑛), and 𝜎1𝜎2𝜎𝑛>0. So 𝑊1,𝑉1 are orthogonal matrices and 𝑊1=𝑉1 since 𝒫11 is symmetric positive definite. Let the transformation matrices for the corresponding first-order system be𝑊𝑊=1𝑊1,𝑉=𝑊,(35)

and the projection matrices be𝑊𝑘=𝑊1𝑟𝑊1𝑟,𝑉𝑘=𝑊𝑘,(36)

where 𝑊1𝑟 consists of the first 𝑟 columns of 𝑊1. For simplicity, we denote 𝑊𝑘 and 𝑉𝑘 by 𝑊 and 𝑉, respectively. Then use 𝑊 and 𝑉 as projection matrices for the corresponding first-order system, that is, perform projections as in (22) and (23). Same results can be applied to 𝒬22. From the dual system, symmetrically 𝒬11 and 𝒫22 are also crucial in weighting the system 2 norms independently. This results in the following algorithms.

Algorithm 5 (Second-order model reduction—SVD on 𝒫11 (or, 𝒫22, 𝒬11, and 𝒬22, resp.)). (1) Compute 𝒫11 (or 𝒫22, 𝒬11, and 𝒬22, resp.), and denote it as 𝑃 for simplicity.(2) Take SVD on 𝑃 and get 𝑃=𝑊1𝑆1𝑊𝑇1. Form the orthogonal projection matrices 𝑊1𝑟 which consists of the first 𝑟 columns of 𝑊1.(3) Perform projection to get the reduced system (𝐵𝑀,𝐺,𝐾,0,𝐶0,𝐷0): 𝑀=𝐼𝑟,𝐾=𝑊𝑇1𝑟𝑀1𝐾𝑊1𝑟,𝐺=𝑊𝑇1𝑟𝑀1𝐺𝑊1𝑟,𝐵0=𝑊𝑇1𝑟𝑀1𝐵0,𝐶0=𝐶0𝑊1𝑟,𝐷0=𝐷0𝑊1𝑟.(37)

Besides eliminating the positions which are difficult to reach in reduced system, it is also desirable to delete the velocities that are difficult to reach. From (9), these velocities are spanned by the eigenvectors of 𝒫22 corresponding to small eigenvalues. This motivates the following procedure for model order reduction. Suppose SVD of 𝒫11 is 𝒫11=𝑊1𝑆1𝑉𝑇1, where 𝑆1=diag(𝜎1,𝜎2,,𝜎𝑛), 𝜎1𝜎2𝜎𝑛>0. And SVD of 𝒫22 is 𝒫22=𝑊2𝑆2𝑉𝑇2, where 𝑆2=diag(𝛿1,𝛿2,,𝛿𝑛), 𝛿1𝛿2𝛿𝑛>0. So 𝑊1, 𝑉1, 𝑊2, and 𝑉2 are all orthogonal matrices, and 𝑊1=𝑉1, 𝑊2=𝑉2. Let the transformation matrices for the corresponding first-order system be𝑊𝑊=1𝑊2,𝑉=𝑊,(38)

and the projection matrices be𝑊𝑘=𝑊1𝑟𝑊2𝑟,𝑉𝑘=𝑊𝑘,(39)

where 𝑊1𝑟 and 𝑊2𝑟 consist of the first 𝑟 columns of 𝑊1 and 𝑊2, respectively. By using similar idea to Algorithm 2, in order to let the reduced system 𝐴 have the companion form of second-order model, we adopt a matrix 𝐻 such that 𝐻1=𝑊𝑇1𝑟𝑉2𝑟, and then perform the projections as in (26) and (27). Similar results can be applied to 𝒬11 and 𝒬22. This results in the following algorithm.

Algorithm 6 (Second-order model reduction—SVD on 𝒫11 and 𝒫22 (or 𝒬11 and 𝒬22, resp.)). (1) Compute 𝒫11 and 𝒫22.(2) Take SVD on 𝒫11 and 𝒫22 to get 𝒫11=𝑊1𝑆1𝑊𝑇1 and 𝒫22=𝑊2𝑆2𝑊𝑇2, respectively.(3) Perform projection to get reduced system (𝐵𝑀,𝐺,𝐾,0,𝐶0,𝐷0): 𝑀=𝐼𝑟,𝐾=𝐻1𝑊𝑇2𝑟𝑀1𝐾𝑊1𝑟,𝐺=𝐻1𝑊𝑇2𝑟𝑀1𝐺𝑊2𝑟𝐵𝐻,0=𝐻1𝑊𝑇2𝑟𝑀1𝐵0,𝐶0=𝐶0𝑊1𝑟,𝐷0=𝐷0𝑊2𝑟𝐻,(40) where 𝐻1=𝑊𝑇1𝑟𝑊2𝑟.

Now, consider Algorithm 5 and SISO second-order systems. Note that by taking SVD,  𝒫11=𝑊1𝑆1𝑊𝑇1. Denote 𝑆11=diag(𝜎1,𝜎2,,𝜎𝑟) and 𝑆12=diag(𝜎𝑟+1,𝜎𝑟+2,,𝜎𝑛). So 𝑆1=diag(𝑆11,𝑆12). Suppose 𝑧=𝑊𝑇1𝐶𝑇0, that is,  𝑧1𝑧2𝑧𝑛=𝑊𝑇1𝑐1𝑐2𝑐𝑛.(41)Then from Proposition 4, Σ22=𝐶0𝒫11𝐶𝑇0=𝐶0𝑊1𝑆𝑊𝑇1𝐶𝑇0=𝑧𝑇𝑆𝑧=𝜎1𝑧21+𝜎2𝑧22++𝜎𝑛𝑧2𝑛.(42)

In reduced system, 𝐶0=𝐶0𝑊1𝑟, 𝒫11=𝑆1. Therefore, Σred22=𝐶0𝒫11𝐶𝑇0=𝐶0𝑊1𝑟𝑆1𝑊𝑇1𝑟𝐶𝑇0=𝜎1𝑧21+𝜎2𝑧22++𝜎𝑟𝑧2𝑟.(43)

In many cases, the eigenvalues of 𝒫11 decay rapidly. Therefore, the reduced system produced by Algorithm 5 takes the main part of the system 2 norm. Actually we also have, Σ22Σred22=𝜎𝑟+1𝑧2𝑟+1+𝜎𝑟+2𝑧2𝑟+2++𝜎𝑛𝑧2𝑛𝜎𝑟+1𝑧22=𝜎𝑟+1𝐶022.(44)

This is not the 2-norm of error system ΣΣred2, but it gives some information to support Algorithm 5. This result for SISO second-order systems can easily extends to MIMO systems. Similar results can be applied to methods in Algorithm 6.

Now, consider complexity of the algorithms. In our SVD method on 𝒫11 in Algorithm 5, after computing 𝒫11 which takes 𝑂(𝑛3) time, it then uses one-time 𝑂(𝑛3) operation which is SVD on 𝒫11 to get the projection matrices. In our Algorithm 6, after computing 𝒫11 and 𝒫22 each of which takes 𝑂(𝑛3) time, it then uses twice 𝑂(𝑛3) operations which are SVD on 𝒫11 and SVD on 𝒫22 to get the projection matrices. In balanced truncation method in Algorithm 1, after computing 𝒫11 and 𝒬11, it also has one-time SVD. Besides this, it uses two Cholesky factorizations 𝒫11=𝑈𝐶𝑈𝑇𝐶 and 𝒬11=𝐿𝑂𝐿𝑇𝑂 which all take time 𝑂(𝑛3) in order to get the projection matrices. Note that Algorithm 2 has more flops than Algorithm 1. Therefore, our SVD algorithms are simpler in time cost than existing balanced truncation methods. There are some other techniques in getting the balanced truncation systems, for example, balance-free truncation [8]. But it is not hard to see that they are also more complex than our SVD methods. This makes our SVD methods more suitable to solve large-scale problems.

4. Computational Results

In this section, we apply the model reduction methods to three numerical examples: the building model, the clamped beam model, and international space station model, see [9] for detailed descriptions. In Table 1, four model reduction methods are compared: SVD method on 𝒫11 in Algorithm 5 (“SVD1”), SVD method on 𝒫11 and 𝒫22 in Algorithm 6 (“SVD12”), balanced truncation on 𝒫11 and 𝒬11 in Algorithm 1 (“BT1”), and balanced truncation on 𝒫11 and 𝒬11, 𝒫22 and 𝒬22 in Algorithm 5 (“BT12”). The comparison is based on the relative error of 2 norm, where 𝜀22Σ=Σ2Σ2,(45)

𝑛 is the size of 𝑀, 𝑟 is the size of 𝑀 in reduced model, and 𝑚 is the size of input vector 𝑢(𝑡), 𝑝 is the size of output vector 𝑦(𝑡). Figures 1 and 2 show the amplitude of frequency response of original and reduced systems for clamped beam model for 𝑟=20 and 𝑟=29, respectively.

5. Error Bounds and Discussions

Error bounds for first-order model reduction do not apply for second-order systems. Consider second-order system (1) with 𝐷0=0Σ𝑀̈𝑞(𝑡)+𝐺̇𝑞(𝑡)+𝐾𝑞(𝑡)=𝐵0𝑢(𝑡)𝑦(𝑡)=𝐶0𝑞(𝑡).(46)

Suppose the reduced system Σ is obtained by keeping 𝑟 largest eigenvalues of the orthogonal eigenspace of 𝒫11. Sorensen and Antoulas in [10] provided a priori error bound showing that the 2 norm of error system is bounded by a constant times the summation of neglected singular values of 𝒫11.

It considers structured systems in [10]. Let 𝑄(𝑠) and 𝑃(𝑠) be polynomial matrices in 𝑠: 𝑄(𝑠)=𝑟𝑗=1𝑄𝑗𝑠𝑗,𝑄𝑗𝑛×𝑛,𝑃(𝑠)=𝑟1𝑗=1𝑃𝑗𝑠𝑗,s𝑃𝑗𝑛×𝑚,(47)

where 𝑄 is invertible, and 𝑄1𝑃 is a strictly proper rational matrix. Denote by 𝑄(𝑑/𝑑𝑡), 𝑃(𝑑/𝑑𝑡) the differential operators 𝑄𝑑=𝑑𝑡𝑟𝑗=1𝑄𝑗𝑑𝑗𝑑𝑡𝑗𝑑,𝑃=𝑑𝑡𝑟1𝑗=1𝑃𝑗𝑑𝑗𝑑𝑡𝑗.(48)

The structured systems are defined by the following equations:𝑄𝑑Σ𝑑𝑑𝑡𝑥(𝑡)=𝑃𝑑𝑡𝑢(𝑡)𝑦(𝑡)=𝐶𝑥(𝑡),(49)where 𝐶𝑝×𝑛. The Gramian is defined as the Gramian of 𝑥(𝑡) when the input is an impulse: 𝒫=0𝑥(𝑡)𝑥(𝑡)𝑑𝑡.(50)

Let 𝒫=𝑉𝑆𝑉(51)

be the eigensystem of 𝒫, where 𝑉 is orthogonal, the diagonal elements of 𝑆 are in decreasing order: 𝜎𝑆=1𝜎2𝜎𝑛,(52)

where 𝜎1𝜎2𝜎𝑛. Partition 𝑉=[𝑉1,𝑉2] where 𝑉1 consists of the first 𝑟 columns of 𝑉. The reduced model is derived from𝑄𝑗=𝑉1𝑄𝑗𝑉1,𝑃𝑗=𝑉1𝑃𝑗,𝐶=𝐶𝑉1,(53)and the reduced system is then𝑄𝑑Σ𝑃𝑑𝑑𝑡̂𝑥(𝑡)=𝑑𝑡𝑢(𝑡)̂𝑦(𝑡)=𝐶̂𝑥(𝑡).(54) Partition accordingly, 𝐶1,𝐶2𝑆=𝐶𝑉,𝑆=1𝑆2.(55)

Theorem 7 (see[10]). The reduced model Σ derived from the dominant eigenspace of the Gramian 𝒫 for Σ as described above satisfies ΣΣ22𝐶trace2𝑆2𝐶2𝑆+𝒦trace2,(56)where 𝒦 is a modest constant depending on Σ and Σ, and the diagonal elements of 𝑆2 are the neglected smallest eigenvalues of 𝒫.

Specially for second-order system (1) with 𝐷0=0:Σ𝑀̈𝑞(𝑡)+𝐺̇𝑞(𝑡)+𝐾𝑞(𝑡)=𝐵0𝑢(𝑡)𝑦(𝑡)=𝐶0𝑞(𝑡),(57)it is easy to see that it can be described by (49) with 𝑄(𝑠)=𝑀𝑠2+𝐺𝑠+𝐾 and 𝑃(𝑠)=𝐵0. So there is the following corollary.

Corollary 8. For second-order system (57), suppose 𝒫11=𝑆1𝑆2 is diagonal, where 𝑆1=diag(𝜎1,𝜎2,,𝜎𝑟), 𝑆2=diag(𝜎𝑟+1,𝜎𝑟+2,,𝜎𝑛), and 𝜎1𝜎2𝜎𝑛. Suppose the reduced system Σ is obtained by keeping 𝑟 largest eigenvalues of 𝒫11 as described in (53) and (54). Then Σ𝑒𝑟𝑟22=ΣΣ22𝑐0𝑆trace2=𝑐0𝜎𝑟+1+𝜎𝑟+2++𝜎𝑛,(58)where 𝑐0 is a modest constant depending on Σ and Σ.

Now, we explain the corollary and what the constant 𝑐0 is. Note 𝑄(𝑠)=𝑀𝑠2+𝐺𝑠+𝐾𝑛×𝑛 and 𝑄(𝑠)=𝑄11(𝑠)𝑄12𝑄(𝑠)21(𝑠)𝑄22(𝑠) with 𝑄11(𝑠)𝑟×𝑟. Let 𝑊(𝑠)=𝑄11(𝑠)1𝑄12(𝑠). Partition accordingly 𝐶0=[𝐶01,𝐶02]. Then, 𝑐0=sup𝜔𝐶01𝑊(𝑖𝜔)𝐶01𝑊(𝑖𝜔)2𝐶022+𝐶0222.(59)

If the reduced system has no poles on the imaginary axis, sup𝜔𝑊(𝑖𝜔)2 is finite, and so sup𝜔(𝐶01𝑊(𝑖𝜔))(𝐶01𝑊(𝑖𝜔)2𝐶02)2 is finite. Therefore, 𝑐0 is a constant depending on Σ and Σ. There is a similar result for 𝒬11 by duality.

Consider Algorithm 5, SVD method on 𝒫11. Note that 𝒫11=0𝑞(𝑡)𝑞(𝑡)𝑑𝑡,(60)

and 𝒫11=𝑉𝑆𝑉, [𝑉1,𝑉2]=𝑉. So the reduction rules in Algorithm 5 are exactly the same as (53). Similar result applies for 𝒬11. Therefore, the error systems produced by our Algorithm 5, SVD method on 𝒫11 or 𝒬11, all have good global error bounds. So far, there is no priori error bound for error systems based on 𝒫22 or 𝒬22.

Even though balanced truncation methods in Algorithms 1 and 2, and SVD method on both 𝒫11 and 𝒫22 (or 𝒬11 and 𝒬22) in Algorithm 6, all have good performance in those three numerical examples, so far no priori error bounds are provided for these methods.

For second-order system (1) with 𝐷0=0, the following theorem provides some observations on structures of gramians 𝒫 and 𝒬 in corresponding first-order system.

Theorem 9. Given second-order system (1) with 𝐷0=0, then in any coordinate system, there is no state space transformation which gives a block diagonal 𝒬.

Proof. From the observability Lyapunov equation 𝐴𝑇𝒬+𝒬𝐴+𝐶𝑇𝐶=0,(61) one can get 0𝐾𝑇𝑀𝐼𝐺𝑇𝑀𝒬11𝒬12𝒬𝑇12𝒬22+𝒬11𝒬12𝒬𝑇12𝒬220𝐼𝐾𝑀𝐺𝑀+𝐶𝑇00𝐶00=0.(62)Equating each block gives 𝐾𝑇𝑀𝒬𝑇12𝒬12𝐾𝑀+𝐶𝑇0𝐶0=0,𝐾𝑇𝑀𝒬22+𝒬11𝒬12𝐺𝑀𝒬=0,12𝐺𝑇𝑀𝒬22+𝒬𝑇12𝒬22𝐺𝑀=0.(63) From the first equation in (63), one can get 𝒬12 is not zero, otherwise 𝐶0=0, and then from (61) the solution for this Lyapunov equation is 𝒬=0𝑒𝐴𝑇𝜏𝐶𝑇𝐶𝑒𝐴𝜏𝑑𝜏=0. This contradicts with 𝒬 being symmetric positive definite. So 𝒬 is not block diagonal.

Algorithm 2 [3], that is, balanced truncation method on (𝒫11,𝒫22) and (𝒬11,𝒬22), is to balance both 𝒫11 and 𝒬11 to get the projection matrices 𝑊1𝑟 and 𝑉1𝑟 by keeping 𝑟 largest eigenvalues of 𝒫11𝒬11, balance both 𝒫22 and 𝒬22 to get the projection matrices 𝑊2𝑟 and 𝑉2𝑟, then use 𝑊=𝑊1𝑟𝑊2𝑟 and 𝑉=𝑉1𝑟𝑉2𝑟 as projection matrices for the corresponding first-order system. This is equivalent to assuming both 𝒫 and 𝒬 are block diagonal in corresponding first-order system, that is, 𝒫𝒫=11𝒫22𝒬,𝒬=11𝒬22,(64)

in order to keep 𝑟 largest eigenvalues of 𝒫11𝒬11 and 𝒫22𝒬22 in reduced corresponding first-order system. Similar results also apply for methods in Algorithm 6, SVD method on 𝒬11 and 𝒬22. From Theorem 9, 𝒬 can not be block diagonal, therefore the reduced system obtained from Algorithm 2 or Algorithm 6 may lose some information. But Algorithm 1 [5] (balanced truncation based on (𝒫11,𝒬11) (or (𝒫22,𝒬22) resp.)) and Algorithm 5 (SVD on 𝒫11 (or 𝒫22, 𝒬11, 𝒬22 resp.)), avoid this drawback in assuming 𝒬 being block diagonal, but they only work for either position Gramians or velocity Gramians, while Algorithms 2 and 6 have the advantage in taking into account both position Gramians and velocity Gramians.

6. Conclusions

Two Gramian based second-order model reduction methods are proposed in this paper based on SVD on diagonal blocks of Gramians. Even though they are the well known SVD methods, to our knowledge, they are the simplest ways in Gramian based second-order model reduction. Two existing balanced truncation methods [3, 5] were presented in 2006 and 2008, respectively. When compared to these two existing techniques, our SVD methods are competitive in numerical examples and more suitable for large-scale setting problems. By using the result given in [10], there is global error bound for one method, SVD on 𝒫11 or 𝒬11. Except this, so far no priori error bounds are provided for other methods on second-order model reduction, even through they work good in numerical examples. Interesting future work would be to provide global error bounds for these methods.

Acknowledgment

The paper is supported in part by NSFC research Project 61070230, Shandong Social Science Program of China 11CKJJ16, and Grant no. 20110205 by Jinan Science and Technology Bureau of China for Returned Oversea Scholars.