Table of Contents Author Guidelines Submit a Manuscript
Journal of Control Science and Engineering
Volumeย 2012, Article IDย 302498, 9 pages
http://dx.doi.org/10.1155/2012/302498
Research Article

Second-Order Model Reduction Based on Gramians

School of Mathematics and Quantitative Economics, Shandong University of Finance and Economics, Jinan, Shandong 250014, China

Received 22 August 2011; Accepted 2 January 2012

Academic Editor: S.ย Skogestad

Copyright ยฉ 2012 Cong Teng. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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 [2โ€“4], it is shown that the optimal value to optimization problemminฬ‡๐‘ž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๐ผโˆ’๐‘Š1โˆ’1๐‘Ÿ๐‘€โˆ’1๐พ๐‘‰1๐‘Ÿโˆ’๐‘Š1โˆ’1๐‘Ÿ๐‘€โˆ’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โŽžโŽŸโŽŸโŽ ๎๐ดโŽ›โŽœโŽœโŽ๐ผ๐ปโŽžโŽŸโŽŸโŽ =โŽ›โŽœโŽœโŽ๐ผ๐ปโˆ’1โŽžโŽŸโŽŸโŽ โŽ›โŽœโŽœโŽ0๐‘Š๐‘‡1๐‘Ÿ๐‘‰2๐‘Ÿโˆ’๐‘Š๐‘‡2๐‘Ÿ๐‘€โˆ’1๐พ๐‘‰1๐‘Ÿโˆ’๐‘Š๐‘‡2๐‘Ÿ๐‘€โˆ’1๐บ๐‘‰2๐‘ŸโŽžโŽŸโŽŸโŽ โŽ›โŽœโŽœโŽ๐ผ๐ปโŽžโŽŸโŽŸโŽ =โŽ›โŽœโŽœโŽ0๐ผโˆ’๐ปโˆ’1๐‘Š๐‘‡2๐‘Ÿ๐‘€โˆ’1๐พ๐‘‰1๐‘Ÿโˆ’๐ปโˆ’1๐‘Š๐‘‡2๐‘Ÿ๐‘€โˆ’1๐บ๐‘‰2๐‘Ÿ๐ปโŽžโŽŸโŽŸโŽ ,๎‚โŽ›โŽœโŽœโŽ๐ผ๐ป๐ต=โˆ’1โŽžโŽŸโŽŸโŽ ๎โŽ›โŽœโŽœโŽ๐ผ๐ป๐ต=โˆ’1โŽžโŽŸโŽŸโŽ โŽกโŽขโŽขโŽฃ0๐‘Š๐‘‡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 โ€–ฮฃโ€–2โ„‹2๎€ฝ=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: โ€–ฮฃโ€–2โ„‹2๎€ฝ๐ต=trace๐‘‡0๐‘€โˆ’๐‘‡๐’ฌ22๐‘€โˆ’1๐ต0๎€พ๎€ฝ๐ถ=trace0๐’ซ11๐ถ๐‘‡0๎€พ.(32)

Proof. From Lemma 3, โ€–ฮฃโ€–2โ„‹2๎€ฝ๐ต=trace๐‘‡๎€พโŽงโŽชโŽจโŽชโŽฉ๎‚ƒ0๎€ท๐‘€๐’ฌ๐ต=traceโˆ’1๐ต0๎€ธ๐‘‡๎‚„โŽกโŽขโŽขโŽฃ๐’ฌ11๐’ฌ12๐’ฌ๐‘‡12๐’ฌ22โŽคโŽฅโŽฅโŽฆโŽกโŽขโŽขโŽฃ0๐‘€โˆ’1๐ต0โŽคโŽฅโŽฅโŽฆโŽซโŽชโŽฌโŽชโŽญ๎€ฝ๐ต=trace๐‘‡0๐‘€โˆ’๐‘‡๐’ฌ22๐‘€โˆ’1๐ต0๎€พ,(33)โ€–ฮฃโ€–2โ„‹2โŽงโŽชโŽจโŽชโŽฉ๎‚ƒ๐ถ=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, โ€–ฮฃโ€–2โ„‹2=๐ถ0๐’ซ11๐ถ๐‘‡0=๐ถ0๐‘Š1๐‘†๐‘Š๐‘‡1๐ถ๐‘‡0=๐‘ง๐‘‡๐‘†๐‘ง=๐œŽ1๐‘ง21+๐œŽ2๐‘ง22+โ‹ฏ+๐œŽ๐‘›๐‘ง2๐‘›.(42)

In reduced system, ๎๐ถ0=๐ถ0๐‘Š1๐‘Ÿ, ๎๐’ซ11=๐‘†1. Therefore, โ€–โ€–ฮฃredโ€–โ€–2โ„‹2=๎๐ถ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, โ€–ฮฃโ€–2โ„‹2โˆ’โ€–โ€–ฮฃredโ€–โ€–2โ„‹2=๐œŽ๐‘Ÿ+1๐‘ง2๐‘Ÿ+1+๐œŽ๐‘Ÿ+2๐‘ง2๐‘Ÿ+2+โ‹ฏ+๐œŽ๐‘›๐‘ง2๐‘›โ‰ค๐œŽ๐‘Ÿ+1โ€–๐‘งโ€–22=๐œŽ๐‘Ÿ+1โ€–โ€–๐ถ0โ€–โ€–22.(44)

This is not the โ„‹2-norm of error system โ€–ฮฃโˆ’ฮฃredโ€–โ„‹2, 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 โ€–๐œ€โ€–2โ€–โ„‹โ€–2โ€–โ€–๎ฮฃโ€–โ€–=โˆถฮฃโˆ’โ„‹2โ€–ฮฃโ€–โ„‹2,(45)

tab1
Table 1: Errors of reduced models produced by four different methods.

๐‘› 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.

302498.fig.001
Figure 1: Frequency response on clamped beam model with size ๐‘›=174 and the reduced model of size ๐‘Ÿ=20 by four different methods.
302498.fig.002
Figure 2: Frequency response on clamped beam model with size ๐‘›=174 and the reduced model of size ๐‘Ÿ=29 by four different methods.

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 โ€–โ€–๎ฮฃโ€–โ€–ฮฃโˆ’2โ„‹2๎€ฝ๐ถโ‰ค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 โ€–โ€–ฮฃ๐‘’๐‘Ÿ๐‘Ÿโ€–โ€–2โ„‹2=โ€–โ€–๎ฮฃโ€–โ€–ฮฃโˆ’2โ„‹2โ‰ค๐‘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๐ถ02๎€ธโ€–โ€–2+โ€–โ€–๐ถ02โ€–โ€–22.(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๐’ฌ22โŽžโŽŸโŽŸโŽ โŽ›โŽœโŽœโŽ0๐ผโˆ’๐พ๐‘€โˆ’๐บ๐‘€โŽžโŽŸโŽŸโŽ +โŽกโŽขโŽขโŽฃ๐ถ๐‘‡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.

References

  1. K. Glover, โ€œAll optimal Hankel-norm approximations of linear multivariable systems and their Lโˆž error bounds,โ€ International Journal of Control, vol. 39, no. 6, pp. 1115โ€“1193, 1984. View at Publisher ยท View at Google Scholar ยท View at Scopus
  2. D. G. Meyer and S. Srinivasan, โ€œBalancing and model reduction for second-order form linear systems,โ€ IEEE Transactions on Automatic Control, vol. 41, no. 11, pp. 1632โ€“1644, 1996. View at Google Scholar ยท View at Scopus
  3. Y. Chahlaoui, D. Lemonnier, A. Vandendorpe, and P. Van Dooren, โ€œSecond-order balanced truncation,โ€ Linear Algebra and Its Applications, vol. 415, no. 2-3, pp. 373โ€“384, 2006. View at Publisher ยท View at Google Scholar ยท View at Scopus
  4. Y. Chahlaoui, K. Gallivan, A. Vandendorpe, and P. Van Dooren, โ€œModel reduction of second-order systems,โ€ in Dimension Reduction of Large-Scale Systems, P. Benner, V. Mehrmann, and D. Sorensen, Eds., vol. 45 of Lecture Notes in Computational Science and Engineering, pp. 149โ€“170, Springer, Berlin, Germany, 2005. View at Google Scholar
  5. T. Reis and T. Stykel, โ€œBalanced truncation model reduction of second-order systems,โ€ Mathematical and Computer Modelling of Dynamical Systems, vol. 14, no. 5, pp. 391โ€“406, 2008. View at Publisher ยท View at Google Scholar ยท View at Scopus
  6. A. J. Laub, M. T. Heath, C. C. Paige, and R. C. Ward, โ€œComputation of system balancing transformations and other applications of simultaneous diagonalization algorithms,โ€ IEEE Transactions on Automatic Control, vol. 32, pp. 115โ€“122, 1987. View at Google Scholar
  7. A. C. Antoulas, Approximation of Large-Scale Dynamical Systems, Advances in Design and Control, SIAM, Philadelphia, PA, USA, 2005.
  8. M. G. Safanov and R. Y. Chiang, โ€œA Schur method for balance-trucation model reduction,โ€ IEEE Transactions on Automatic Control, vol. 34, no. 7, pp. 729โ€“733, 1989. View at Publisher ยท View at Google Scholar
  9. Y. Chahlaoui and P. Van Dooren, โ€œBenchmark examples for model reduction of linear time-invariant dynamical systems,โ€ in Dimension Reduction of Large-Scale Systems, P. Benner, V. Mehrmann, and D. Sorensen, Eds., vol. 45 of Lecture Notes in Computational Science and Engineering, pp. 381โ€“395, Springer, Berlin, Germany, 2005. View at Google Scholar
  10. D. C. Sorensen and A. C. Antoulas, โ€œGramians of structured systems and an error bound for structure-preserving model reduction,โ€ in Dimension Reduction of Large-Scale Systems, P. Benner, V. Mehrmann, and D. Sorensen, Eds., vol. 45 of Lecture Notes in Computational Science and Engineering, pp. 117โ€“130, Springer, Berlin, Germany, 2005. View at Google Scholar