Second-Order Model Reduction Based on Gramians

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.


Introduction
Models of linear dynamical systems are often given in second-order form M q(t) + G q(t) + Kq(t) = B 0 u(t) y(t) = C 0 q(t) + D 0 q(t), (1) where M, G, K ∈ R n×n , B 0 ∈ R n×m , C 0 , and D 0 ∈ R p×n are given matrices, u(t) ∈ R m is the given vector of inputs, y(t) ∈ R p is the unknown vector of outputs, q(t) ∈ R n is the unknown vector of internal variables, n is the dimension of the system, and M is assumed to be invertible.Models of mechanical systems in particular are usually of second-order form (1).For such a system, M = M T , G and K = K T are respectively the mass, damping, and stiffness matrices, u(t) ∈ R m is the input, B 0 u(t) = f (t) ∈ R n is the vector of external forces, and q(t) ∈ R n is the vector of internal variables.
Many applications lead to very large models where n is very large, while m, p n. 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: M q(t) + G q(t) + Kq(t) = B 0 u(t) y(t) = C 0 q(t) + D 0 q(t), (2) such that (1) y − y / u < tol, that is, outputs corresponding to the same inputs are close.
(2) The reduced system preserves important system properties such as stability and passivity.
By letting x = q q , second-order system (1) is then transformed to a corresponding 2n-dimensional first-order system Σ : where Since algorithms and theory are well developed for firstorder model reduction, it is natural to use these techniques to develop algorithms for second-order systems.
Under assumption that first-order system is stable, it is well known that the reachability and observability Gramians P and Q are the unique solutions to Lyapunov equation (5), where both P and Q are symmetric positive definite.This Gramian P has also the following variational interpretation proposed in [1].
The optimal value to optimization problem is x * 0 P −1 x 0 , that is, the minimal energy required to steer the state of the system from t = −∞ to state x 0 at time t = 0 is x * 0 P −1 x 0 .Similarly, the optimal value to its dual system is Partition P and Q into four equal blocks: In [2][3][4], it is shown that the optimal value to optimization problem is q * 0 P −1 11 q 0 , that is, the minimal energy required to reach the given position q 0 over all past inputs and initial velocities.And the optimal value to optimization problem min is q * 0 P −1 22 q0 , that is, the minimal energy required to reach the given velocity q0 over all past inputs and initial positions.In [3,4], P 11 and P 22 are defined to be position Gramian and velocity Gramian, respectively.By duality, Q 11 and Q 22 are position Gramian and velocity Gramian, respectively corresponding to Q.
It is well known that Suppose x(t) is the solution of first-order system for impulse input u(t) = δ(t).Then For second-order system and its corresponding first-order system, Therefore, where P 11 and P 22 are position and velocity Gramians, respectively.Similar results can be applied to Q, Q 11 , and Q 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 H 2 -norm.Compared to existing techniques, they are simple and more suitable for largescale 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.

Balanced Truncation Method
For first-order system Σ : balanced truncation method is to transform the system to another coordinate system such that both P and Q are equal and diagonal: where It then truncates the model by keeping the states corresponding to k largest eigenvalues of the product P Q, that is, k 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 P and Q are is where U C is upper triangular and L O is lower triangular matrices.Take SVD of Let Σ k be the first k by k principle submatrix of Σ, U k , and V k consist of the first k columns of U and V , respectively.Then the balanced projection matrices are The reduced system is obtained in: Balanced truncation method has a global error bound [1]: For second-order system there are two balanced truncation methods.One was given in [5] which is to balance both P 11 and Q 11 to form the projection matrices W 1r and V 1r and then get the reduced system by keeping r largest eigenvalues of P 11 Q 11 .It is equivalent to performing projections to its corresponding first-order system as following: The reduction rules are then Similar results can be applied to P 22 and Q 22 .This gives the following algorithm.
(2) Compute the balanced truncation matrices W 1r , are the r largest eigenvalues of P 11 Q 11 .
(3) Perform projection to (M, G, K, B 0 , C 0 , D 0 ) and get ( M, G, K, B 0 , C 0 , D 0 ): (24) In [5], it also gave two other similar methods.One is to balance P 11 and Q 22 to get the projection matrices, the other is to balance P 22 and Q 11 to get the projection matrices.
Another second-order balanced truncation method is given in [3] which is to balance both P 11 and Q 11 to get the projection matrices W 1r and V 1r by keeping r largest eigenvalues of P 11 Q 11 , balance both P 22 and Q 22 to get the projection matrices W 2r and V 2r , then use as projection matrices for the corresponding first-order system: In order to let reduced system have the companion form of second-order system, it then takes the following transformation by letting The corresponding algorithm is as follows: Algorithm 2 (Balanced truncation method based on P11 P22 and Q11 Q22 [3]).(1) Compute P 11 , P 22 , Q 11 , and Q 22 .
(2) Compute the balanced truncation matrices for P 11 and Q 11 to get W 1r , V 1r ∈ R n×r , and balanced transformation matrices for P 22 and Q 22 to get W 2r , V 2r ∈ R n×r .

Some New Algorithms Based on SVD
In this section, we propose some new algorithms based on SVD method directly.First, H 2 norm is an important quantity for bounding linear systems.
Lemma 3 (see [7]).Given a first-order system suppose P and Q are reachability and observability Gramians.
Then, H 2 norm of the system can be expressed as H 2 norm for second-order system (1) with D 0 = 0 can then be stated as in the following proposition.Proposition 4. Given a second-order system (1) with D 0 = 0, suppose is the transformed first-order system, where Assume reachability and observability Gramians P and Q are divided into four equal blocks, and . Then the following holds: Proof.From Lemma 3, From ( 8), we know the positions which are difficult to reach are spanned by the eigenvectors corresponding to small eigenvalues of P 11 .We would like to eliminate those positions in reduced systems.Proposition 4 shows that P 11 can independently decide the system H 2 norm.Therefore, we may get the reduced system by keeping the positions corresponding to r largest eigenvalues of P 11 .This motivates the following procedure for model order reduction.Suppose SVD of P 11 is P 11 = W 1 SV T  1 , where S = diag(σ 1 , σ 2 , . . ., σ n ), and are orthogonal matrices and W 1 = V 1 since P 11 is symmetric positive definite.Let the transformation matrices for the corresponding first-order system be and the projection matrices be where W 1r consists of the first r columns of W 1 .For simplicity, we denote W k and V k by W and V , respectively.Then use W and V as projection matrices for the corresponding firstorder system, that is, perform projections as in ( 22) and ( 23).Same results can be applied to Q 22 .From the dual system, symmetrically Q 11 and P 22 are also crucial in weighting the system H 2 norms independently.This results in the following algorithms.
(2) Take SVD on P and get P = W 1 S 1 W T 1 .Form the orthogonal projection matrices W 1r which consists of the first r columns of W 1 .
(3) Perform projection to get the reduced system ( M, G, K, B 0 , C 0 , D 0 ): 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 P 22 corresponding to small eigenvalues.This motivates the following procedure for model order reduction.Suppose SVD of P 11 is and V 2 are all orthogonal matrices, and Let the transformation matrices for the corresponding firstorder system be and the projection matrices be where W 1r and W 2r consist of the first r columns of W 1 and W 2 , respectively.By using similar idea to Algorithm 2, in order to let the reduced system A have the companion form of second-order model, we adopt a matrix H such that H −1 = W T 1r V 2r , and then perform the projections as in ( 26) and ( 27).Similar results can be applied to Q 11 and Q 22 .This results in the following algorithm.
(2) Take SVD on P 11 and P 22 to get P 11 = W 1 S 1 W T 1 and P 22 = W 2 S 2 W T 2 , respectively.(3) Perform projection to get reduced system ( M, G, K, B 0 , C 0 , D 0 ): where Now, consider Algorithm 5 and SISO second-order systems.Note that by taking SVD, Then from Proposition 4, (42) In reduced system, C 0 = C 0 W 1r , P 11 = S 1 .Therefore, In many cases, the eigenvalues of P 11 decay rapidly.Therefore, the reduced system produced by Algorithm 5 takes the main part of the system H 2 norm.Actually we also have, This is not the H 2 -norm of error system Σ − Σ red H2 , 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 P 11 in Algorithm 5, after computing P 11 which takes O(n 3 ) time, it then uses one-time O(n 3 ) operation which is SVD on P 11 to get the projection matrices.In our Algorithm 6, after computing P 11 and P 22 each of which takes O(n 3 ) time, it then uses twice O(n 3 ) operations which are SVD on P 11 and SVD on P 22 to get the projection matrices.In balanced truncation method in Algorithm 1, after computing P 11 and Q 11 , it also has one-time SVD.Besides this, it uses two Cholesky factorizations O which all take time O(n 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.

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 P 11 in Algorithm 5 ("SVD1"), SVD method on P 11 and P 22 in Algorithm 6 ("SVD12"), balanced truncation on P 11 and Q 11 in Algorithm 1 ("BT1"), and balanced truncation on P 11 and Q 11 , P 22 and Q 22 in Algorithm 5 ("BT12").The comparison is based on the relative error of H 2 norm, where n is the size of M, r is the size of M in reduced model, and m is the size of input vector u(t), p is the size of output vector y(t).
Figures 1 and 2 show the amplitude of frequency response of original and reduced systems for clamped beam model for r = 20 and r = 29, respectively.

Error Bounds and Discussions
Error bounds for first-order model reduction do not apply for second-order systems.Consider second-order system (1) with D 0 = 0 : (46) Suppose the reduced system Σ is obtained by keeping r largest eigenvalues of the orthogonal eigenspace of P 11 .Sorensen and Antoulas in [10] provided a priori error bound showing that the H 2 norm of error system is bounded by a constant times the summation of neglected singular values of P 11 .It considers structured systems in [10].Let Q(s) and P(s) be polynomial matrices in s: where Q is invertible, and Q −1 P is a strictly proper rational matrix.Denote by Q(d/dt), P(d/dt) the differential operators The structured systems are defined by the following equations: where C ∈ R p×n .The Gramian is defined as the Gramian of x(t) when the input is an impulse: Let be the eigensystem of P , where V is orthogonal, the diagonal elements of S are in decreasing order: where where V 1 consists of the first r columns of V .The reduced model is derived from and the reduced system is then Partition accordingly, Theorem 7 (see [10]).The reduced model Σ derived from the dominant eigenspace of the Gramian P for Σ as described above satisfies where K is a modest constant depending on Σ and Σ, and the diagonal elements of S 2 are the neglected smallest eigenvalues of P .
Specially for second-order system (1) with D 0 = 0: it is easy to see that it can be described by (49) with Q(s) = Ms 2 +Gs+K and P(s) = B 0 .So there is the following corollary.

Now, we explain the corollary and what the constant
If the reduced system has no poles on the imaginary axis, sup ω∈R W(iω) 2 is finite, and so sup ω∈R (C 01 W(iω)) * (C 01 W(iω) − 2C 02 ) 2 is finite.Therefore, c 0 is a constant depending on Σ and Σ.There is a similar result for Q 11 by duality.Consider Algorithm 5, SVD method on P 11 .Note that and For second-order system (1) with D 0 = 0, the following theorem provides some observations on structures of gramians P and Q in corresponding first-order system.Theorem 9. Given second-order system (1) with D 0 = 0, then in any coordinate system, there is no state space transformation which gives a block diagonal Q.

Proof. From the observability Lyapunov equation
one can get Equating each block gives From the first equation in (63), one can get Q 12 is not zero, otherwise C 0 = 0, and then from (61) the solution for this Lyapunov equation is Q = ∞ 0 e A T τ C T Ce Aτ dτ = 0.This contradicts with Q being symmetric positive definite.So Q is not block diagonal.
Algorithm 2 [3], that is, balanced truncation method on (P 11 , P 22 ) and (Q 11 , Q 22 ), is to balance both P 11 and Q 11 to get the projection matrices W 1r and V 1r by keeping r largest eigenvalues of P 11 Q 11 , balance both P 22 and Q 22 to get the projection matrices W 2r and V 2r , then use W = as projection matrices for the corresponding first-order system.This is equivalent to assuming both P and Q are block diagonal in corresponding first-order system, that is, in order to keep r largest eigenvalues of P 11 Q 11 and P 22 Q 22 in reduced corresponding first-order system.Similar results also apply for methods in Algorithm 6, SVD method on Q 11 and Q 22 .From Theorem 9, Q 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 (P 11 , Q 11 ) (or (P 22 , Q 22 ) resp.)) and Algorithm 5 (SVD on P 11 (or P 22 , Q 11 , Q 22 resp.)),avoid this drawback in assuming Q 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.

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 P 11 or Q 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.

Figure 1 :Figure 2 :
Figure 1: Frequency response on clamped beam model with size n = 174 and the reduced model of size r = 20 by four different methods.

Table 1 :
Errors of reduced models produced by four different methods.
So the reduction rules in Algorithm 5 are exactly the same as (53).Similar result applies for Q 11 .Therefore, the error systems produced by our Algorithm 5, SVD method on P 11 or Q 11 , all have good global error bounds.So far, there is no priori error bound for error systems based on P 22 or Q 22 .Even though balanced truncation methods in Algorithms 1 and 2, and SVD method on both P 11 and P 22 (or Q 11 and Q 22 ) in Algorithm 6, all have good performance in those three numerical examples, so far no priori error bounds are provided for these methods.