Model Reduction for Second‐Order Dynamical Systems Revisited

Second‐order control systems are used to describe the dynamics of mechanical and vibrational systems, and in particular, their response to excitations. In this article, we discuss model order reduction (MOR) of such systems. A particular focus is on preserving the second‐order structure and physical properties such as stability and passivity.


Introduction
In this article, we consider linear mechanical structures including some external forcing. These can often be modeled as linear second-order systems of the form with M, D, K ∈ R n×n , B ∈ R n×m , and C p , C v ∈ R p×n , where M and K are symmetric and positive definite. In many real-life problems, n can become unfeasibly large for control, optimization, or simulations, which causes a demand for a good approximation of the system. On the one hand, the aim is that the reduced system should inherit the second-order structurẽ withM ,D,K ∈ R r×r ,B ∈ R r×m , andC p ,C v ∈ R p×r , as well as preserving physically meaningful system properties such as asymptotic stability and passivity. On the other hand, one wishes to have a quality measure for the approximation. This can be an upper bound on some norm of the error system. In systems theory, norms classically used for this purpose are the H ∞ -norm, forcing the L 2 /L 2 -gain of the error system to be small, or the H 2 -norm. Another error measure (which is not a norm) is the so-called gap metric [1]. For the latter, a small upper bound guarantees a good approximation of the closed subspace of the stable trajectories of the system, i. e., trajectories ( u y ) ∈ L 2 . Considering damping optimization problems, prominent reduction techniques used in the literature are based on modal truncation [2]. These, however, fail to guarantee the passivity of the reduced-order model. Moreover, modal truncation does not take into account the input/output structure of the system. Some progress on that was made in [3,4] and preservation of physical properties in the reduced-order model is possible for collocated inputs and outputs. Further second-order reduction methods have been developed in [5][6][7][8][9][10]. However, all the approaches mentioned above are lacking of theoretical results. Neither preservation of stability nor error bounds have been established. Here, we try to address these issues, at least partially, with new methods. Section 2 shows an extension of the first-order limited balanced truncation approaches to second-order systems, while in Section 3 the recovery of the system structure is discussed. Finally, in Section 4, we take a look on a greedy structure-preserving interpolation method.

Second-order frequency-and time-limited balanced truncation
The practice has shown, that very often a global approximation of the system behavior in either the frequency or the time domain is not required. Restricting to specified frequency and time ranges leads to the frequency-and time-limited balanced truncation methods. Rewriting the second-order system (1) into the first-order form allows to combine the ideas for first-order frequency-and time-limited balanced truncation [11] with the second-order balancing approaches [3,6,9]. To describe the approach, we consider the realization (3). The corresponding frequency-limited Gramians P Ω and Q Ω are the unique positive semi-definite solutions of the indefinite Lyapunov equations where the frequency-dependent right-hand sides are given as can also be described as the solutions of where the time-dependent right-hand sides are Taking the sub-Gramian matrices as in [3,6,9] gives the possibility to get a reduced-order system in secondorder form. Implementations of the different resulting second-order frequency-and time-limited balanced truncation methods can be found in MORLAB [12] for medium-size dense systems and in [13] for large-scale sparse systems using Krylov subspace methods.
As a numerical example, we use the triple chain oscillator from [14]. For clarity of presentation, only the results for the approach based on [6] are shown and we restrict ourselves to the frequency-limited method for the approximation between 10 −3 and 10 −2 Hz. The system has been set up with 400 masses per chain, which results in a state-space dimension of n = 1201. The damping behavior between the masses is modeled by Rayleigh damping, i. e., D = αM + βK where we use α = β = 0.003, and the viscosity of additional dampers at the end of each chain is set to ν = 5.
The experiments reported here have been executed with version 4.0 of the MORLAB toolbox [12]. Fig. 1 shows, that in the approximated frequency range, the accuracy of the frequency-limited method (SOFLBT) is higher than for the global approach (SOBT), while both reduced-order models are of order r = 43.

Positive real balanced truncation
In this section, we consider an asymptotically stable second-order system of the form (1), where M, K > 0, D ≥ 0, and the inputs and outputs are collocated, which means C p = 0 and B = C T v have full column rank. Using the Cholesky factorization K = GG T , the system can be rewritten in first-order form aṡ which is internally passive [15], since A + A T ≤ 0 and B = C T . Hence, we can apply positive real balanced truncation [16] to (4). We get the following result. Theorem 3.1 Let a system (4) with minimal realization (A, B, C) be given with its transfer function denoted by G. Then G has a positive real balanced realization A b , B b , B T b and a positive real Gramian Σ with block structure where A 22 ∈ R 2 ×2 has only purely imaginary eigenvalues, To obtain a reduced first-order model of size 2r with r := q+m+ , we truncate the blocks corresponding to σ ± q+1 , . . . , σ ± k . Then we get the reduced system Ã ,B,B T with which has the same block structure as (5). Since we perform positive real balanced truncation, its transfer functionG(s) := B T sI 2r −Ã −1B is also positive real. Further, an a priori gap metric error bound together with an a posteriori H ∞ -error bound for positive real balanced truncation are established in [17].
It turns out that the reduced system does not always preserve the specific structure that is needed to be able to transform it into second-order form. A necessary condition on the so called sign characteristics of the nonzero real zeros of the reduced transfer function has to be fulfilled. Note that, by the block structure (6) of our reduced system, the zeros ofG with negative real part are exactly the eigenvalues of . If they are semi-simple, we can apply [18, Chap. I.5, Thm. 5.3] that states the existence of an invertible matrix T ∈ R 2q×2q such that where Λ 11 ∈ R nc×nc has no real eigenvalues and Λ ± = diag(λ ± 1 , . . . , λ ± nr ) with 0 < λ ± 1 ≤ . . . ≤ λ ± nr . Theorem 3.2 Let the reduced system Ã ,B,B T with transfer functionG be given as above. In order that it has a secondorder realization of the form (2) withB =C T v ,C p = 0,M ,K > 0, andD =D T , it is necessary and sufficient that Λ − < Λ + is satisfied in (7). This condition is always fulfilled if the original system is overdamped, i. e., If the condition of Theorem 3.2 is fulfilled, then using the block structure ofÃ,B and the normal form in (7), we are able to construct matricesB ∈ R r×m and X, Y T ∈ R r×2r such that These matrices allow us to deriveM ,D, andK as in [19] from the moments One can show thatM ,K > 0 and thatD =D T has at most m negative eigenvalues. In particular, for overdamped systems we get thatD ≤ 0. Now, if the reduced system does not fulfill the necessary condition in Theorem 3.2, we can always construct an H ∈ H m×m ∞ such that this condition is fulfilled forĜ := G −1 + H −1 , whereĜ is again passive, asymptotically stable, and the gap metric δ G ,Ĝ is arbitrarily small. Here, the number of additional poles ofĜ is bounded from above by the number of negative zeros ofG.

Greedy H ∞ interpolation
In this section, we briefly describe an interpolatory H ∞ model order reduction scheme for systems with collocated inputs and outputs, i. e., we have C p = 0 and B = C T v in order to be able to preserve symmetry and asymptotic stability by an appropriate www.gamm-proceedings.com choice of the projection spaces. More precisely, we construct a sequence of reduced transfer functions of the form . . , for appropriately chosen projection matrices V j . Our method aims to iteratively reduce the H ∞ -norm of the error transfer function E j := G −G j . To do so, we evaluate E j H∞ := max ω∈R∪{∞} E j (iω) 2 = E j (iω j ) 2 and choose V j+1 := . This choice guarantees Hermite interpolation properties between the original and reduced transfer functions G andG j at the interpolation points iω 1 , iω 2 , . . . , iω j , such that the error near these points becomes small. This procedure is repeated until a specified error tolerance is met.
The bottleneck of this algorithm is the repeated computation of the H ∞ -norm of the error transfer function, which is expensive to evaluate since the error system is of large dimension. However, the repeated computation of the H ∞ -norm of transfer functions of such large-scale systems has been made possible by the methods presented in [20,21] which we use here.  In Figure 2, we illustrate the effectiveness of our algorithm using a triple chain oscillator example [14] of order n = 1000 and an error bound of 10 −6 . It can be observed that the error at the single interpolation points is zero. Compared with two different second-order balanced truncation (SOBT) approaches from [3], the greedy interpolation approach has a slightly larger maximal error with the same reduced order r = 29. However, in contrast to SOBT, we obtain full information on the current error and may terminate whenever the reduced-order model satisfies a prescribed error bound. Moreover, the method also allows for an easy adaptation to frequency-limited reduction. We currently investigate post-processing strategies to improve the performance of the greedy interpolation such as an additional optimization of the interpolation points. Furthermore, our approach may be combined with the subspaces obtained from (SO)BT to initialize the first projection space.