About this Journal Submit a Manuscript Table of Contents 
Journal of Control Science and Engineering
Volume 2012 (2012), Article ID 302498, 9 pages
doi:10.1155/2012/302498
Research Article

Second-Order Model Reduction Based on Gramians

Cong Teng

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) ‖ 𝑦 − ̂ 𝑦 ‖ / ‖ 𝑢 ‖ < t o l , 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 m i n 𝑢 𝐽 ( 𝑢 , − ∞ , 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: ⎛ ⎜ ⎜ ⎝ 𝒫 𝒫 = 1 1 𝒫 1 2 𝒫 2 1 𝒫 2 2 ⎞ ⎟ ⎟ ⎠ ⎛ ⎜ ⎜ ⎝ 𝒬 , 𝒬 = 1 1 𝒬 1 2 𝒬 2 1 𝒬 2 2 ⎞ ⎟ ⎟ ⎠ . ( 7 )

In [2–4], it is shown that the optimal value to optimization problem m i n ̇ 𝑞 0 ∈ ℝ 𝑛 m i n 𝑢 𝐽 ( 𝑢 , − ∞ , 0 ) 𝑠 . 𝑡 . 𝑀 ̈ 𝑞 ( 𝑡 ) + 𝐺 ̇ 𝑞 ( 𝑡 ) + 𝐾 𝑞 ( 𝑡 ) = 𝐵 𝑢 ( 𝑡 ) , 𝑞 ( 0 ) = 𝑞 0 ( 8 ) is 𝑞 ∗ 0 𝒫 − 1 1 1 𝑞 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 problem m i n 𝑞 0 ∈ ℝ 𝑛 m i n 𝑢 𝐽 ( 𝑢 , − ∞ , 0 ) s . t . 𝑀 ̈ 𝑞 ( 𝑡 ) + 𝐺 ̇ 𝑞 ( 𝑡 ) + 𝐾 𝑞 ( 𝑡 ) = 𝐵 𝑢 ( 𝑡 ) , ̇ 𝑞 ( 0 ) = ̇ 𝑞 0 ( 9 ) is ̇ 𝑞 ∗ 0 𝒫 − 1 2 2 ̇ 𝑞 0 , that is, the minimal energy required to reach the given velocity ̇ 𝑞 0 over all past inputs and initial positions. In [3, 4], 𝒫 1 1 and 𝒫 2 2 are defined to be position Gramian and velocity Gramian, respectively. By duality, 𝒬 1 1 and 𝒬 2 2 are position Gramian and velocity Gramian, respectively corresponding to 𝒬 .

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

Suppose 𝑥 ( 𝑡 ) is the solution of first-order system for impulse input 𝑢 ( 𝑡 ) = 𝛿 ( 𝑡 ) . Then 𝒫 = ∞ 0 𝑥 ( 𝑡 ) 𝑥 ∗ ( 𝑡 ) 𝑑 𝑡 . ( 1 1 ) For second-order system and its corresponding first-order system, 𝒫 = ∞ 0 𝑥 ( 𝑡 ) 𝑥 ∗ = ( 𝑡 ) 𝑑 𝑡 ∞ 0 ⎛ ⎜ ⎜ ⎝ ⎞ ⎟ ⎟ ⎠ 𝑞 𝑞 ( 𝑡 ) ̇ 𝑞 ( 𝑡 ) ∗ ( 𝑡 ) ̇ 𝑞 ∗ = ( 𝑡 ) 𝑑 𝑡 ∞ 0 ⎛ ⎜ ⎜ ⎝ 𝑞 ( 𝑡 ) 𝑞 ∗ ( 𝑡 ) 𝑞 ( 𝑡 ) ̇ 𝑞 ∗ ( 𝑡 ) ̇ 𝑞 ( 𝑡 ) 𝑞 ∗ ( 𝑡 ) ̇ 𝑞 ( 𝑡 ) ̇ 𝑞 ∗ ⎞ ⎟ ⎟ ⎠ ( 𝑡 ) 𝑑 𝑡 . ( 1 2 )

Therefore, 𝒫 1 1 = ∞ 0 𝑞 ( 𝑡 ) 𝑞 ∗ 𝒫 ( 𝑡 ) 𝑑 𝑡 , 2 2 = ∞ 0 ̇ 𝑞 ( 𝑡 ) ̇ 𝑞 ∗ ( 𝑡 ) 𝑑 𝑡 , ( 1 3 ) where 𝒫 1 1 and 𝒫 2 2 are position and velocity Gramians, respectively. Similar results can be applied to 𝒬 , 𝒬 1 1 , and 𝒬 2 2 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 Σ ∶ ̇ 𝑥 ( 𝑡 ) = 𝐴 𝑥 ( 𝑡 ) + 𝐵 𝑢 ( 𝑡 ) , 𝑦 ( 𝑡 ) = 𝐶 𝑥 ( 𝑡 ) . ( 1 4 ) balanced truncation method is to transform the system to another coordinate system such that both 𝒫 and 𝒬 are equal and diagonal: 𝜎 𝒫 = 𝒬 = d i a g 1 , 𝜎 2 , … , 𝜎 𝑁 , ( 1 5 )

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 𝒫 = 𝑈 𝐶 𝑈 𝑇 𝐶 , 𝒬 = 𝐿 𝑂 𝐿 𝑇 𝑂 , ( 1 6 )

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

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 . ( 1 8 )

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

Balanced truncation method has a global error bound [1]: ‖ ‖ Σ e r r ‖ ‖ ℋ ∞ = ‖ ‖ Σ − Σ r e d ‖ ‖ ℋ ∞ 𝜎 ≤ 2 𝑘 + 1 + ⋯ + 𝜎 𝑛 . ( 2 0 )

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

there are two balanced truncation methods. One was given in [5] which is to balance both 𝒫 1 1 and 𝒬 1 1 to form the projection matrices 𝑊 1 𝑟 and 𝑉 1 𝑟 and then get the reduced system by keeping 𝑟 largest eigenvalues of 𝒫 1 1 𝒬 1 1 . 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 𝑟 . ( 2 2 ) The reduction rules are then 𝑀 = 𝐼 𝑟 , 𝐾 = 𝑊 𝑇 1 𝑟 𝑀 − 1 𝐾 𝑉 1 𝑟 , 𝐺 = 𝑊 𝑇 1 𝑟 𝑀 − 1 𝐺 𝑉 1 𝑟 , 𝐵 0 = 𝑊 𝑇 1 𝑟 𝑀 − 1 𝐵 0 , 𝐶 0 = 𝐶 0 𝑉 1 𝑟 , 𝐷 0 = 𝐷 0 𝑉 1 𝑟 . ( 2 3 ) Similar results can be applied to 𝒫 2 2 and 𝒬 2 2 . This gives the following algorithm.

Algorithm 1 (Balanced truncation method based on 𝒫 1 1 and 𝒬 1 1 (or 𝒫 2 2 and 𝒬 2 2 , resp.) [5]). (1) Compute 𝒫 1 1 and 𝒬 1 1 .(2) Compute the balanced truncation matrices 𝑊 1 𝑟 , 𝑉 1 𝑟 ∈ ℝ 𝑛 × 𝑟 such that 𝑊 𝑇 1 𝑟 𝑉 1 𝑟 = 𝐼 𝑟 and 𝑊 𝑇 1 𝑟 𝒫 1 1 𝑊 1 𝑟 = 𝑉 𝑇 1 𝑟 𝒬 1 1 𝑉 1 𝑟 = d i a g ( 𝜎 1 , 𝜎 2 , … , 𝜎 𝑟 ) , where 𝜎 2 1 , 𝜎 2 2 , … , 𝜎 2 𝑟 are the 𝑟 largest eigenvalues of 𝒫 1 1 𝒬 1 1 . (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 𝑟 . ( 2 4 )

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

Another second-order balanced truncation method is given in [3] which is to balance both 𝒫 1 1 and 𝒬 1 1 to get the projection matrices 𝑊 1 𝑟 and 𝑉 1 𝑟 by keeping 𝑟 largest eigenvalues of 𝒫 1 1 𝒬 1 1 , balance both 𝒫 2 2 and 𝒬 2 2 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 𝑟 . ( 2 5 )

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 𝑟 𝐻 . ( 2 6 ) The corresponding algorithm is as follows:

Algorithm 2 (Balanced truncation method based on 𝒫 1 1 𝒫 2 2 and 𝒬 1 1 𝒬 2 2 [3]). (1) Compute 𝒫 1 1 , 𝒫 2 2 , 𝒬 1 1 , and 𝒬 2 2 .(2) Compute the balanced truncation matrices for 𝒫 1 1 and 𝒬 1 1 to get 𝑊 1 𝑟 , 𝑉 1 𝑟 ∈ ℝ 𝑛 × 𝑟 , and balanced transformation matrices for 𝒫 2 2 and 𝒬 2 2 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 𝑟 𝐻 , ( 2 7 ) 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 ⎛ ⎜ ⎜ ⎝ 𝐴 Σ = 𝐵 𝐶 ⎞ ⎟ ⎟ ⎠ , ( 2 8 ) suppose 𝒫 and 𝒬 are reachability and observability Gramians. Then, ℋ 2 norm of the system can be expressed as ‖ Σ ‖ 2 ℋ 2 = t r a c e 𝐶 𝒫 𝐶 𝑇 𝐵 = t r a c e 𝑇 . 𝒬 𝐵 ( 2 9 )

ℋ 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 ⎛ ⎜ ⎜ ⎝ 𝐴 Σ = 𝐵 𝐶 ⎞ ⎟ ⎟ ⎠ ( 3 0 ) is the transformed first-order system, where ⎡ ⎢ ⎢ ⎣ 𝐴 = 0 𝐼 − 𝑀 − 1 𝐾 − 𝑀 − 1 𝐺 ⎤ ⎥ ⎥ ⎦ , ⎡ ⎢ ⎢ ⎣ 0 𝑀 𝐵 = − 1 𝐵 0 ⎤ ⎥ ⎥ ⎦ , 𝐶 𝐶 = 0 0 . ( 3 1 ) Assume reachability and observability Gramians 𝒫 and 𝒬 are divided into four equal blocks, and 𝒫 = 𝒫 1 1 𝒫 1 2 𝒫 𝑇 1 2 𝒫 2 2 , 𝒬 = 𝒬 1 1 𝒬 1 2 𝒬 𝑇 1 2 𝒬 2 2 . Then the following holds: ‖ Σ ‖ 2 ℋ 2 𝐵 = t r a c e 𝑇 0 𝑀 − 𝑇 𝒬 2 2 𝑀 − 1 𝐵 0 𝐶 = t r a c e 0 𝒫 1 1 𝐶 𝑇 0 . ( 3 2 )

Proof. From Lemma 3, ‖ Σ ‖ 2 ℋ 2 𝐵 = t r a c e 𝑇 ⎧ ⎪ ⎨ ⎪ ⎩ 0 𝑀 𝒬 𝐵 = t r a c e − 1 𝐵 0 𝑇 ⎡ ⎢ ⎢ ⎣ 𝒬 1 1 𝒬 1 2 𝒬 𝑇 1 2 𝒬 2 2 ⎤ ⎥ ⎥ ⎦ ⎡ ⎢ ⎢ ⎣ 0 𝑀 − 1 𝐵 0 ⎤ ⎥ ⎥ ⎦ ⎫ ⎪ ⎬ ⎪ ⎭ 𝐵 = t r a c e 𝑇 0 𝑀 − 𝑇 𝒬 2 2 𝑀 − 1 𝐵 0 , ( 3 3 ) ‖ Σ ‖ 2 ℋ 2 ⎧ ⎪ ⎨ ⎪ ⎩ 𝐶 = t r a c e 0 0 ⎡ ⎢ ⎢ ⎣ 𝒫 1 1 𝒫 1 2 𝒫 𝑇 1 2 𝒫 2 2 ⎤ ⎥ ⎥ ⎦ ⎡ ⎢ ⎢ ⎣ 𝐶 𝑇 0 0 ⎤ ⎥ ⎥ ⎦ ⎫ ⎪ ⎬ ⎪ ⎭ 𝐶 = t r a c e 0 𝒫 1 1 𝐶 𝑇 0 . ( 3 4 )

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

and the projection matrices be 𝑊 𝑘 = ⎛ ⎜ ⎜ ⎝ 𝑊 1 𝑟 𝑊 1 𝑟 ⎞ ⎟ ⎟ ⎠ , 𝑉 𝑘 = 𝑊 𝑘 , ( 3 6 )

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 𝒬 2 2 . From the dual system, symmetrically 𝒬 1 1 and 𝒫 2 2 are also crucial in weighting the system ℋ 2 norms independently. This results in the following algorithms.

Algorithm 5 (Second-order model reduction—SVD on 𝒫 1 1 (or, 𝒫 2 2 , 𝒬 1 1 , and 𝒬 2 2 , resp.)). (1) Compute 𝒫 1 1 (or 𝒫 2 2 , 𝒬 1 1 , and 𝒬 2 2 , 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 𝑟 . ( 3 7 )

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 𝒫 2 2 corresponding to small eigenvalues. This motivates the following procedure for model order reduction. Suppose SVD of 𝒫 1 1 is 𝒫 1 1 = 𝑊 1 𝑆 1 𝑉 𝑇 1 , where 𝑆 1 = d i a g ( 𝜎 1 , 𝜎 2 , … , 𝜎 𝑛 ) , 𝜎 1 ≥ 𝜎 2 ≥ ⋯ ≥ 𝜎 𝑛 > 0 . And SVD of 𝒫 2 2 is 𝒫 2 2 = 𝑊 2 𝑆 2 𝑉 𝑇 2 , where 𝑆 2 = d i a g ( 𝛿 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 ⎞ ⎟ ⎟ ⎠ , 𝑉 = 𝑊 , ( 3 8 )

and the projection matrices be 𝑊 𝑘 = ⎛ ⎜ ⎜ ⎝ 𝑊 1 𝑟 𝑊 2 𝑟 ⎞ ⎟ ⎟ ⎠ , 𝑉 𝑘 = 𝑊 𝑘 , ( 3 9 )

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 𝒬 1 1 and 𝒬 2 2 . This results in the following algorithm.

Algorithm 6 (Second-order model reduction—SVD on 𝒫 1 1 and 𝒫 2 2 (or 𝒬 1 1 and 𝒬 2 2 , resp.)). (1) Compute 𝒫 1 1 and 𝒫 2 2 .(2) Take SVD on 𝒫 1 1 and 𝒫 2 2 to get 𝒫 1 1 = 𝑊 1 𝑆 1 𝑊 𝑇 1 and 𝒫 2 2 = 𝑊 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 𝑟 𝐻 , ( 4 0 ) where 𝐻 − 1 = 𝑊 𝑇 1 𝑟 𝑊 2 𝑟 .

Now, consider Algorithm 5 and SISO second-order systems. Note that by taking SVD,   𝒫 1 1 = 𝑊 1 𝑆 1 𝑊 𝑇 1 . Denote 𝑆 1 1 = d i a g ( 𝜎 1 , 𝜎 2 , … , 𝜎 𝑟 ) and 𝑆 1 2 = d i a g ( 𝜎 𝑟 + 1 , 𝜎 𝑟 + 2 , … , 𝜎 𝑛 ) . So 𝑆 1 = d i a g ( 𝑆 1 1 , 𝑆 1 2 ) . Suppose 𝑧 ∶ = 𝑊 𝑇 1 𝐶 𝑇 0 , that is,  ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ 𝑧 1 𝑧 2 ⋮ 𝑧 𝑛 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ = 𝑊 𝑇 1 ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ 𝑐 1 𝑐 2 ⋮ 𝑐 𝑛 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ . ( 4 1 ) Then from Proposition 4, ‖ Σ ‖ 2 ℋ 2 = 𝐶 0 𝒫 1 1 𝐶 𝑇 0 = 𝐶 0 𝑊 1 𝑆 𝑊 𝑇 1 𝐶 𝑇 0 = 𝑧 𝑇 𝑆 𝑧 = 𝜎 1 𝑧 2 1 + 𝜎 2 𝑧 2 2 + ⋯ + 𝜎 𝑛 𝑧 2 𝑛 . ( 4 2 )

In reduced system, 𝐶 0 = 𝐶 0 𝑊 1 𝑟 , 𝒫 1 1 = 𝑆 1 . Therefore, ‖ ‖ Σ r e d ‖ ‖ 2 ℋ 2 = 𝐶 0 𝒫 1 1 𝐶 𝑇 0 = 𝐶 0 𝑊 1 𝑟 𝑆 1 𝑊 𝑇 1 𝑟 𝐶 𝑇 0 = 𝜎 1 𝑧 2 1 + 𝜎 2 𝑧 2 2 + ⋯ + 𝜎 𝑟 𝑧 2 𝑟 . ( 4 3 )

In many cases, the eigenvalues of 𝒫 1 1 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 − ‖ ‖ Σ r e d ‖ ‖ 2 ℋ 2 = 𝜎 𝑟 + 1 𝑧 2 𝑟 + 1 + 𝜎 𝑟 + 2 𝑧 2 𝑟 + 2 + ⋯ + 𝜎 𝑛 𝑧 2 𝑛 ≤ 𝜎 𝑟 + 1 ‖ 𝑧 ‖ 2 2 = 𝜎 𝑟 + 1 ‖ ‖ 𝐶 0 ‖ ‖ 2 2 . ( 4 4 )

This is not the ℋ 2 -norm of error system ‖ Σ − Σ r e d ‖ ℋ 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 𝒫 1 1 in Algorithm 5, after computing 𝒫 1 1 which takes 𝑂 ( 𝑛 3 ) time, it then uses one-time 𝑂 ( 𝑛 3 ) operation which is SVD on 𝒫 1 1 to get the projection matrices. In our Algorithm 6, after computing 𝒫 1 1 and 𝒫 2 2 each of which takes 𝑂 ( 𝑛 3 ) time, it then uses twice 𝑂 ( 𝑛 3 ) operations which are SVD on 𝒫 1 1 and SVD on 𝒫 2 2 to get the projection matrices. In balanced truncation method in Algorithm 1, after computing 𝒫 1 1 and 𝒬 1 1 , it also has one-time SVD. Besides this, it uses two Cholesky factorizations 𝒫 1 1 = 𝑈 𝐶 𝑈 𝑇 𝐶 and 𝒬 1 1 = 𝐿 𝑂 𝐿 𝑇 𝑂 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 𝒫 1 1 in Algorithm 5 (“SVD1”), SVD method on 𝒫 1 1 and 𝒫 2 2 in Algorithm 6 (“SVD12”), balanced truncation on 𝒫 1 1 and 𝒬 1 1 in Algorithm 1 (“BT1”), and balanced truncation on 𝒫 1 1 and 𝒬 1 1 , 𝒫 2 2 and 𝒬 2 2 in Algorithm 5 (“BT12”). The comparison is based on the relative error of ℋ 2 norm, where ‖ 𝜀 ‖ 2 ‖ ℋ ‖ 2 ‖ ‖ Σ ‖ ‖ = ∶ Σ − ℋ 2 ‖ Σ ‖ ℋ 2 , ( 4 5 )

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 𝑟 = 2 0 and 𝑟 = 2 9 , respectively.

302498.fig.001
Figure 1: Frequency response on clamped beam model with size 𝑛 = 1 7 4 and the reduced model of size 𝑟 = 2 0 by four different methods.
302498.fig.002
Figure 2: Frequency response on clamped beam model with size 𝑛 = 1 7 4 and the reduced model of size 𝑟 = 2 9 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 𝑞 ( 𝑡 ) . ( 4 6 )

Suppose the reduced system Σ is obtained by keeping 𝑟 largest eigenvalues of the orthogonal eigenspace of 𝒫 1 1 . 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 𝒫 1 1 .

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

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

The structured systems are defined by the following equations: 𝑄 𝑑 Σ ∶ 𝑑 𝑑 𝑡 𝑥 ( 𝑡 ) = 𝑃 𝑑 𝑡 𝑢 ( 𝑡 ) 𝑦 ( 𝑡 ) = 𝐶 𝑥 ( 𝑡 ) , ( 4 9 ) where 𝐶 ∈ ℝ 𝑝 × 𝑛 . The Gramian is defined as the Gramian of 𝑥 ( 𝑡 ) when the input is an impulse: 𝒫 ∫ ∶ = ∞ 0 𝑥 ( 𝑡 ) 𝑥 ( 𝑡 ) ∗ 𝑑 𝑡 . ( 5 0 )

Let 𝒫 = 𝑉 𝑆 𝑉 ∗ ( 5 1 )

be the eigensystem of 𝒫 , where 𝑉 is orthogonal, the diagonal elements of 𝑆 are in decreasing order: ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ 𝜎 𝑆 = 1 𝜎 2 ⋱ 𝜎 𝑛 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ , ( 5 2 )

where 𝜎 1 ≥ 𝜎 2 ≥ ⋯ ≥ 𝜎 𝑛 . Partition 𝑉 = [ 𝑉 1 , 𝑉 2 ] where 𝑉 1 consists of the first 𝑟 columns of 𝑉 . The reduced model is derived from 𝑄 𝑗 = 𝑉 ∗ 1 𝑄 𝑗 𝑉 1 , 𝑃 𝑗 = 𝑉 ∗ 1 𝑃 𝑗 , 𝐶 = 𝐶 𝑉 1 , ( 5 3 ) and the reduced system is then 𝑄 𝑑 Σ ∶ 𝑃 𝑑 𝑑 𝑡 ̂ 𝑥 ( 𝑡 ) = 𝑑 𝑡 𝑢 ( 𝑡 ) ̂ 𝑦 ( 𝑡 ) = 𝐶 ̂ 𝑥 ( 𝑡 ) . ( 5 4 ) Partition accordingly, 𝐶 1 , 𝐶 2 ⎛ ⎜ ⎜ ⎝ 𝑆 = 𝐶 𝑉 , 𝑆 = 1 𝑆 2 ⎞ ⎟ ⎟ ⎠ . ( 5 5 )

Theorem 7 (see[10]). The reduced model Σ derived from the dominant eigenspace of the Gramian 𝒫 for Σ as described above satisfies ‖ ‖ Σ ‖ ‖ Σ − 2 ℋ 2 𝐶 ≤ t r a c e 2 𝑆 2 𝐶 ∗ 2 𝑆 + 𝒦 t r a c e 2 , ( 5 6 ) 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 𝑞 ( 𝑡 ) , ( 5 7 ) 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 𝒫 1 1 = 𝑆 1 𝑆 2 is diagonal, where 𝑆 1 = d i a g ( 𝜎 1 , 𝜎 2 , … , 𝜎 𝑟 ) , 𝑆 2 = d i a g ( 𝜎 𝑟 + 1 , 𝜎 𝑟 + 2 , … , 𝜎 𝑛 ) , and 𝜎 1 ≥ 𝜎 2 ≥ ⋯ ≥ 𝜎 𝑛 . Suppose the reduced system Σ is obtained by keeping 𝑟 largest eigenvalues of 𝒫 1 1 as described in (53) and (54). Then ‖ ‖ Σ 𝑒 𝑟 𝑟 ‖ ‖ 2 ℋ 2 = ‖ ‖ Σ ‖ ‖ Σ − 2 ℋ 2 ≤ 𝑐 0 𝑆 t r a c e 2 = 𝑐 0 𝜎 𝑟 + 1 + 𝜎 𝑟 + 2 + … + 𝜎 𝑛 , ( 5 8 ) where 𝑐 0 is a modest constant depending on Σ and Σ .

Now, we explain the corollary and what the constant 𝑐 0 is. Note 𝑄 ( 𝑠 ) = 𝑀 𝑠 2 + 𝐺 𝑠 + 𝐾 ∈ ℂ 𝑛 × 𝑛 and 𝑄 ( 𝑠 ) = 𝑄 1 1 ( 𝑠 ) 𝑄 1 2 𝑄 ( 𝑠 ) 2 1 ( 𝑠 ) 𝑄 2 2 ( 𝑠 ) with 𝑄 1 1 ( 𝑠 ) ∈ ℂ 𝑟 × 𝑟 . Let 𝑊 ( 𝑠 ) = 𝑄 1 1 ( 𝑠 ) − 1 𝑄 1 2 ( 𝑠 ) . Partition accordingly 𝐶 0 = [ 𝐶 0 1 , 𝐶 0 2 ] . Then, 𝑐 0 = s u p 𝜔 ∈ ℝ ‖ ‖ 𝐶 0 1 𝑊 ( 𝑖 𝜔 ) ∗ 𝐶 0 1 𝑊 ( 𝑖 𝜔 ) − 2 𝐶 0 2 ‖ ‖ 2 + ‖ ‖ 𝐶 0 2 ‖ ‖ 2 2 . ( 5 9 )

If the reduced system has no poles on the imaginary axis, s u p 𝜔 ∈ ℝ ‖ 𝑊 ( 𝑖 𝜔 ) ‖ 2 is finite, and so s u p 𝜔 ∈ ℝ ‖ ( 𝐶 0 1 𝑊 ( 𝑖 𝜔 ) ) ∗ ( 𝐶 0 1 𝑊 ( 𝑖 𝜔 ) − 2 𝐶 0 2 ) ‖ 2 is finite. Therefore, 𝑐 0 is a constant depending on Σ and Σ . There is a similar result for 𝒬 1 1 by duality.

Consider Algorithm 5, SVD method on 𝒫 1 1 . Note that 𝒫 1 1 = ∞ 0 𝑞 ( 𝑡 ) 𝑞 ∗ ( 𝑡 ) 𝑑 𝑡 , ( 6 0 )

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

Even though balanced truncation methods in Algorithms 1 and 2, and SVD method on both 𝒫 1 1 and 𝒫 2 2 (or 𝒬 1 1 and 𝒬 2 2 ) 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 , ( 6 1 ) one can get ⎛ ⎜ ⎜ ⎝ 0 − 𝐾 𝑇 𝑀 𝐼 − 𝐺 𝑇 𝑀 ⎞ ⎟ ⎟ ⎠ ⎛ ⎜ ⎜ ⎝ 𝒬 1 1 𝒬 1 2 𝒬 𝑇 1 2 𝒬 2 2 ⎞ ⎟ ⎟ ⎠ + ⎛ ⎜ ⎜ ⎝ 𝒬 1 1 𝒬 1 2 𝒬 𝑇 1 2 𝒬 2 2 ⎞ ⎟ ⎟ ⎠ ⎛ ⎜ ⎜ ⎝ 0 𝐼 − 𝐾 𝑀 − 𝐺 𝑀 ⎞ ⎟ ⎟ ⎠ + ⎡ ⎢ ⎢ ⎣ 𝐶 𝑇 0 0 ⎤ ⎥ ⎥ ⎦ 𝐶 0 0 = 0 . ( 6 2 ) Equating each block gives − 𝐾 𝑇 𝑀 𝒬 𝑇 1 2 − 𝒬 1 2 𝐾 𝑀 + 𝐶 𝑇 0 𝐶 0 = 0 , − 𝐾 𝑇 𝑀 𝒬 2 2 + 𝒬 1 1 − 𝒬 1 2 𝐺 𝑀 𝒬 = 0 , 1 2 − 𝐺 𝑇 𝑀 𝒬 2 2 + 𝒬 𝑇 1 2 − 𝒬 2 2 𝐺 𝑀 = 0 . ( 6 3 ) From the first equation in (63), one can get 𝒬 1 2 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 ( 𝒫 1 1 , 𝒫 2 2 ) and ( 𝒬 1 1 , 𝒬 2 2 ) , is to balance both 𝒫 1 1 and 𝒬 1 1 to get the projection matrices 𝑊 1 𝑟 and 𝑉 1 𝑟 by keeping 𝑟 largest eigenvalues of 𝒫 1 1 𝒬 1 1 , balance both 𝒫 2 2 and 𝒬 2 2 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, ⎛ ⎜ ⎜ ⎝ 𝒫 𝒫 = 1 1 𝒫 2 2 ⎞ ⎟ ⎟ ⎠ ⎛ ⎜ ⎜ ⎝ 𝒬 , 𝒬 = 1 1 𝒬 2 2 ⎞ ⎟ ⎟ ⎠ , ( 6 4 )

in order to keep 𝑟 largest eigenvalues of 𝒫 1 1 𝒬 1 1 and 𝒫 2 2 𝒬 2 2 in reduced corresponding first-order system. Similar results also apply for methods in Algorithm 6, SVD method on 𝒬 1 1 and 𝒬 2 2 . 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 ( 𝒫 1 1 , 𝒬 1 1 ) (or ( 𝒫 2 2 , 𝒬 2 2 ) resp.)) and Algorithm 5 (SVD on 𝒫 1 1 (or 𝒫 2 2 , 𝒬 1 1 , 𝒬 2 2 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 𝒫 1 1 or 𝒬 1 1 . 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 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.
  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.
  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.
  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.