Balanced truncation model order reduction in limited time intervals for large systems

Balanced truncation model order reduction in limited time intervals for large systems Adv Comput Math (2018) 44:1821–1844 https://doi.org/10.1007/s10444-018-9608-6 Balanced truncation model order reduction in limited time intervals for large systems Patrick Kurschner ¨ Received: 7 July 2017 / Accepted: 16 April 2018 / Published online: 5 June 2018 © The Author(s) 2018 Abstract In this article we investigate model order reduction of large-scale systems using time-limited balanced truncation, which restricts the well known balanced trun- cation framework to prescribed finite time intervals. The main emphasis is on the efficient numerical realization of this model reduction approach in case of large sys- tem dimensions. We discuss numerical methods to deal with the resulting matrix exponential functions and Lyapunov equations which are solved for low-rank approx- imations. Our main tool for this purpose are rational Krylov subspace methods. We also discuss the eigenvalue decay and numerical rank of the solutions of the Lyapunov equations. These results, and also numerical experiments, will show that depending on the final time horizon, the numerical rank of the Lyapunov solutions in time- limited balanced truncation can be smaller compared to standard balanced truncation. In numerical experiments we test the approaches for computing low-rank factors of the involved Lyapunov solutions and illustrate that time-limited balanced truncation can generate reduced order models having a higher accuracy in the considered time region. Keywords Lyapunov equation · Rational Krylov subspaces · Model order reduction · Balanced truncation · Matrix exponential Mathematics Subject Classification 2010 15A16 · 15A18 · 15A24 · 65F60 · 93A15 · 93C Communicated by: Peter Benner Patrick Kurschner ¨ kuerschner@mpi-magdeburg.mpg.de Computational Methods in Systems and Control Theory, Max Planck Institute for Dynamics of Complex Technical Systems, Sandtorstr. 1, Magdeburg, 39106, Germany 1822 P. Kurschner ¨ 1 Introduction 1.1 Model reduction of linear systems Consider continuous-time, linear, time-invariant (LTI) systems x( ˙ t) = Ax(t ) + Bu(t ), x(0) = 0, (1a) y(t) = Cx(t) (1b) n×n n×m p×n with A ∈ R , B ∈ R ,and C ∈ R . Typically, the vector functions x(t) ∈ n m p R , u(t ) ∈ R ,and y(t) ∈ R are referred to as state, control, and, output vector, respectively. We assume that m, p  n and A is Hurwitz, i.e., (A) ⊂ C ,such that (1) is asymptotically stable. Given a large state space dimension n, model order reduction aims toward finding a reduced order model ˙ ˜ ˜ x( ˜ t) = Ax( ˜ t) + Bu(t ), x( ˜ 0) = 0, (2a) y( ˜ t) = Cx( ˜ t) (2b) r×r r×m p×r r with A ∈ R , B ∈ R , C ∈ R , x( ˜ t) ∈ R and a drastically reduced state dimension r  n. The smaller reduced system (2) should approximate the input- output behavior or the original system (1) well. Moreover, simulating the system, i.e., solving the differential equations in (2) for many different control functions u(t ) should be numerically much less expensive compared to the original system (1). For the approximation of (1) regarding the actual time domain behavior, it is desired that for all feasible input functions u(t ), y( ˜ t) ≈ y(t) for t ≥ 0, (3a) in other words, y(t)−˜ y(t) should be small ∀t ≥ 0 in some norm ·. With the help of the Laplace transformation, one can also formulate the approximation problem in the frequency domain, e.g., via H(ıω) ≈ H(ıω) for ω ∈ R,ı =−1, −1 −1 ˜ ˜ ˜ ˜ where H(s) = C(sI − A) B, H(s) = C(sI − A) B (3b) are the transfer function matrices of (1)and(2). There exist different model order reduction technologies and here we focus on balanced truncation (BT) [36] model order reduction. The backbone of BT are the infinite controllability and observabil- ity Gramians of (1) which are the the symmetric, positive semidefinite solutions P ,Q of the continuous-time, algebraic Lyapunov equations ∞ ∞ T T T T AP + P A =−BB ,A Q + Q A =−C C. (4) ∞ ∞ ∞ ∞ The Hankel singular values (HSV) of (1) are the square roots of the eigenvalues of the product P Q and are system invariants under state spac transformations. ∞ ∞ The magnitude of the HSV enables to identify components that are weakly con- trollable and observable. In BT this is achieved by first transforming (1)intoa balanced realization such that P = Q =  = diag (σ ,...,σ ). Dropping all ∞ ∞ ∞ 1 n states corresponding to small σ gives the reduced order model. Solving (4)forthe j Balanced truncation model order reduction in limited time... 1823 Algorithm 1: Square-root balanced truncation with low-rank factors Input : System matrices A, B, C defining an asymptotically stable dynamical system (1). ˜ ˜ ˜ Output: Matrices A, B, C of the reduced system. 1 Compute Z ,Z (e.g., with the methods described in [9, 46]), such that P Q T T Z Z ≈ P , Z Z ≈ Q in (4). P ∞ Q ∞ P Q 2 Compute thin singular value decomposition T T X X Y Y Z Z = XY = diag ( , ) P 1 2 1 2 1 2 with  = diag (σ ,...,σ ) containing the largest r (approximate) HSV. 1 1 r 1 1 − − 2 2 3 Construct T := Z Y  and S := Z X  . P 1 Q 1 1 1 4 Generate reduced order model T T ˜ ˜ ˜ A := S AT , B := S B, C := CT . (5) Gramians P ,Q is the computationally most challenging part of balanced trun- ∞ ∞ cation. For large-scale systems one therefore uses low-rank approximations of the T n×k Gramians instead, e.g., Z Z ≈ P with low-rank solution factors Z ∈ R , P ∞ P rank(Z ) = k  n, and likewise for Q . This strategy is backed up by the often P P ∞ numerically observed and theoretically explained rapid eigenvalue decay of solutions of Lyapunov equations [2, 3, 27, 37, 48]. The computation of the low-rank factors Z ,Z of P ,Q can be done efficiently by state-of-the-art numerical algorithms P Q ∞ ∞ for solving large Lyapunov equations [9, 46]. BT using low-rank factors Z ,Z of P Q the Gramians (4) is illustrated in Algorithm 1. The H -norm H  = sup(H(ıω) ). H 2 ω∈R of a stable system (1)isthe L induced norm of the convolution operator. It con- nects to the time-domain behavior via y ≤H  u . With exact Gramian L L 2 ∞ 2 T T factors, i.e., Z Z = P , Z Z = Q , BT is known to always generate a P ∞ Q ∞ P Q asymptotically stable reduced system for which the error bound H − H  ≤ 2 σ (6) H j j =r+1 holds. 1.2 Model reduction in limited time- and frequency intervals The approximation paradigms (3) enforce that the reduced system (2) is accurate for all times t ∈ R and frequencies ω ∈ R. From a practical point of view, achieving (3) might overshoot a realistic objective. For instance, if (1) models a mechanical or 1824 P. Kurschner ¨ electrical system, practitioners working with this model (and its approximation (2)) might only be interested in certain finite frequency intervals 0 ≤ ω <ω < ∞. 1 2 Likewise, when the final goal is to carry out simulations of (1), i.e., acquire time- domain solutions for various controls u(t ), one is usually only interested in y(t) for t smaller than some final time t < ∞. Hence, we consider time- and frequency restricted versions of (3) of the form y( ˜ t) ≈ y(t) for t ∈ T ⊂ R , (7a) H(iω) ≈ H(iω) for ω ∈  ⊂ R, =−, (7b) where the time- and frequency regions T , of interest should be provided by the underlying application. The expressions (7) demand that the reduced order model (2) is only accurate in T , but allow larger errors outside these regions. Compared to MOR approaches for the unrestricted setting (3), it is desired that MOR approaches for (7) achieve smaller approximation errors in T , with the same reduced order r. Alternatively, one demands that time- and frequency restricted MOR leads to com- parable approximation errors in T , with reduced systems of smaller order r.A secondary question is how the added time- and frequency restrictions influence the computational effort compared to an unrestricted MOR method of the same type. This issue will be in our particular focus. Typically, the time region in (7a)has the form T =[0,t ],t < ∞ (8a) e e which will also be the main situation in this work, but the more general restriction T =[t ,t ], 0 <t <t < ∞ (8b) s e s e will also be briefly discussed. Regarding the frequency restricted setting (7b), the typical regions are ˆ ˆ ˆ :=  ∪−,  := [ω ,ω ] with 0 ≤ ω < ... < ω <ω < ∞. i i+1 1 h h+1 i=1 Introducing time- or frequency restrictions into balanced truncation MOR has been originally proposed in [25] and further studied in, e.g., [8, 18, 21, 29]. In certain applications, e.g., circuit design, only single frequencies ω ∈ R might be of interest and an associated version of balanced truncation is addressed in [19]. H -MOR with limitations or weights on the frequency and time domain is investigated in, e.g., [11, 26, 32, 38, 39, 47, 49]. As continuation of our work in [8] regarding frequency-limited BT, we consider in this paper the numerically efficient realization of time-limited balanced truncation (TLBT) for large-scale systems. 1.3 Overview of this article We start in Section 2 by reviewing the general concept of time-limited BT from [25], mainly for restrictions of the form (8a). This includes the associated time-limited Gramians as well as the respective Lyapunov equations. Similar to standard BT, Balanced truncation model order reduction in limited time... 1825 executing TLBT for large systems heavily relies on how well the time-limited Grami- ans can be approximated by low-rank factorizations. This issue is investigated in Section 3.1 where we have a particular interest in the question when t induces sig- nificant differences between the infinite and time-limited Gramians. The actual issue of numerically dealing with the arising matrix functions and computing the low-rank factors of the time-limited Gramians is topic of Sections 3.2 and 3.3, respectively. Motivated by the promising results for frequency-limited BT studied in [8], we again employ rational Krylov subspace methods for this task. Section 4 collects different generalizations of TLBT including general state-space and certain differential alge- braic systems, the time restriction (8b), and stability preservation. In Section 5,the proposed numerical approach is tested numerically with respect to the approximation of the Gramians as well as the reduction of the dynamical system. 1.4 Notation The real and complex numbers are denoted by, respectively, R and C, R ,(R ), − + C (C ) are the sets of strictly negative (positive) real numbers and the open left − + (right) half plane. The space of real (complex) Matrices of dimension n × m is n×m n×m R (C ). For any complex quantity X = Re (X) + ı Im (X), we denote by Re (X) and Im (X) its real and, respectively, imaginary parts with ı being the imag- inary unit, and its the complex conjugate is X = Re (X) − ı Im (X). The absolute value of any real or complex scalar is denoted by |z|. Unless stated otherwise, · is the Euclidean vector- or subordinate matrix norm (spectral norm). By A and A = A we indicate the transpose and complex conjugate transpose of a real and n×n −1 complex matrix, respectively. If A ∈ C is a nonsingular, its inverse is A ,and −H H −1 −1 A = (A ) . Expressions of the form x = A b are always to be understood as solving a linear system Ax = b for x. The identity matrix of dimension n is indicated T m by I , and the vector of ones is denoted by 1 := (1,..., 1) ∈ R . The notation n m A 0(≺ 0) indicates symmetric positive (negative) definiteness of a symmetric or Hermitian matrix A,  () refers to semi-definiteness, and A  ()B means A − B  ()0. The spectrum of a matrix A is denoted by (A). 2 Gramians and balanced truncation for finite time horizons Since A is assumed to be Hurwitz, the infinite Gramians P ,Q in (4) can be ∞ ∞ represented in integral form as ∞ ∞ T T At T A t A t T At P = e BB e dt, Q = e C C e dt. (9) ∞ ∞ 0 0 Restricting the integration limits in the integrals (9) to a time interval [t ,t ] s e immediately yields the time-limited Gramians P , Q . T T 1826 P. Kurschner ¨ Definition 1 (Time-limited Gramians [25]) The time-limited reachability and observability Gramians of (1) with respect to the time-interval T =[t ,t ], 0 ≤ t < s e e t < ∞ are defined by t t e e At T A t T T At P = e BB e dt, Q = e C C e dt. (10) T T t t s s Theorem 1 (Lyapunov equations for the time-limited Gramians [25]) The time- limited Gramians P and Q according to Definition 1 are equivalently given in the T T following ways. a) The finite time Gramians P ,Q from (10) are given by T T T T At A t At A t s s e e P = e P e − e P e , (11a) T ∞ ∞ T T A t At A t At s s e e Q = e Q e − e Q e , (11b) ∞ ∞ where P and Q are the infinite reachability and observability Gramians (4). ∞ ∞ b) The time-limited Gramians P ,Q satisfy the time-limited reachability and T T observability Lyapunov equations T T T At At s e AP +P A =−B B + B B ,B := e B, B := e B (12a) t t t t T T s t e t s e s e T T T At At s e A Q +Q A =−C C + C C ,C := C e ,C := C e . (12b) t t t t T T t s t e s e s e Note that the time-limited Gramians (10) also exist for unstable A and if (A) ∩ (−A) =∅ they still solve the Lyapunov (12), see [40]. Hence, TLBT might be one possible approach to reduce unstable systems. Except for one brief numerical experi- ment, we will not pursue this topic any further. In analogy to the infinite time horizon case, the square roots of the eigenvalues of the product P Q are called time-limited T T Hankel singular values. A basic calculation shows that the time-limited Hankel sin- gular values are, as the standard Hankel singular values, invariant with respect to state-space transformations. By comparing the Lyapunov equations for the infinite Gramians (4) with (12), one immediately sees that the only differences are the inho- mogeneities, while the left hand sides are the same unchanged (adjoint) Lyapunov operators. This raises the question how much different the time-limited Gramians are from the infinite ones, and how this depends on the time interval of interest. We will pursue this in the next section, especially regarding the numerically important issue of how well the time-limited Gramians can be approximated by low-rank factoriza- T T tions P ≈ Z Z , Q ≈ Z Z . Of course, before approximately solving T P T Q T P T Q T T At the time-limited Lyapunov equations, the actions of the matrix exponentials e to B A t T (aswellase to C ) have to be dealt with numerically. Numerical approaches for handling the matrix exponential and computing low-rank solution factors Z ,Z P Q T T are topic of Section 3. A square-root version of TLBT is simply carried out by substituting Step 1 of Algorithm 1 by the code snippet shown in Algorithm 2 below and using the low-rank Balanced truncation model order reduction in limited time... 1827 solution factors Z , Z in the remaining steps. Depending on the time region of P Q T T interest, TLBT might in general not be a stability preserving method and, thus, there is also no H -error bound similar to (6). This can be regained by modifying TLBT further [29] and we discuss this issue in Section 4. Without this modification it is possible to establish an error bound in the H -norm [40]. Algorithm 2: Required changes in Algorithm 1 for TLBT At A t i i 1a Compute (approximation of) B := e B, C = C e , i ∈{s, e}. t t i i T T 1b Compute Z ,Z , such that P ≈ Z Z , Q ≈ Z Z in (12). P Q T P T Q T T T T P Q T T 3 Numerical computation of low-rank factors of time-limited Gramians This section is concerned with the actual numerical computation of low-rank factors of the time-limited Gramians. Before algorithms for dealing with the matrix exponen- tials and the Lyapunov equations are discussed, we briefly investigate the numerical ranks of the time-limited Gramians. For the sake of brevity, all considerations will be mostly restricted to the reachability Gramian because the observability Gramian can T T be dealt with similarly by replacing A, B with A ,C , respectively. Moreover, only the situation (8a) will be discussed, i.e., t = 0and 0 <t < ∞. s e 3.1 The difference of the infinite and time-limited Gramians In order to approximate P , Q by low-rank factorizations, it is desirable that their T T eigenvalues decay rapidly. For investigating this decay we assume from now that eigenvalues of symmetric (positive definite) matrices are given in a non-increasing order λ ≥ ... ≥ λ . The inhomogeneities of the time-limited Lyapunov (12)have 1 n up to twice the rank of their unlimited counterparts (4). Hence, by the available theory on the eigenvalue decay of solutions of matrix equations [2, 27, 37, 45] one expects that the eigenvalues of the time-limited Gramians decay somewhat slower than those of the infinite ones. We will see in the numerical experiments that, similar to the frequency-limited Gramians, in most situations it is the opposite case: the eigenvalues of the time-limited Gramians exhibit a faster decay and, consequently, have smaller numerical ranks. Assuming that (A, B) is controllable (rank[A − λI, B]= n, ∀λ ∈ C) it holds P 0 and, using (10), the relation (11a) yields 0  P = P − ∞ T ∞ T T At A t At A t e e e e e P e .Since E(t ) := e P e 0 it holds ∞ e ∞ P  P and also λ (P ) =P  P  = λ (P ). ∞ T 1 ∞ ∞ 2 T 2 1 T Due to the stability of A, the difference E(t ) = P − P will decay for increas- e ∞ ing values of t . Tight bounds for the eigenvalue behavior of P and E(t ) with e T e respect to the parameter t are difficult to derive and is a research topic of its own. Here, for simplicity we restrict the discussion to the case m = 1 and present a basic investigation of how the decay of E(t ) depends on t . e e 1828 P. Kurschner ¨ n −1 Lemma 1 Let B ∈ R , (A, B) controllable, A be diagonalizable, i.e., A = XX , −1 = diag (λ ,...,λ ), and define w := X B, X := X diag (w), N(t ) := 1 n B e t At A t H e e e X e diag (w).Then E(t ) = e P e = N(t )CN(t ) , e ∞ e e H t  t H e e P = P − E(t ) = P − N(t )CN(t ) = X C − e C e X , (13) ∞ e ∞ e e B −1 and C := is a Hermitian positive definite Cauchy matrix. λ +λ i j i,j =1 Proof Apply the eigendecomposition of A and P = X CX from [2, Lemma 3.2] ∞ B to (11a). Consider the impulse response of (1), At t y (t ) = Cη(t ), η(t ) := e B = X e w = N(t)1 , δ n indicating that the impulse-to-state-map η(t) and N(t) decay at a similar rate. With the spectral abscissa R := max Re (λ), the basic point wise bounds λ∈(A) Rt Rt η(t)≤ e Xw, N(t)≤ e Xw make this more visible. Using this and (13), we can conclude that for increasing t, E(t) is getting smaller at a similar speed as η(t). Consequently, a significant differ- ence between P and P might be only observed when t is chosen small enough T ∞ e in the sense that η(t ) has not reached an almost stationary state. The handling of the case m> 1 can be carried out similarly. Moreover, with a similar argumentation the decay rate of the time-limited Hankel singular values can be roughly connected to the decay of y (t ). To conclude, TLBT might be only practicable for small time horizons or for weakly damped systems. A similar investigation regarding TLBT for unstable systems is given in [47]. 3.2 Approximating the products with the matrix exponential There are several numerical approaches available to approximate the action of the matrix exponential to (a couple of) vectors, see, e.g., [1, 6, 12–15, 23, 24, 30, 31, 34, 35, 44]. Here, we are mostly interested in projection methods using (block) rational Krylov subspaces of the form ⎡ ⎤ −1 ⎣ ⎦ RK = span {q ,...,q },q = (A − s I) B (14) k 1 k k j j =1 where s ∈ C ∪ ıR ∪{∞} are called shifts. They represent the poles of a rational k + approximation r = ψ /φ of e with polynomials ψ ,φ of degree at most k k−1 k−1 k−1 k−1 k − 1. Balanced truncation model order reduction in limited time... 1829 n×km Let Q ∈ R have orthonormal columns and range (Q ) = RK .Thena k k k At Galerkin approximation [44]ofe B takes the form H t k e ˆ ˆ B ≈ B := Q B , B := e B , (15) t t ,k k t ,k t ,k k e e e e T T H := Q AQ ,B := Q B. k k k k k Note that for m = 1, the rational function r underlying the rational Krylov H t T zt k e approximation Q e Q B interpolates f(z) = e at (H ) (rational Ritz val- k k ues) [31, Theorem 3.3]. Further information on the approximation properties can be found in, e.g., [6, 10, 14–16, 30]. In the remainder we assume s =∞ s.t. q = B, 1 1 T T km×m m×m range B ⊆ range Q and B =[β , 0] ∈ R ,where β ∈ R .The ( ) ( ) k k orthonormal basis matrix Q and the restriction H can be efficiently computed by k k H t k e a (block) rational Arnoldi process [41]. Since H is low-dimensional, e can be computed by standard dense methods for the matrix exponential [34]. The choice of shifts s (k> 1) is crucial for a rapid convergence and several strategies exists for this purpose [31]. In this work we exclusively use adaptive shift generation tech- niques [14–16] because this appeared to be the safest strategy in the majority of experiments. For the situation m = 1 and after step k of the rational Arnoldi process, the next shift s is selected via k+1 s−z s = argmax |r (s)|,r (s) = , (16) k+1 k k s−s ∂S j =1 where z are the eigenvalues of H (Ritz values of A); s ,j = 1,...,k are the j k j previously used shifts; and ∂S marks a set of discrete points from the boundary of S approximating (A). We follow [16] and use S as the convex hull of (H ). k k k In the symmetric case, one can also use the spectral interval given by the largest and smallest eigenvalue of A [14, 15, 31]. For m> 1 we simply use each s in the denominator of r in (16) m times as in [16]. A different technique to deal with m> 1 includes, e.g., tangential rational Krylov methods [17]. The rational Krylov subspace method was chosen not only because of the good approximation properties, but also because the generated subspace is also a good candidate to acquire low-rank solution factors of (12). 3.3 Computing the low-rank Gramian factors Using (14), (15), a Galerkin approximation of the time-limited Gramian is P = T ,k Q Y Q ≈ P ,where Y solves the projected time-limited Lyapunov equation k k T k T T T ˆ ˆ H Y + Y H =−B B + B B , (17) k k k k t ,k k k t ,k H t ˆ k s (cf. [43]). It holds B = B = e B for the special case t = 0 considered here, k 0,k k s indicating already how t = 0 can be included as well. Hence, after B of sufficient s t ,k accuracy is found, we follow the idea in [8] by recycling the rational Krylov basis and continuing the rational Arnoldi process for (12a) until P leads to a sufficiently T ,k 1830 P. Kurschner ¨ Algorithm 3: Rational Krylov subspace method for time-limited Lyapunov (12a) Input : A, B, t as in (12a), tolerances 0 <τ ,τ  1. e f P n× T Output: Z ∈ R such that Z Z ≈ P with  ≤ mk  n. k k T 1 B = q β s.t. q q = I , Q = q . 1 1 m 1 1 2 for k = 1, 2,... do 3 Get new shift s via, e.g., (16). k+1 4 Solve (A − s I)g = q for g. k+1 k 5 Real, orthogonal expansion of Q : g = Re g − Q (Q Re g ). ( ) ( ) k + k 6 q = g β s.t. q q = I , Q =[Q ,q ]. k+1 + k k+1 m k+1 k k+1 k+1 7 if Im (s ) = 0 then k+1 8 k = k + 1, g = Im g − Q (Q Im g ) ( ) ( ) + k 9 q = g β s.t. q q = I , Q =[Q ,q ]. j +1 + k k+1 m k+1 k k+1 k+1 T T 10 H = Q AQ , B = Q B. k k k k k H t ˆ k e ˆ 11 B = e B , B = Q B t ,k k t ,k k t ,k e e e 12 if B − B /B  <τ then t ,k t ,k−1 t ,k f e e e T T T ˆ ˆ 13 Solve H Y + Y H + B B − B B = 0for Y . k k k k t ,k k k k t ,k T T T T T A(Q Y Q )+(Q Y Q )A +BB −B B k k k k t ,k k k t ,k 14 Set μ := . T T ˆ ˆ B B −B B k t ,k k t ,k 15 if μ <τ then k P T T 16 Y = SS , S S = I ,  = diag (γ ,...,γ ). k km 1 km 17 Truncate if necessary:  = diag (γ ,...,γ ), S = S(:, 1 : ), ≤ mk. 18 return low-rank solution factor Z = Q S . k k small norm of the Lyapunov residual. For the Lyapunov stage, the same adaptive shift generation (16) can be used [16]. The rational Krylov subspace method for generating a low-rank approximation Z Z ≈ P is illustrated in Algorithm 3. k T In the presence of a complex shift s , it is implicitly assumed that s is the k+1 k+1 subsequent shift. For this situation, the orthogonal expansion in Steps 5–9 of the already computed basis matrix Q by real basis vectors goes back to [42]. The pro- jected matrix H in Step 10, as well as the Lyapunov residual norm in Step 14 can be computed without accessing the large matrix A, see, e.g., [10, 16, 31, 41]. The small Lyapunov equation in Step 13 can be solved by standard methods for dense matrix equations, e.g., the Bartels-Stewart [4] method which we employ here. Once the scaled Lyapunov residual norm falls below the desired threshold, the rational Arnoldi process is terminated and the Steps 16–18 bring the computed low-rank Gramian approximation in the desired form Z Z and allow a rank truncation. Note that typi- cally, once the approximation of B is found, the generated subspace is already good enough to acquire a low-rank Gramian approximation without many additional iter- ations of the rational Arnoldi process. Often, the criteria in Steps 12–15 are satisfied in the same iteration step. Balanced truncation model order reduction in limited time... 1831 4 Extensions and further problems 4.1 Multiple time values Computing low-rank factors of the time-limited Gramians (12) for the more general approximation setting (8b) with a nonzero start time t , i.e., T =[t ,t ], can also be s s e done using Algorithm 3 with minor adjustments. Having computed a rational Krylov At basis for approximating e B for some time value t , the same basis typically also provides good approximations for any other time values [31]. Consequently, Algo- H t ˆ k s ˆ rithm 3 has to be simply changed by adding B = e B , B = Q B t ,k k t ,k k t ,k s s s and appropriately adjusting the steps regarding the Gramian approximation. Also the methods relying on Taylor approximations [1] can be easily modified to handle multi- ple time values. However, even if a nonzero t does not yield additional computational complications, this does not imply that TLBT will produce a accurate approximation of the transient behavior of (1)in [t ,t ]. Already the original TLBT paper [25] states s e that TLBT in [t ,t ] is expected to only give good approximations of the impulse s e m At response (u(t ) = vδ(t), v ∈ R ) y (t ) = C e Bv and numerical experiments confirm this. For an intuitive explanation assume for simplicity that no truncation is done in TLBT, i.e. in Step 2 of Algorithm 1. Then it holds range (Q ) = range (T ), At At At s e s e B ≈ B ∈ range (Q ),e B ≈ B ∈ range (Q ), and likewise for C e , k,t k k,t k s e At C e and range (S). Hence, the impulse response is accurately approximated at the relevant times. A(t −τ) For the response of an arbitrary input u(t ), y (t ) = C e Bu(τ )dτ,such argumentation clearly does not automatically hold. The value of y (t ) with respect to a general u(t ) depends, in general, on the values at the times before t. In the present form TLBT restricts the approximation to the time frame T and, thus, y (t ) will be poorly approximated for t ≤ t which, in turn, makes it difficult to acquire good approximations in T . One approach is to apply a time translation to the underlying system (1) such that the requested time-interval is transformed to [0,t − t ]. However, this time trans- e s lation will also introduce an inhomogeneous initial value x , which is an additional difficulty for model order reduction. Some strategies to cope with nonzero initial val- ues are given in [5, 33] and we plan to investigate the incorporation to time-limited BT in the future. 4.2 Generalized state-space and index-one descriptor systems Consider generalized state-space systems Mx( ˙ t) = Ax(t ) + Bu(t ), x(0) = 0, (18a) y(t) = Cx(t) (18b) n×n with M ∈ R nonsingular. Simple manipulations reveal that, similar to unre- stricted [7] and frequency-limited BT [8], the time-limited Gramians of (18)are P , T 1832 P. Kurschner ¨ M Q M,where P , Q solve the time-limited generalized Lyapunov equations T T T −1 T T T T M At −1 AP M + MP A =−B B + B B ,B := M e M B, (19a) T T t t t s t e t i s e −1 T T T T M At A Q M + M Q A =−C C + C C ,C := C e (19b) T T t t t t s t e i s e −1 AM t for t = t ,t . Note that B := e B. Low-rank factors of P , Q can be i s e t T T computed by methods for generalized Lyapunov equations. In particular, the ratio- nal Krylov subspace methods which we employ for approximating B will implicitly −1 −1 T −1 −T −1 work on M A, M B or alternatively, if M = LL 0, on L AL , L B. T T With low-rank approximations P ≈ Z Z , Q ≈ Z Z ,theSVDof T P T Q T P T Q T T Z MZ has to be used in Step 2 of Algorithm 1. In some of the numerical experiments we will encounter the situation M = A A B M 0 n ×n 1 1 2 1 f f C C ,A = ,B = ,C = [ ] with M ∈ R ,A ∈ 1 2 1 4 A A B 3 4 2 n−n ×n−n f f R nonsingular, s.t. (18) becomes a semi-explicit index-one descriptor system. Eliminating the algebraic constraints leads to a general state-space system −1 −1 −1 ˆ ˆ ˆ defined by M , A := A −A A A and B := B −A A B , C := C −C A A , 1 1 2 3 1 2 2 1 2 3 4 4 4 −1 and an additional feed-through term −C A B ,see [22]. Time-limited and unre- 2 2 ˆ ˆ ˆ stricted BT can be applied right away to this system via (19) defined by M , A, B, C. However, the matrix A will in general be dense and, thus, solving the linear sys- tems in Algorithm 3 can be very expensive. In the context of unrestricted BT, the authors of [22] exploit in a LR-ADI iteration for the Gramians that the arising ˆ ˆ ˆ dense linear systems (A − sM )V = W are equivalent to the sparse linear systems T T T T ˆ ˆ (A − sM)V = W with V =[V ,] , W =[W , 0] which are easier to solve numerically. We use the same trick within Step 4 of Algorithm 3. 4.3 Stability preservation and modified TLBT Because of the altered and sometimes indefinite right hand sides of (12), (19), TLBT is in general not stability preserving. Only when the used time interval is long enough such that the right hand sides are negative semi-definite, TLBT will produce an asymptotically stable reduced order model [25, Condition 1]. For the general sit- uation, a stability preserving modification of TLBT is proposed in [29]usingthe Lyapunov equations mod mod T T T AP M + MP A =−B B , mod T T mod (20) T mod T mod T A Q M + M Q A =−C C , mod T T mod 1 1 2 2 B B C C T B := Q diag |λ |,..., |λ | ,C := diag |λ |,..., |λ | Q . mod B mod 1 r 1 r C B C n×r n×r B C where Q ∈ R , Q ∈ R contain the eigenvectors corresponding to the B C r ≤ 2m, r ≤ 2p nonzero eigenvalues λ ,λ of the right hand sides of (12a), (19a) B C i i and, respectively, (12b), (19b). The rational Krylov subspace method in Algorithm 3 for M = I, t = 0 can be easily modified for (20) by replacing Step 13 with the steps shown in Algorithm 4. Balanced truncation model order reduction in limited time... 1833 Algorithm 4: Changes in Algorithm 3 for modified time-limited Gramians 13a Compute partial eigendecomposition T B B T B ββ 0 ˆ ˆ − B B Q = Q diag λ ,...,λ , Q Q = I , λ = 0. t ,k B B B r e r B t ,k 1 B i e B B B 2 13b Factor of projected modified rhs B := Q diag |λ |,..., |λ | . mod,k B 1 B T T 13c Solve H Y + Y H + B B for Y . k k k mod,k k k mod,k Generalization for M = I and 0 <t <t are straightforward. Note that because s e the modified time-limited Gramians do not fulfill a relation of the form (10)or (11), we cannot expect a fast singular value decay similar to P ,Q . As observed T T in [8, 29] for modified frequency-limited BT, modified TLBT might also lead to less accurate approximations in the considered time region compared to unmodified TLBT. 5 Numerical experiments Here, we illustrate numerically the results of Section 3 regarding the numerical rank of the frequency-limited Gramians as well as the numerical method for computing low-rank factors of P ,Q . Afterwards, the quality of the approximations obtained T T by TLBT with the low-rank factors is evaluated and compared against unlimited BT. Further topics like nonzero starting times and the modified TLBT scheme are also examined along the way. All experiments are done in MATLAB 8.0.0.783 ® ® (R2012b) on a Intel Xeon CPU X5650 @ 2.67GHz with 48 GB RAM. Table 1 list the used examples and their properties. For time-domain simulation of (1), (18)and the reduced order models for a given input function u(t ), an implicit midpoint rule with a fixed small time step t is used. Because the impulse response of (1)for m At u(t ) = δ(t)v, v ∈ R , x = 0 can be expressed as y (t ) = C e Bv, it is computed 0 δ by the same integrator applied to the uncontrolled system (u(t ) = 0) with initial con- dition x = Bv. For the impulse, or step responses, the vector distributing the control to the columns of B is set to v = 1 . Table 1 Dimensions, final integration time t , step size t, matrix properties, and source of the test systems Example nmpt t properties source bips 606 7135 4 4 20 0.04 index-1, n = 606, M = I morwiki ,[22] f 1 n bips 3078 21128 4 4 20 0.04 index-1, n = 3078, M = I morwiki ,[22] f 1 n vertstand 16626 6 6 600 0.6 A ≺ 0, M 0, C random morwiki ,[28] rail 79841 7 6 400 0.4 A ≺ 0, M 0 Oberwolfach Collection , ID=38881 https://morwiki.mpi-magdeburg.mpg.de/morwiki http://portal.uni-freiburg.de/imteksimulation/downloads/benchmark 1834 P. Kurschner ¨ Fig. 1 Scaled eigenvalues of P (top left), Hankel and time-limited Hankel singular values λ(P Q ) T T T (bottom left), norm of impulse response y (t ) (top right), and numerical ranks of infinite, (modified) time-limited Gramians against varying final times t (bottom right) for the bips 606 system Because the bips systems have eigenvalues very close to the imaginary axis which caused difficulties for all considered numerical methods, the shifted matrix A − 0.08M is used instead as in [22]. The generalized state-space systems vertstand and rail are dealt with Cholesky factorizations M = LL as explained in Section 4.2. There, the sparse Cholesky factors L are computed by the MATLAB command chol(M,’vector’). Matrix exponentials and Lya- punov equations defined by smaller dense matrices (including the projected ones in Algorithm 3) are solved directly by the expm and lyap routines. 5.1 Decay of the eigenvalues of the Gramians and the time-limited singular values At first we investigate how the end time t influences the eigenvalue decay of the time-limited Gramians. The index one descriptor system bips 606 is used for this experiment and reformulated into an equivalent state space system of dimension n = 606 as explained in Section 4.2. This comparatively small size allows a direct computation of the matrix exponentials and the Gramians. The top left plot in Fig. 1 shows the scaled and ordered eigenvalues λ /λ of the infinite reachability Gramian j 1 P and the time-limited one P for t = 1, 3, 10. Obviously, a distinctly faster ∞ T e eigenvalue decay of P is only observed for small time values t = 1, 3. As the final T e time increases, the eigenvalues move closer to the ones of P . The eigenvalues of ∞ Balanced truncation model order reduction in limited time... 1835 the observability Gramians exhibit a similar behavior. This observation is even more drastic for the decay of the time-limited Hankel singular values shown in the bottom left plot. For the largest value t = 10, hardly any difference to the Hankel singular values is visible. In the top right plot the point wise norm of the impulse response y (t ) shows that after t = 10, y (t ) has already almost reached its stationary phase. δ e δ This confirms the expectation that significant differences between infinite and time- limited Gramians occur only for times t which are small with respect to the behavior of the impulse response. The bottom right plots shows the numerical rank of infinite, time-limited and modified time-limited Gramians against t . The numerical ranks of P ,Q clearly move toward the numerical ranks of P ,Q as t increases. In T T ∞ ∞ e mod mod contrast, the numerical ranks of P ,Q (Section 4.3) are always close to the T T ones of P ,Q and appear to be largely unaffected by different values of t .This ∞ ∞ e is a very similar behavior as observed for the modified frequency-limited Gramians in [8]. 5.2 Computing low-rank factors of the time-limited Gramians We proceed by testing the computation of low-rank factors of the infinite and (mod- ified) time-limited reachability Gramians by the rational Krylov subspace method in Algorithm 3. The stopping criteria for matrix function and Gramian approximations −8 use the thresholds τ = τ = 10 . To save some computational cost, the projected f P matrix exponentials and Lyapunov equations (Steps 11 and 13 in Algorithm 3) are only dealt with every 5th iteration step. Table 2 summarizes the used time values t , the produced subspace dimensions d, the ranks r of the low-rank approximations, and the total computing times t for rk mod the Gramians P , P ,and P . Apparently, larger rational Krylov subspaces and ∞ T approximately twice as long computation times are needed to obtain low-rank factors of the time-limited Gramians, but the final ranks of the low-rank approximations of the time-limited Gramians P are in all cases smaller compared to P . The modified T ∞ mod time-limited Gramians P do not show this behavior as they have ranks similar to the infinite Gramians. For the rail example, increasing the time horizon t does not change the required subspace dimension d for the approximation of the time-limited Table 2 Results of the numerical numerical computation of low-rank factors of the different Gramians: time horizon t , generated subspace dimension d,rank r of the low-rank approximations, and computing time t in seconds rk mod P P P ∞ T Example t dr t dr t dr t e rk rk rk bips 3078 3 308 245 18.33 664 131 105.89 664 241 106.17 vertstand 300 180 164 6.18 300 140 13.42 300 183 13.55 rail 10 245 224 34.71 420 168 74.06 420 228 76.06 100 245 224 34.71 420 198 76.53 420 227 74.68 1836 P. Kurschner ¨ Gramian P , but its rank r is clearly larger. In all cases, no additional steps of the At rational Krylov method were necessary once the approximation of e B was found. The results for the observability Gramians were largely similar. For the bips 3078 example, the observability Gramians appeared to be more demanding for Algorithm 3 than the reachability Gramians. 5.3 Model reduction results Now we execute BT [36]aswellas(modified)TLBT[25, 29] using the square-root method (Algorithm 1) with the low-rank factors of the reachability and observability Gramians computed in the previous section. For this purpose, we restrict ourselves to the reduction to fixed specified orders r. The approximation quality of the con- structed reduced order models is assessed via the point wise and maximal relative error norms y(t)−y (t ) r 2 E (t ) := ,t ≤ t , E := max E (t ) f T y(t) t ∈[0,t ] of the output responses y(t), y( ˜ t) of original and reduced order models. We consider the response y (t ) for the impulse input u(t ) = δ(t)1 for all examples. Moreover, δ m for each used test system, the transient response with respect to an additional input signal u(t ) is also considered. We use step like input signals u(t ) = 1 and u(t ) = 501 for the bips 3078 and, respectively, rail example, and u := [5 · 10 · m ∗ 2 T 0.198(sin(t π/100) ), 4, 2, 1, 3, 1] for vertstand. The results are given in Table 3 listing the largest relative error E in [0,t ] and the overall computation time t , i.e., the computation time for computing low-rank mor factors of both reachability and observability Gramians by Algorithm 3 plus the time for Algorithm 1 to execute the BT variants. We also indicate whether the produced reduced order models are asymptotically stable (s = 1) or unstable (s = 0).For some selected settings of u(t ) and t , the system responses and point wise relative Table 3 Model reduction results by BT, TLBT, and modified TLBT using different t and u(t ): reduced order r, largest relative error E in [0,t ], overall computation time t ,and s ∈{0, 1} indicates if T e mor reduced system is asymptotically stable or unstable BT TLBT mod. TLBT Example ut E t s E t s E t s e T mor T mor T mor bips 3078, δ 3 5.10e-04 291.3 1 1.08e-06 331.9 0 5.15e-04 399.6 1 r = 100 1 3 6.90e-06 291.3 1 6.33e-09 331.9 0 1.20e-05 399.6 1 vertstand, δ 300 8.90e-02 13.2 1 2.44e-02 30.4 1 6.48e-02 31.6 1 r = 20 u 300 8.26e-03 13.2 1 9.18e-04 30.4 1 9.95e-03 31.6 1 rail, δ 10 6.60e-01 62.2 1 5.66e-04 130.8 0 6.59e-01 137.4 1 r = 50 δ 100 6.60e-01 62.2 1 7.90e-03 135.7 1 6.60e-01 136.4 1 50 10 1.78e-03 62.2 1 6.08e-07 130.8 0 2.61e-03 137.4 1 50 100 1.78e-03 62.2 1 2.57e-05 135.7 1 2.15e-03 136.4 1 Balanced truncation model order reduction in limited time... 1837 Fig. 2 Responses of original and reduced systems obtained by different BT versions: (top left) bips 3078, t = 3, u(t ) =impulse, r = 100, (top right) vertstand, t = 300, u(t ) = u , r = 20, e e ∗ (bottom left) rail, t = 10, u(t ) =impulse, r = 50, (bottom left) rail, t = 100, u(t ) = 50, r = 50 e e errors E (t ) are plotted against the time t in Figs. 2 and 3, respectively. Figure 4 shows the behavior of E as the reduced order r increases. Fig. 3 Relative errors E (t ) obtained by different BT versions: (top left) bips 3078, t = 3, u(t ) =impulse, r = 100, (top right) vertstand, t = 300, u(t ) = u , r = 20, (bottom left) rail, e ∗ t = 10, u(t ) =impulse, r = 50, (bottom left) rail, t = 100, u(t ) = 50, r = 50 e e 1838 P. Kurschner ¨ Fig. 4 Maximum relative errors E in [0,t ] against increasing reduced orders r for different u(t ) and T e BT variants: (top left) bips 3078, t = 3, u(t ) =impulse, (top right) vertstand, t = 300, u(t ) = u , e e ∗ (bottom left) rail79k, t = 10, u(t ) =impulse, (bottom left) rail79k, t = 100, u(t ) = 50 e e Apparently, for the chosen orders r and in the time regions of interest, the largest relative errors E of the reduced order models returned by TLBT are in most experi- ments more than one order of magnitude smaller compared to standard and modified time-limited BT. The plots in Fig. 4 also indicate that much larger reduced order models are needed for BT and modified TLBT to achieved the same accuracy as unmodified TLBT. Figure 3 shows that after leaving the time region T , TLBT deliv- ers larger errors. However, Table 3 also reveals that executing TLBT and its modified version is more time consuming than standard BT. This is a direct consequence of the higher computation times for getting the low-rank factors of the (modified) time-limited Gramians which was pointed out in the previous subsection (Table 2). Using the concept of angles between subspaces or the modal assurance crite- T 2 2 2 rion (MAC) (MAC(x, y) =|y x| /(x y ),[20]) indicated that the spaces spanned by the projection matrices T , S in TLBT and T , S in unre- TLBT TLBT BT BT stricted BT, respectively, are different. For example, computing the MAC for the right n ×100 projection matrices T , T ∈ R for the bips 3078 system showed that TLBT BT only a few of the columns of both matrices are well correlated to each other (i.e., MAC(T e ,T e ) ≈ 1 for very few i, j ∈{1,..., 100}). TLBT i BT j Moreover, albeit the higher accuracy in [0,t ], in some cases TLBT produces unstable reduced order models. This is especially visible in the upper left plot of Fig. 2 showing the impulse responses of bips 3078. For times t ≥ t outside the interval [0,t ], the impulse response of reduced order model generated by TLBT exhibits an exponential growth and departs from the original response. As illustrated Balanced truncation model order reduction in limited time... 1839 Table 4 Results of BT, TLBT, and modified TLBT reduction to r = 15 of vertstand example with respect to time frame T =[t ,t ]=[50, 100] s e BT TLBT mod. TLBT input u E t s E t s E t s T mor T mor T mor δ 7.95e-04 14.2 1 3.07e-04 37.6 0 9.21e-04 32.5 1 u 6.34e-03 14.2 1 1.78e-02 37.6 0 6.85e-03 32.5 1 with the rail examples and proven in [25], using a higher end time t can already cure this. Modified TLBT does not generate unstable reduced systems, but its approx- imation quality in [0,t ] is very close to standard BT without time restrictions. Taking also into account the higher computational costs of modified TLBT given in Table 3, the introduction of the time restriction is rendered essentially redundant because no smaller errors are achieved in the targeted time region. Hence, if stability preservation in the reduced order model is crucial, we recommend to stick to standard BT. To conclude, TLBT fulfills the goal to acquire smaller errors in the desired time interval [0,t ], but at the price of somewhat larger execution times because the computation of required low-rank Gramian factors is currently more costly. Now we carry out one experiment to evaluate the approximation qualities of TLBT for nonzero t .The vertstand example is used and the reduced order model should approximate the output y(t) with respect to u(t ) = δ(t)v and u(t ) = u =[5 · 10 · 2 T 0.198(sin(t π/100) ), 4, 2, 1, 3, 1] in the time window T =[t ,t ]=[50, 100]. s e The low-rank factors of the Gramians are computed as before using Algorithm 3 with the required small extensions mentioned in Section 4.1. The results are summarized in Table 4 and Fig. 5 illustrates the relative errors plots as well as the largest relative error in T against the reduced order r. Compared to stan- dard BT, TLBT delivers, as expected, more accurate results of the impulse response, but fails for u(t ) = u . Moreover, the bottom left plot Fig. 5 shows that the decay of E with respect to the impulse response is less monotonic as in the case t = 0s.t. T s for some reduced orders r, especially larger ones, standard BT outperforms TLBT. The failure of TLBT to approximate the transient response for arbitrary inputs u(t ) was also observed in other experiments with different time intervals and examples, and occasionally even for the impulse response. In experiments with smaller systems, using exact matrix exponentials and Gramian factors did not yield any improvement of the reduction results such that the problems with TLBT are most likely not a result of inaccurate Gramian approximations. In summary, although TLBT in the current form does indeed deliver significantly more accurate reduced order models in the time intervals [0,t ], the same cannot be said when time intervals [t ,t ], t > 0are e s e s considered. Improving the performance of TLBT in this scenario is therefore subject of further research. As the final experiment we briefly test the application of TLBT for the reduction of unstable systems. For this purpose a modification of the bips 606 example is ˆ ˆ used with A := A + 0.2M leading to max Re (λ) = 0.323, λ ∈ (A). The final 1840 P. Kurschner ¨ Fig. 5 Results for the reduction of the vertstand example in [t ,t ]=[50, 100]. The top plots show s e the relative error norms against t for u(t ) = δ(t)v (left) and u(t ) = u (right) after a reduction to r = 15. The behavior of E in [t ,t ] for increasing reduced dimensions r is illustrated in the bottom plots T s e time is set to t = 5 and a reduced order model with r = 100 is generated. Due to the comparatively small system dimension (n = 606), the matrix exponentials and Gramians could be computed by direct, dense methods for these tasks. Moreover, the employed rational Krylov methods from the experiments before were not able to compute the required low-rank Gramians factors. The impulse response of exact and reduced system and the relative error are illustrated in Fig. 6. Apparently, TLBT was able to reduce this mildly unstable system to a reduced −2 order model with maximal relative error E ≈ 1.2 · 10 in [0,t ]. Hence, provided T e t is chosen in a reasonable way, e.g., with a sufficient distance to the exponential Fig. 6 Results for the reduction of the unstable variation of the bips 606 example for t = 5, r = 100, u = δ(t)v (left: system response, right: relative error) Balanced truncation model order reduction in limited time... 1841 growth of y(t), TLBT appears to be a potential candidate from model order reduc- tion of unstable systems. However, the applications to large-scale unstable systems is currently still difficult because algorithms for computing low-rank solutions of Lya- punov equations usually require that A is (anti)stable. Advances in this direction are, therefore, necessary to pursue this type of reduction further. 6 Conclusions BT model order reduction for large-scale systems restricted to finite time inter- vals [25] was investigated. The resulting Lyapunov equations that have to be solved numerically also include the matrix exponential in their inhomogeneities. We first showed that the difference of time-limited and infinite Gramians decays for increas- ing times in a similar way as the impulse response of the underlying system. Hence, for small time intervals, a reduced numerical rank of the time-limited Gramians can be observed. Future research should further investigate the influence of the chosen end time on the eigenvalue decay of the Gramians. As in frequency-limited BT [8], we proposed to handle the matrix exponen- tials and the Lyapunov equations by an efficient rational Krylov subspace method incorporating a subspace recycling idea. While this numerical approach already led to satisfactory results, there is some room for improvement, especially regarding the approximation of the action of the matrix exponential. In this context, further work could include, for instance, enhanced strategies like tangential directions [17], different inner products [23], or using different methods [1, 12] altogether. The numerical reduction experiments indicated that TLBT is able to acquire sev- eral orders of magnitude more accurate reduced order models in time intervals of the form [0,t ] at a somewhat higher, although still comparable, numerical effort. Similar techniques were applied for stability preserving modified TLBT [29]which, however, could not keep up with standard or time-limited BT in terms of efficiency and accuracy of the reduced order models. Hence, we recommend to use standard BT if the preservation of stability is an irrevocable goal. For keeping the high accuracy in the time-interval of interest, a different way to make TLBT a stability preserving method has to be found. TLBT for time regions [t ,t ] with nonzero start times t > 0 s e s provided much worse results than with t = 0, except for approximating the impulse response. In the form introduced in [25], TLBT appears to be incapable of produc- ing good reduced order models with respect to [t ,t ] and, thus, further research is s e necessary in this direction. The results in [5, 33] might be one possible ingredient for this. A brief experiment also indicated that TLBT can be employed to reduce unsta- ble systems. The efficient computation of low-rank Gramian factors in this case is currently not as advanced as in the stable situation, making this a further interesting research topic. Acknowledgments Open access funding provided by Max Planck Society. I thank the referees for their helpful comments. Moreover, I am grateful for the constructive discussions with Maria Cruz Varona, Serkan Gugercin, and Stefan Guettel. 1842 P. Kurschner ¨ Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, dis- tribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. References 1. Al-Mohy, A.H., Higham, N.J.: Computing the action of the matrix exponential, with an application to exponential integrators. SIAM J. Sci. Comput. 33(2), 488–511 (2011) 2. Antoulas, A.C., Sorensen, D.C., Zhou, Y.: On the decay rate of Hankel singular values and related issues. Syst. Cont. Lett. 46(5), 323–342 (2002) 3. Baker, J., Embree, M., Sabino, J.: Fast singular value decay for Lyapunov solutions with nonnormal coefficients. SIAM J. Matrix Anal. Appl. 36(2), 656–668 (2015) 4. Bartels, R.H., Stewart, G.W.: Solution of the matrix equation AX + XB = C: Algorithm 432. Comm. ACM 15, 820–826 (1972) 5. Beattie, C., Gugercin, S., Mehrmann, V.: Model reduction for systems with inhomogeneous initial conditions. Syst. Control Lett. 99, 99–106 (2017). https://doi.org/10.1016/j.sysconle.2016.11.007 6. Beckermann, B., Reichel, L.: Error estimates and evaluation of matrix functions via the faber transform. SIAM J. Numer. Anal. 47(5), 3849–3883 (2009). https://doi.org/10.1137/080741744 7. Benner, P.: Solving large-scale control problems. IEEE Control Syst. Mag. 14(1), 44–59 (2004) 8. Benner, P., Kurschner, ¨ P., Saak, J.: Frequency-limited balanced truncation with low-rank approxima- tions. SIAM J. Sci. Comput. 38(1), A471–A499 (2016). https://doi.org/10.1137/15M1030911 9. Benner, P., Saak, J.: Numerical solution of large and sparse continuous time algebraic matrix Ric- cati and Lyapunov equations: a state of the art survey. GAMM Mitteilungen 36(1), 32–52 (2013). https://doi.org/10.1002/gamm.201310003 10. Berljafa, M., Guttel, ¨ S.: Generalized rational krylov decompositions with an application to rational approximation. SIAM J. Matrix Anal. Appl. 36(2), 894–916 (2015) 11. Breiten, T., Beattie, C., Gugercin, S.: Near-optimal frequency-weighted interpolatory model reduc- tion. Sys. Control Lett. 78, 8–18 (2015) 12. Caliari, M., Kandolf, P., Ostermann, A., Rainer, S.: Comparison of software for computing the action of the matrix exponential. BIT 54(1), 113–128 (2014) 13. Davies, P.I., Higham, N.J.: Computing f(A)b for matrix functions f . In: Boric¸i, A., Frommer, A., Joo, ´ B., Kennedy, A., Pendleton, B. (eds.) QCD and Numerical Analysis III, Lect. Notes Comput. Sci. Eng., vol. 47, pp. 15–24. Springer, Berlin (2005). https://doi.org/10.1007/3-540-28504-0 2 14. Druskin, V., Knizhnerman, L., Zaslavsky, M.: Solution of large scale evolutionary problems using rational krylov subspaces with optimized shifts. SIAM J. Sci. Comput. 31(5), 3760–3780 (2009) 15. Druskin, V., Lieberman, C., Zaslavsky, M.: On adaptive choice of shifts in rational krylov subspace reduction of evolutionary problems. SIAM J. Sci. Comput. 32(5), 2485–2496 (2010) 16. Druskin, V., Simoncini, V.: Adaptive rational Krylov subspaces for large-scale dynamical systems. Syst. Control Lett. 60(8), 546–560 (2011) 17. Druskin, V., Simoncini, V., Zaslavsky, M.: Adaptive tangential interpolation in rational Krylov subspaces for MIMO dynamical systems. SIAM J. Matrix Anal. Appl. 35(2), 476–498 (2014). https://doi.org/10.1137/120898784 18. Du, X., Benner, P.: Finite-Frequency Model Order Reduction of Linear Systems via Parameterized Frequency-dependent Balanced Truncation. Tech. Rep. 1602.04408. ArXiv e-prints (2016) 19. Du, X., Benner, P., Yang, G., Ye, D.: Balanced truncation of linear time-invariant systems at a sin- gle frequency. Preprint MPIMD/13-02, Max Planck Institute Magdeburg. Available from http://www. mpi-magdeburg.mpg.de/preprints/ (2013) 20. Ewins, D.: Modal testing: theory, practice, and application (2nd edn) Research Study Press LTD (2000) 21. Fehr, J., Fischer, M., Haasdonk, B., Eberhard, P.: Greedy-based approximation of frequency-weighted Gramian matrices for model reduction in multibody dynamics. Z. Angew. Math. Mech. 93(8), 501– 519 (2013) 22. Freitas, F., Rommes, J., Martins, N.: Gramian-based reduction method applied to large sparse power system descriptor models. IEEE Trans. Power Syst. 23(3), 1258–1270 (2008) Balanced truncation model order reduction in limited time... 1843 23. Frommer, A., Lund, K., Szyld, D.B.: Block Krylov subspace methods for computing functions of matrices applied to multiple vectors. Tech. Rep. 17-03-21, Department of Mathematics Temple University (2017) 24. Frommer, A., Simoncini, V.: Matrix functions. In: Schilders, W., van der Vorst, H., Rommes, J. (eds.) Model Order Reduction: Theory, Research Aspects and Applications, Mathematics in Industry, vol. 13, pp. 275–303. Springer, Berlin (2008). https://doi.org/10.1007/978-3-540-78841-6 13 25. Gawronski, W., Juang, J.: Model reduction in limited time and frequency intervals. Int. J. Syst. Sci. 21(2), 349–376 (1990). https://doi.org/10.1080/00207729008910366 26. Goyal, P., Redmann, M.: Towards time-limited H -optimal model order reduction. Tech. Rep. WIAS Preprint No. 2441 (2017) 27. Grasedyck, L.: Existence of a low rank or H -matrix approximant to the solution of a Sylvester equation. Numer. Lin. Alg. Appl. 11, 371–389 (2004) 28. Großmann, K., Stadel, C., Galant, A., Muhl, ¨ A.: Berechnung von Temperaturfeldern an Werkzeug- maschinen. Zeitschrift fu,r ¨ Wirtschaftlichen Fabrikbetrieb 107(6), 452–456 (2012) 29. Gugercin, S., Antoulas, A.C.: A survey of model reduction by balanced truncation and some new results. Internat. J. Control 77(8), 748–766 (2004). https://doi.org/10.1080/00207170410001713448 30. Guttel, ¨ S.: Rational Krylov Methods for Operator Functions. Ph.D. thesis, Technische Uni- versitat ¨ Bergakademie Freiberg, Germany. Available from http://nbn-resolving.de/urn:nbn:de:bsz: 105-qucosa-27645 (2010) 31. Guttel, ¨ S.: Rational Krylov approximation of matrix functions: Numerical methods and optimal pole selection. GAMM-Mitteilungen 36(1), 8–31 (2013). https://doi.org/10.1002/gamm.201310002 32. Halevi, Y.: Frequency weighted model reduction via optimal projection. IEEE Trans. Automat. Control 37(10), 1537–1542 (1992) 33. Heinkenschloss, M., Reis, T., Antoulas, A.C.: Balanced truncation model reduction for systems with inhomogeneous initial conditions. Automatica 47(3), 559–564 (2011). https://doi.org/10.1016/j.auto matica.2010.12.002 34. Higham, N.: Functions of matrices: Theory and Computation. Society for Industrial and Applied Mathematics, Philadelphia (2008) 35. Knizhnerman, L.A.: Calculation of functions of unsymmetric matrices using Arnoldi’s method. Comput. Math. Math. Phys. 31(1), 1–9 (1992) 36. Moore, B.C.: Principal component analysis in linear systems: controllability, observability, and model reduction. IEEE Trans. Autom. Control AC–26(1), 17–32 (1981). https://doi.org/10.1109/TAC.1981. 37. Penzl, T.: Eigenvalue decay bounds for solutions of Lyapunov equations: the symmetric case. Syst. Cont. Lett. 40, 139–144 (2000). https://doi.org/10.1016/S0167-6911(00)00010-4 38. Petersson, D.: A Nonlinear Optimization Approach to H2-Optimal Modeling and Control. Ph.D. thesis, Linkoping ¨ University. Available from http://www.diva-portal.org/smash/get/diva2:647068/ FULLTEXT01.pdf (2013) 39. Petersson, D., Lofber ¨ g, J.: Model reduction using a frequency-limited H -cost. Sys. Control Lett. 67, 32–39 (2014) 40. Redmann, M., Kurschner, ¨ P.: An H -Type Error Bound for Time-Limited Balanced Truncation. Tech. Rep. 1710.07572v1. ArXiv e-prints (2017) 41. Ruhe, A.: Rational Krylov sequence methods for eigenvalue computation. Linear Algebra Appl. 58, 391–405 (1984) 42. Ruhe, A.: The rational Krylov algorithm for nonsymmetric Eigenvalue problems. III: complex shifts for real matrices. BIT 34, 165–176 (1994) 43. Saad, Y.: Numerical Solution of Large Lyapunov Equation. In: Kaashoek, M.A., van Schuppen, J.H., Ran, A.C.M. (eds.) Signal Processing, Scattering, Operator Theory and Numerical Methods, pp. 503– 511, Birkhauser (1990) 44. Saad, Y.: Analysis of some Krylov subspace approximations to the matrix exponential operator. SIAM J. Numer. Anal. 29(1), 209–228 (1992) 45. Sabino, J.: Solution of large-scale Lyapunov equations via the block modified smith method. Ph.D. thesis, Rice University, Houston, Texas. http://www.caam.rice.edu/tech reports/2006/TR06-08.pdf (2007) 46. Simoncini, V.: Analysis of the rational Krylov subspace projection method for large-scale algebraic Riccati equations. SIAM J. Matrix Anal. Appl. 37(4), 1655–1674 (2016). https://doi.org/10.1137/16M 1059382 1844 P. Kurschner ¨ 47. Sinani, K., Gugercin, S.: Iterative Rational Krylov Algorithms for Unstable Dynamical Systems and Optimality Conditions for a Finite-Time Horizon (2017). http://meetings.siam.org/sess/dsp talk.cfm? p=81168 SIAM CSE (2017) 48. Truhar, N., Veselic, ´ K.: Bounds on the trace of a solution to the Lyapunov equation with a general sta- ble matrix. Syst. Cont. Lett. 56(7–8), 493–503 (2007). https://doi.org/10.1016/j.sysconle.2007.02.003 49. Vuillemin, P.: Frequency-limited model approximation of large-scale dynamical models. Ph.D. thesis, Universite de Toulouse. https://hal.archives-ouvertes.fr/tel-01092051 (2014) http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Advances in Computational Mathematics Springer Journals

Balanced truncation model order reduction in limited time intervals for large systems

Free
24 pages

Loading next page...
 
/lp/springer_journal/balanced-truncation-model-order-reduction-in-limited-time-intervals-0Sdr5Kna0I
Publisher
Springer Journals
Copyright
Copyright © 2018 by The Author(s)
Subject
Mathematics; Computational Mathematics and Numerical Analysis; Mathematical Modeling and Industrial Mathematics; Mathematical and Computational Biology; Computational Science and Engineering; Visualization
ISSN
1019-7168
eISSN
1572-9044
D.O.I.
10.1007/s10444-018-9608-6
Publisher site
See Article on Publisher Site

Abstract

Adv Comput Math (2018) 44:1821–1844 https://doi.org/10.1007/s10444-018-9608-6 Balanced truncation model order reduction in limited time intervals for large systems Patrick Kurschner ¨ Received: 7 July 2017 / Accepted: 16 April 2018 / Published online: 5 June 2018 © The Author(s) 2018 Abstract In this article we investigate model order reduction of large-scale systems using time-limited balanced truncation, which restricts the well known balanced trun- cation framework to prescribed finite time intervals. The main emphasis is on the efficient numerical realization of this model reduction approach in case of large sys- tem dimensions. We discuss numerical methods to deal with the resulting matrix exponential functions and Lyapunov equations which are solved for low-rank approx- imations. Our main tool for this purpose are rational Krylov subspace methods. We also discuss the eigenvalue decay and numerical rank of the solutions of the Lyapunov equations. These results, and also numerical experiments, will show that depending on the final time horizon, the numerical rank of the Lyapunov solutions in time- limited balanced truncation can be smaller compared to standard balanced truncation. In numerical experiments we test the approaches for computing low-rank factors of the involved Lyapunov solutions and illustrate that time-limited balanced truncation can generate reduced order models having a higher accuracy in the considered time region. Keywords Lyapunov equation · Rational Krylov subspaces · Model order reduction · Balanced truncation · Matrix exponential Mathematics Subject Classification 2010 15A16 · 15A18 · 15A24 · 65F60 · 93A15 · 93C Communicated by: Peter Benner Patrick Kurschner ¨ kuerschner@mpi-magdeburg.mpg.de Computational Methods in Systems and Control Theory, Max Planck Institute for Dynamics of Complex Technical Systems, Sandtorstr. 1, Magdeburg, 39106, Germany 1822 P. Kurschner ¨ 1 Introduction 1.1 Model reduction of linear systems Consider continuous-time, linear, time-invariant (LTI) systems x( ˙ t) = Ax(t ) + Bu(t ), x(0) = 0, (1a) y(t) = Cx(t) (1b) n×n n×m p×n with A ∈ R , B ∈ R ,and C ∈ R . Typically, the vector functions x(t) ∈ n m p R , u(t ) ∈ R ,and y(t) ∈ R are referred to as state, control, and, output vector, respectively. We assume that m, p  n and A is Hurwitz, i.e., (A) ⊂ C ,such that (1) is asymptotically stable. Given a large state space dimension n, model order reduction aims toward finding a reduced order model ˙ ˜ ˜ x( ˜ t) = Ax( ˜ t) + Bu(t ), x( ˜ 0) = 0, (2a) y( ˜ t) = Cx( ˜ t) (2b) r×r r×m p×r r with A ∈ R , B ∈ R , C ∈ R , x( ˜ t) ∈ R and a drastically reduced state dimension r  n. The smaller reduced system (2) should approximate the input- output behavior or the original system (1) well. Moreover, simulating the system, i.e., solving the differential equations in (2) for many different control functions u(t ) should be numerically much less expensive compared to the original system (1). For the approximation of (1) regarding the actual time domain behavior, it is desired that for all feasible input functions u(t ), y( ˜ t) ≈ y(t) for t ≥ 0, (3a) in other words, y(t)−˜ y(t) should be small ∀t ≥ 0 in some norm ·. With the help of the Laplace transformation, one can also formulate the approximation problem in the frequency domain, e.g., via H(ıω) ≈ H(ıω) for ω ∈ R,ı =−1, −1 −1 ˜ ˜ ˜ ˜ where H(s) = C(sI − A) B, H(s) = C(sI − A) B (3b) are the transfer function matrices of (1)and(2). There exist different model order reduction technologies and here we focus on balanced truncation (BT) [36] model order reduction. The backbone of BT are the infinite controllability and observabil- ity Gramians of (1) which are the the symmetric, positive semidefinite solutions P ,Q of the continuous-time, algebraic Lyapunov equations ∞ ∞ T T T T AP + P A =−BB ,A Q + Q A =−C C. (4) ∞ ∞ ∞ ∞ The Hankel singular values (HSV) of (1) are the square roots of the eigenvalues of the product P Q and are system invariants under state spac transformations. ∞ ∞ The magnitude of the HSV enables to identify components that are weakly con- trollable and observable. In BT this is achieved by first transforming (1)intoa balanced realization such that P = Q =  = diag (σ ,...,σ ). Dropping all ∞ ∞ ∞ 1 n states corresponding to small σ gives the reduced order model. Solving (4)forthe j Balanced truncation model order reduction in limited time... 1823 Algorithm 1: Square-root balanced truncation with low-rank factors Input : System matrices A, B, C defining an asymptotically stable dynamical system (1). ˜ ˜ ˜ Output: Matrices A, B, C of the reduced system. 1 Compute Z ,Z (e.g., with the methods described in [9, 46]), such that P Q T T Z Z ≈ P , Z Z ≈ Q in (4). P ∞ Q ∞ P Q 2 Compute thin singular value decomposition T T X X Y Y Z Z = XY = diag ( , ) P 1 2 1 2 1 2 with  = diag (σ ,...,σ ) containing the largest r (approximate) HSV. 1 1 r 1 1 − − 2 2 3 Construct T := Z Y  and S := Z X  . P 1 Q 1 1 1 4 Generate reduced order model T T ˜ ˜ ˜ A := S AT , B := S B, C := CT . (5) Gramians P ,Q is the computationally most challenging part of balanced trun- ∞ ∞ cation. For large-scale systems one therefore uses low-rank approximations of the T n×k Gramians instead, e.g., Z Z ≈ P with low-rank solution factors Z ∈ R , P ∞ P rank(Z ) = k  n, and likewise for Q . This strategy is backed up by the often P P ∞ numerically observed and theoretically explained rapid eigenvalue decay of solutions of Lyapunov equations [2, 3, 27, 37, 48]. The computation of the low-rank factors Z ,Z of P ,Q can be done efficiently by state-of-the-art numerical algorithms P Q ∞ ∞ for solving large Lyapunov equations [9, 46]. BT using low-rank factors Z ,Z of P Q the Gramians (4) is illustrated in Algorithm 1. The H -norm H  = sup(H(ıω) ). H 2 ω∈R of a stable system (1)isthe L induced norm of the convolution operator. It con- nects to the time-domain behavior via y ≤H  u . With exact Gramian L L 2 ∞ 2 T T factors, i.e., Z Z = P , Z Z = Q , BT is known to always generate a P ∞ Q ∞ P Q asymptotically stable reduced system for which the error bound H − H  ≤ 2 σ (6) H j j =r+1 holds. 1.2 Model reduction in limited time- and frequency intervals The approximation paradigms (3) enforce that the reduced system (2) is accurate for all times t ∈ R and frequencies ω ∈ R. From a practical point of view, achieving (3) might overshoot a realistic objective. For instance, if (1) models a mechanical or 1824 P. Kurschner ¨ electrical system, practitioners working with this model (and its approximation (2)) might only be interested in certain finite frequency intervals 0 ≤ ω <ω < ∞. 1 2 Likewise, when the final goal is to carry out simulations of (1), i.e., acquire time- domain solutions for various controls u(t ), one is usually only interested in y(t) for t smaller than some final time t < ∞. Hence, we consider time- and frequency restricted versions of (3) of the form y( ˜ t) ≈ y(t) for t ∈ T ⊂ R , (7a) H(iω) ≈ H(iω) for ω ∈  ⊂ R, =−, (7b) where the time- and frequency regions T , of interest should be provided by the underlying application. The expressions (7) demand that the reduced order model (2) is only accurate in T , but allow larger errors outside these regions. Compared to MOR approaches for the unrestricted setting (3), it is desired that MOR approaches for (7) achieve smaller approximation errors in T , with the same reduced order r. Alternatively, one demands that time- and frequency restricted MOR leads to com- parable approximation errors in T , with reduced systems of smaller order r.A secondary question is how the added time- and frequency restrictions influence the computational effort compared to an unrestricted MOR method of the same type. This issue will be in our particular focus. Typically, the time region in (7a)has the form T =[0,t ],t < ∞ (8a) e e which will also be the main situation in this work, but the more general restriction T =[t ,t ], 0 <t <t < ∞ (8b) s e s e will also be briefly discussed. Regarding the frequency restricted setting (7b), the typical regions are ˆ ˆ ˆ :=  ∪−,  := [ω ,ω ] with 0 ≤ ω < ... < ω <ω < ∞. i i+1 1 h h+1 i=1 Introducing time- or frequency restrictions into balanced truncation MOR has been originally proposed in [25] and further studied in, e.g., [8, 18, 21, 29]. In certain applications, e.g., circuit design, only single frequencies ω ∈ R might be of interest and an associated version of balanced truncation is addressed in [19]. H -MOR with limitations or weights on the frequency and time domain is investigated in, e.g., [11, 26, 32, 38, 39, 47, 49]. As continuation of our work in [8] regarding frequency-limited BT, we consider in this paper the numerically efficient realization of time-limited balanced truncation (TLBT) for large-scale systems. 1.3 Overview of this article We start in Section 2 by reviewing the general concept of time-limited BT from [25], mainly for restrictions of the form (8a). This includes the associated time-limited Gramians as well as the respective Lyapunov equations. Similar to standard BT, Balanced truncation model order reduction in limited time... 1825 executing TLBT for large systems heavily relies on how well the time-limited Grami- ans can be approximated by low-rank factorizations. This issue is investigated in Section 3.1 where we have a particular interest in the question when t induces sig- nificant differences between the infinite and time-limited Gramians. The actual issue of numerically dealing with the arising matrix functions and computing the low-rank factors of the time-limited Gramians is topic of Sections 3.2 and 3.3, respectively. Motivated by the promising results for frequency-limited BT studied in [8], we again employ rational Krylov subspace methods for this task. Section 4 collects different generalizations of TLBT including general state-space and certain differential alge- braic systems, the time restriction (8b), and stability preservation. In Section 5,the proposed numerical approach is tested numerically with respect to the approximation of the Gramians as well as the reduction of the dynamical system. 1.4 Notation The real and complex numbers are denoted by, respectively, R and C, R ,(R ), − + C (C ) are the sets of strictly negative (positive) real numbers and the open left − + (right) half plane. The space of real (complex) Matrices of dimension n × m is n×m n×m R (C ). For any complex quantity X = Re (X) + ı Im (X), we denote by Re (X) and Im (X) its real and, respectively, imaginary parts with ı being the imag- inary unit, and its the complex conjugate is X = Re (X) − ı Im (X). The absolute value of any real or complex scalar is denoted by |z|. Unless stated otherwise, · is the Euclidean vector- or subordinate matrix norm (spectral norm). By A and A = A we indicate the transpose and complex conjugate transpose of a real and n×n −1 complex matrix, respectively. If A ∈ C is a nonsingular, its inverse is A ,and −H H −1 −1 A = (A ) . Expressions of the form x = A b are always to be understood as solving a linear system Ax = b for x. The identity matrix of dimension n is indicated T m by I , and the vector of ones is denoted by 1 := (1,..., 1) ∈ R . The notation n m A 0(≺ 0) indicates symmetric positive (negative) definiteness of a symmetric or Hermitian matrix A,  () refers to semi-definiteness, and A  ()B means A − B  ()0. The spectrum of a matrix A is denoted by (A). 2 Gramians and balanced truncation for finite time horizons Since A is assumed to be Hurwitz, the infinite Gramians P ,Q in (4) can be ∞ ∞ represented in integral form as ∞ ∞ T T At T A t A t T At P = e BB e dt, Q = e C C e dt. (9) ∞ ∞ 0 0 Restricting the integration limits in the integrals (9) to a time interval [t ,t ] s e immediately yields the time-limited Gramians P , Q . T T 1826 P. Kurschner ¨ Definition 1 (Time-limited Gramians [25]) The time-limited reachability and observability Gramians of (1) with respect to the time-interval T =[t ,t ], 0 ≤ t < s e e t < ∞ are defined by t t e e At T A t T T At P = e BB e dt, Q = e C C e dt. (10) T T t t s s Theorem 1 (Lyapunov equations for the time-limited Gramians [25]) The time- limited Gramians P and Q according to Definition 1 are equivalently given in the T T following ways. a) The finite time Gramians P ,Q from (10) are given by T T T T At A t At A t s s e e P = e P e − e P e , (11a) T ∞ ∞ T T A t At A t At s s e e Q = e Q e − e Q e , (11b) ∞ ∞ where P and Q are the infinite reachability and observability Gramians (4). ∞ ∞ b) The time-limited Gramians P ,Q satisfy the time-limited reachability and T T observability Lyapunov equations T T T At At s e AP +P A =−B B + B B ,B := e B, B := e B (12a) t t t t T T s t e t s e s e T T T At At s e A Q +Q A =−C C + C C ,C := C e ,C := C e . (12b) t t t t T T t s t e s e s e Note that the time-limited Gramians (10) also exist for unstable A and if (A) ∩ (−A) =∅ they still solve the Lyapunov (12), see [40]. Hence, TLBT might be one possible approach to reduce unstable systems. Except for one brief numerical experi- ment, we will not pursue this topic any further. In analogy to the infinite time horizon case, the square roots of the eigenvalues of the product P Q are called time-limited T T Hankel singular values. A basic calculation shows that the time-limited Hankel sin- gular values are, as the standard Hankel singular values, invariant with respect to state-space transformations. By comparing the Lyapunov equations for the infinite Gramians (4) with (12), one immediately sees that the only differences are the inho- mogeneities, while the left hand sides are the same unchanged (adjoint) Lyapunov operators. This raises the question how much different the time-limited Gramians are from the infinite ones, and how this depends on the time interval of interest. We will pursue this in the next section, especially regarding the numerically important issue of how well the time-limited Gramians can be approximated by low-rank factoriza- T T tions P ≈ Z Z , Q ≈ Z Z . Of course, before approximately solving T P T Q T P T Q T T At the time-limited Lyapunov equations, the actions of the matrix exponentials e to B A t T (aswellase to C ) have to be dealt with numerically. Numerical approaches for handling the matrix exponential and computing low-rank solution factors Z ,Z P Q T T are topic of Section 3. A square-root version of TLBT is simply carried out by substituting Step 1 of Algorithm 1 by the code snippet shown in Algorithm 2 below and using the low-rank Balanced truncation model order reduction in limited time... 1827 solution factors Z , Z in the remaining steps. Depending on the time region of P Q T T interest, TLBT might in general not be a stability preserving method and, thus, there is also no H -error bound similar to (6). This can be regained by modifying TLBT further [29] and we discuss this issue in Section 4. Without this modification it is possible to establish an error bound in the H -norm [40]. Algorithm 2: Required changes in Algorithm 1 for TLBT At A t i i 1a Compute (approximation of) B := e B, C = C e , i ∈{s, e}. t t i i T T 1b Compute Z ,Z , such that P ≈ Z Z , Q ≈ Z Z in (12). P Q T P T Q T T T T P Q T T 3 Numerical computation of low-rank factors of time-limited Gramians This section is concerned with the actual numerical computation of low-rank factors of the time-limited Gramians. Before algorithms for dealing with the matrix exponen- tials and the Lyapunov equations are discussed, we briefly investigate the numerical ranks of the time-limited Gramians. For the sake of brevity, all considerations will be mostly restricted to the reachability Gramian because the observability Gramian can T T be dealt with similarly by replacing A, B with A ,C , respectively. Moreover, only the situation (8a) will be discussed, i.e., t = 0and 0 <t < ∞. s e 3.1 The difference of the infinite and time-limited Gramians In order to approximate P , Q by low-rank factorizations, it is desirable that their T T eigenvalues decay rapidly. For investigating this decay we assume from now that eigenvalues of symmetric (positive definite) matrices are given in a non-increasing order λ ≥ ... ≥ λ . The inhomogeneities of the time-limited Lyapunov (12)have 1 n up to twice the rank of their unlimited counterparts (4). Hence, by the available theory on the eigenvalue decay of solutions of matrix equations [2, 27, 37, 45] one expects that the eigenvalues of the time-limited Gramians decay somewhat slower than those of the infinite ones. We will see in the numerical experiments that, similar to the frequency-limited Gramians, in most situations it is the opposite case: the eigenvalues of the time-limited Gramians exhibit a faster decay and, consequently, have smaller numerical ranks. Assuming that (A, B) is controllable (rank[A − λI, B]= n, ∀λ ∈ C) it holds P 0 and, using (10), the relation (11a) yields 0  P = P − ∞ T ∞ T T At A t At A t e e e e e P e .Since E(t ) := e P e 0 it holds ∞ e ∞ P  P and also λ (P ) =P  P  = λ (P ). ∞ T 1 ∞ ∞ 2 T 2 1 T Due to the stability of A, the difference E(t ) = P − P will decay for increas- e ∞ ing values of t . Tight bounds for the eigenvalue behavior of P and E(t ) with e T e respect to the parameter t are difficult to derive and is a research topic of its own. Here, for simplicity we restrict the discussion to the case m = 1 and present a basic investigation of how the decay of E(t ) depends on t . e e 1828 P. Kurschner ¨ n −1 Lemma 1 Let B ∈ R , (A, B) controllable, A be diagonalizable, i.e., A = XX , −1 = diag (λ ,...,λ ), and define w := X B, X := X diag (w), N(t ) := 1 n B e t At A t H e e e X e diag (w).Then E(t ) = e P e = N(t )CN(t ) , e ∞ e e H t  t H e e P = P − E(t ) = P − N(t )CN(t ) = X C − e C e X , (13) ∞ e ∞ e e B −1 and C := is a Hermitian positive definite Cauchy matrix. λ +λ i j i,j =1 Proof Apply the eigendecomposition of A and P = X CX from [2, Lemma 3.2] ∞ B to (11a). Consider the impulse response of (1), At t y (t ) = Cη(t ), η(t ) := e B = X e w = N(t)1 , δ n indicating that the impulse-to-state-map η(t) and N(t) decay at a similar rate. With the spectral abscissa R := max Re (λ), the basic point wise bounds λ∈(A) Rt Rt η(t)≤ e Xw, N(t)≤ e Xw make this more visible. Using this and (13), we can conclude that for increasing t, E(t) is getting smaller at a similar speed as η(t). Consequently, a significant differ- ence between P and P might be only observed when t is chosen small enough T ∞ e in the sense that η(t ) has not reached an almost stationary state. The handling of the case m> 1 can be carried out similarly. Moreover, with a similar argumentation the decay rate of the time-limited Hankel singular values can be roughly connected to the decay of y (t ). To conclude, TLBT might be only practicable for small time horizons or for weakly damped systems. A similar investigation regarding TLBT for unstable systems is given in [47]. 3.2 Approximating the products with the matrix exponential There are several numerical approaches available to approximate the action of the matrix exponential to (a couple of) vectors, see, e.g., [1, 6, 12–15, 23, 24, 30, 31, 34, 35, 44]. Here, we are mostly interested in projection methods using (block) rational Krylov subspaces of the form ⎡ ⎤ −1 ⎣ ⎦ RK = span {q ,...,q },q = (A − s I) B (14) k 1 k k j j =1 where s ∈ C ∪ ıR ∪{∞} are called shifts. They represent the poles of a rational k + approximation r = ψ /φ of e with polynomials ψ ,φ of degree at most k k−1 k−1 k−1 k−1 k − 1. Balanced truncation model order reduction in limited time... 1829 n×km Let Q ∈ R have orthonormal columns and range (Q ) = RK .Thena k k k At Galerkin approximation [44]ofe B takes the form H t k e ˆ ˆ B ≈ B := Q B , B := e B , (15) t t ,k k t ,k t ,k k e e e e T T H := Q AQ ,B := Q B. k k k k k Note that for m = 1, the rational function r underlying the rational Krylov H t T zt k e approximation Q e Q B interpolates f(z) = e at (H ) (rational Ritz val- k k ues) [31, Theorem 3.3]. Further information on the approximation properties can be found in, e.g., [6, 10, 14–16, 30]. In the remainder we assume s =∞ s.t. q = B, 1 1 T T km×m m×m range B ⊆ range Q and B =[β , 0] ∈ R ,where β ∈ R .The ( ) ( ) k k orthonormal basis matrix Q and the restriction H can be efficiently computed by k k H t k e a (block) rational Arnoldi process [41]. Since H is low-dimensional, e can be computed by standard dense methods for the matrix exponential [34]. The choice of shifts s (k> 1) is crucial for a rapid convergence and several strategies exists for this purpose [31]. In this work we exclusively use adaptive shift generation tech- niques [14–16] because this appeared to be the safest strategy in the majority of experiments. For the situation m = 1 and after step k of the rational Arnoldi process, the next shift s is selected via k+1 s−z s = argmax |r (s)|,r (s) = , (16) k+1 k k s−s ∂S j =1 where z are the eigenvalues of H (Ritz values of A); s ,j = 1,...,k are the j k j previously used shifts; and ∂S marks a set of discrete points from the boundary of S approximating (A). We follow [16] and use S as the convex hull of (H ). k k k In the symmetric case, one can also use the spectral interval given by the largest and smallest eigenvalue of A [14, 15, 31]. For m> 1 we simply use each s in the denominator of r in (16) m times as in [16]. A different technique to deal with m> 1 includes, e.g., tangential rational Krylov methods [17]. The rational Krylov subspace method was chosen not only because of the good approximation properties, but also because the generated subspace is also a good candidate to acquire low-rank solution factors of (12). 3.3 Computing the low-rank Gramian factors Using (14), (15), a Galerkin approximation of the time-limited Gramian is P = T ,k Q Y Q ≈ P ,where Y solves the projected time-limited Lyapunov equation k k T k T T T ˆ ˆ H Y + Y H =−B B + B B , (17) k k k k t ,k k k t ,k H t ˆ k s (cf. [43]). It holds B = B = e B for the special case t = 0 considered here, k 0,k k s indicating already how t = 0 can be included as well. Hence, after B of sufficient s t ,k accuracy is found, we follow the idea in [8] by recycling the rational Krylov basis and continuing the rational Arnoldi process for (12a) until P leads to a sufficiently T ,k 1830 P. Kurschner ¨ Algorithm 3: Rational Krylov subspace method for time-limited Lyapunov (12a) Input : A, B, t as in (12a), tolerances 0 <τ ,τ  1. e f P n× T Output: Z ∈ R such that Z Z ≈ P with  ≤ mk  n. k k T 1 B = q β s.t. q q = I , Q = q . 1 1 m 1 1 2 for k = 1, 2,... do 3 Get new shift s via, e.g., (16). k+1 4 Solve (A − s I)g = q for g. k+1 k 5 Real, orthogonal expansion of Q : g = Re g − Q (Q Re g ). ( ) ( ) k + k 6 q = g β s.t. q q = I , Q =[Q ,q ]. k+1 + k k+1 m k+1 k k+1 k+1 7 if Im (s ) = 0 then k+1 8 k = k + 1, g = Im g − Q (Q Im g ) ( ) ( ) + k 9 q = g β s.t. q q = I , Q =[Q ,q ]. j +1 + k k+1 m k+1 k k+1 k+1 T T 10 H = Q AQ , B = Q B. k k k k k H t ˆ k e ˆ 11 B = e B , B = Q B t ,k k t ,k k t ,k e e e 12 if B − B /B  <τ then t ,k t ,k−1 t ,k f e e e T T T ˆ ˆ 13 Solve H Y + Y H + B B − B B = 0for Y . k k k k t ,k k k k t ,k T T T T T A(Q Y Q )+(Q Y Q )A +BB −B B k k k k t ,k k k t ,k 14 Set μ := . T T ˆ ˆ B B −B B k t ,k k t ,k 15 if μ <τ then k P T T 16 Y = SS , S S = I ,  = diag (γ ,...,γ ). k km 1 km 17 Truncate if necessary:  = diag (γ ,...,γ ), S = S(:, 1 : ), ≤ mk. 18 return low-rank solution factor Z = Q S . k k small norm of the Lyapunov residual. For the Lyapunov stage, the same adaptive shift generation (16) can be used [16]. The rational Krylov subspace method for generating a low-rank approximation Z Z ≈ P is illustrated in Algorithm 3. k T In the presence of a complex shift s , it is implicitly assumed that s is the k+1 k+1 subsequent shift. For this situation, the orthogonal expansion in Steps 5–9 of the already computed basis matrix Q by real basis vectors goes back to [42]. The pro- jected matrix H in Step 10, as well as the Lyapunov residual norm in Step 14 can be computed without accessing the large matrix A, see, e.g., [10, 16, 31, 41]. The small Lyapunov equation in Step 13 can be solved by standard methods for dense matrix equations, e.g., the Bartels-Stewart [4] method which we employ here. Once the scaled Lyapunov residual norm falls below the desired threshold, the rational Arnoldi process is terminated and the Steps 16–18 bring the computed low-rank Gramian approximation in the desired form Z Z and allow a rank truncation. Note that typi- cally, once the approximation of B is found, the generated subspace is already good enough to acquire a low-rank Gramian approximation without many additional iter- ations of the rational Arnoldi process. Often, the criteria in Steps 12–15 are satisfied in the same iteration step. Balanced truncation model order reduction in limited time... 1831 4 Extensions and further problems 4.1 Multiple time values Computing low-rank factors of the time-limited Gramians (12) for the more general approximation setting (8b) with a nonzero start time t , i.e., T =[t ,t ], can also be s s e done using Algorithm 3 with minor adjustments. Having computed a rational Krylov At basis for approximating e B for some time value t , the same basis typically also provides good approximations for any other time values [31]. Consequently, Algo- H t ˆ k s ˆ rithm 3 has to be simply changed by adding B = e B , B = Q B t ,k k t ,k k t ,k s s s and appropriately adjusting the steps regarding the Gramian approximation. Also the methods relying on Taylor approximations [1] can be easily modified to handle multi- ple time values. However, even if a nonzero t does not yield additional computational complications, this does not imply that TLBT will produce a accurate approximation of the transient behavior of (1)in [t ,t ]. Already the original TLBT paper [25] states s e that TLBT in [t ,t ] is expected to only give good approximations of the impulse s e m At response (u(t ) = vδ(t), v ∈ R ) y (t ) = C e Bv and numerical experiments confirm this. For an intuitive explanation assume for simplicity that no truncation is done in TLBT, i.e. in Step 2 of Algorithm 1. Then it holds range (Q ) = range (T ), At At At s e s e B ≈ B ∈ range (Q ),e B ≈ B ∈ range (Q ), and likewise for C e , k,t k k,t k s e At C e and range (S). Hence, the impulse response is accurately approximated at the relevant times. A(t −τ) For the response of an arbitrary input u(t ), y (t ) = C e Bu(τ )dτ,such argumentation clearly does not automatically hold. The value of y (t ) with respect to a general u(t ) depends, in general, on the values at the times before t. In the present form TLBT restricts the approximation to the time frame T and, thus, y (t ) will be poorly approximated for t ≤ t which, in turn, makes it difficult to acquire good approximations in T . One approach is to apply a time translation to the underlying system (1) such that the requested time-interval is transformed to [0,t − t ]. However, this time trans- e s lation will also introduce an inhomogeneous initial value x , which is an additional difficulty for model order reduction. Some strategies to cope with nonzero initial val- ues are given in [5, 33] and we plan to investigate the incorporation to time-limited BT in the future. 4.2 Generalized state-space and index-one descriptor systems Consider generalized state-space systems Mx( ˙ t) = Ax(t ) + Bu(t ), x(0) = 0, (18a) y(t) = Cx(t) (18b) n×n with M ∈ R nonsingular. Simple manipulations reveal that, similar to unre- stricted [7] and frequency-limited BT [8], the time-limited Gramians of (18)are P , T 1832 P. Kurschner ¨ M Q M,where P , Q solve the time-limited generalized Lyapunov equations T T T −1 T T T T M At −1 AP M + MP A =−B B + B B ,B := M e M B, (19a) T T t t t s t e t i s e −1 T T T T M At A Q M + M Q A =−C C + C C ,C := C e (19b) T T t t t t s t e i s e −1 AM t for t = t ,t . Note that B := e B. Low-rank factors of P , Q can be i s e t T T computed by methods for generalized Lyapunov equations. In particular, the ratio- nal Krylov subspace methods which we employ for approximating B will implicitly −1 −1 T −1 −T −1 work on M A, M B or alternatively, if M = LL 0, on L AL , L B. T T With low-rank approximations P ≈ Z Z , Q ≈ Z Z ,theSVDof T P T Q T P T Q T T Z MZ has to be used in Step 2 of Algorithm 1. In some of the numerical experiments we will encounter the situation M = A A B M 0 n ×n 1 1 2 1 f f C C ,A = ,B = ,C = [ ] with M ∈ R ,A ∈ 1 2 1 4 A A B 3 4 2 n−n ×n−n f f R nonsingular, s.t. (18) becomes a semi-explicit index-one descriptor system. Eliminating the algebraic constraints leads to a general state-space system −1 −1 −1 ˆ ˆ ˆ defined by M , A := A −A A A and B := B −A A B , C := C −C A A , 1 1 2 3 1 2 2 1 2 3 4 4 4 −1 and an additional feed-through term −C A B ,see [22]. Time-limited and unre- 2 2 ˆ ˆ ˆ stricted BT can be applied right away to this system via (19) defined by M , A, B, C. However, the matrix A will in general be dense and, thus, solving the linear sys- tems in Algorithm 3 can be very expensive. In the context of unrestricted BT, the authors of [22] exploit in a LR-ADI iteration for the Gramians that the arising ˆ ˆ ˆ dense linear systems (A − sM )V = W are equivalent to the sparse linear systems T T T T ˆ ˆ (A − sM)V = W with V =[V ,] , W =[W , 0] which are easier to solve numerically. We use the same trick within Step 4 of Algorithm 3. 4.3 Stability preservation and modified TLBT Because of the altered and sometimes indefinite right hand sides of (12), (19), TLBT is in general not stability preserving. Only when the used time interval is long enough such that the right hand sides are negative semi-definite, TLBT will produce an asymptotically stable reduced order model [25, Condition 1]. For the general sit- uation, a stability preserving modification of TLBT is proposed in [29]usingthe Lyapunov equations mod mod T T T AP M + MP A =−B B , mod T T mod (20) T mod T mod T A Q M + M Q A =−C C , mod T T mod 1 1 2 2 B B C C T B := Q diag |λ |,..., |λ | ,C := diag |λ |,..., |λ | Q . mod B mod 1 r 1 r C B C n×r n×r B C where Q ∈ R , Q ∈ R contain the eigenvectors corresponding to the B C r ≤ 2m, r ≤ 2p nonzero eigenvalues λ ,λ of the right hand sides of (12a), (19a) B C i i and, respectively, (12b), (19b). The rational Krylov subspace method in Algorithm 3 for M = I, t = 0 can be easily modified for (20) by replacing Step 13 with the steps shown in Algorithm 4. Balanced truncation model order reduction in limited time... 1833 Algorithm 4: Changes in Algorithm 3 for modified time-limited Gramians 13a Compute partial eigendecomposition T B B T B ββ 0 ˆ ˆ − B B Q = Q diag λ ,...,λ , Q Q = I , λ = 0. t ,k B B B r e r B t ,k 1 B i e B B B 2 13b Factor of projected modified rhs B := Q diag |λ |,..., |λ | . mod,k B 1 B T T 13c Solve H Y + Y H + B B for Y . k k k mod,k k k mod,k Generalization for M = I and 0 <t <t are straightforward. Note that because s e the modified time-limited Gramians do not fulfill a relation of the form (10)or (11), we cannot expect a fast singular value decay similar to P ,Q . As observed T T in [8, 29] for modified frequency-limited BT, modified TLBT might also lead to less accurate approximations in the considered time region compared to unmodified TLBT. 5 Numerical experiments Here, we illustrate numerically the results of Section 3 regarding the numerical rank of the frequency-limited Gramians as well as the numerical method for computing low-rank factors of P ,Q . Afterwards, the quality of the approximations obtained T T by TLBT with the low-rank factors is evaluated and compared against unlimited BT. Further topics like nonzero starting times and the modified TLBT scheme are also examined along the way. All experiments are done in MATLAB 8.0.0.783 ® ® (R2012b) on a Intel Xeon CPU X5650 @ 2.67GHz with 48 GB RAM. Table 1 list the used examples and their properties. For time-domain simulation of (1), (18)and the reduced order models for a given input function u(t ), an implicit midpoint rule with a fixed small time step t is used. Because the impulse response of (1)for m At u(t ) = δ(t)v, v ∈ R , x = 0 can be expressed as y (t ) = C e Bv, it is computed 0 δ by the same integrator applied to the uncontrolled system (u(t ) = 0) with initial con- dition x = Bv. For the impulse, or step responses, the vector distributing the control to the columns of B is set to v = 1 . Table 1 Dimensions, final integration time t , step size t, matrix properties, and source of the test systems Example nmpt t properties source bips 606 7135 4 4 20 0.04 index-1, n = 606, M = I morwiki ,[22] f 1 n bips 3078 21128 4 4 20 0.04 index-1, n = 3078, M = I morwiki ,[22] f 1 n vertstand 16626 6 6 600 0.6 A ≺ 0, M 0, C random morwiki ,[28] rail 79841 7 6 400 0.4 A ≺ 0, M 0 Oberwolfach Collection , ID=38881 https://morwiki.mpi-magdeburg.mpg.de/morwiki http://portal.uni-freiburg.de/imteksimulation/downloads/benchmark 1834 P. Kurschner ¨ Fig. 1 Scaled eigenvalues of P (top left), Hankel and time-limited Hankel singular values λ(P Q ) T T T (bottom left), norm of impulse response y (t ) (top right), and numerical ranks of infinite, (modified) time-limited Gramians against varying final times t (bottom right) for the bips 606 system Because the bips systems have eigenvalues very close to the imaginary axis which caused difficulties for all considered numerical methods, the shifted matrix A − 0.08M is used instead as in [22]. The generalized state-space systems vertstand and rail are dealt with Cholesky factorizations M = LL as explained in Section 4.2. There, the sparse Cholesky factors L are computed by the MATLAB command chol(M,’vector’). Matrix exponentials and Lya- punov equations defined by smaller dense matrices (including the projected ones in Algorithm 3) are solved directly by the expm and lyap routines. 5.1 Decay of the eigenvalues of the Gramians and the time-limited singular values At first we investigate how the end time t influences the eigenvalue decay of the time-limited Gramians. The index one descriptor system bips 606 is used for this experiment and reformulated into an equivalent state space system of dimension n = 606 as explained in Section 4.2. This comparatively small size allows a direct computation of the matrix exponentials and the Gramians. The top left plot in Fig. 1 shows the scaled and ordered eigenvalues λ /λ of the infinite reachability Gramian j 1 P and the time-limited one P for t = 1, 3, 10. Obviously, a distinctly faster ∞ T e eigenvalue decay of P is only observed for small time values t = 1, 3. As the final T e time increases, the eigenvalues move closer to the ones of P . The eigenvalues of ∞ Balanced truncation model order reduction in limited time... 1835 the observability Gramians exhibit a similar behavior. This observation is even more drastic for the decay of the time-limited Hankel singular values shown in the bottom left plot. For the largest value t = 10, hardly any difference to the Hankel singular values is visible. In the top right plot the point wise norm of the impulse response y (t ) shows that after t = 10, y (t ) has already almost reached its stationary phase. δ e δ This confirms the expectation that significant differences between infinite and time- limited Gramians occur only for times t which are small with respect to the behavior of the impulse response. The bottom right plots shows the numerical rank of infinite, time-limited and modified time-limited Gramians against t . The numerical ranks of P ,Q clearly move toward the numerical ranks of P ,Q as t increases. In T T ∞ ∞ e mod mod contrast, the numerical ranks of P ,Q (Section 4.3) are always close to the T T ones of P ,Q and appear to be largely unaffected by different values of t .This ∞ ∞ e is a very similar behavior as observed for the modified frequency-limited Gramians in [8]. 5.2 Computing low-rank factors of the time-limited Gramians We proceed by testing the computation of low-rank factors of the infinite and (mod- ified) time-limited reachability Gramians by the rational Krylov subspace method in Algorithm 3. The stopping criteria for matrix function and Gramian approximations −8 use the thresholds τ = τ = 10 . To save some computational cost, the projected f P matrix exponentials and Lyapunov equations (Steps 11 and 13 in Algorithm 3) are only dealt with every 5th iteration step. Table 2 summarizes the used time values t , the produced subspace dimensions d, the ranks r of the low-rank approximations, and the total computing times t for rk mod the Gramians P , P ,and P . Apparently, larger rational Krylov subspaces and ∞ T approximately twice as long computation times are needed to obtain low-rank factors of the time-limited Gramians, but the final ranks of the low-rank approximations of the time-limited Gramians P are in all cases smaller compared to P . The modified T ∞ mod time-limited Gramians P do not show this behavior as they have ranks similar to the infinite Gramians. For the rail example, increasing the time horizon t does not change the required subspace dimension d for the approximation of the time-limited Table 2 Results of the numerical numerical computation of low-rank factors of the different Gramians: time horizon t , generated subspace dimension d,rank r of the low-rank approximations, and computing time t in seconds rk mod P P P ∞ T Example t dr t dr t dr t e rk rk rk bips 3078 3 308 245 18.33 664 131 105.89 664 241 106.17 vertstand 300 180 164 6.18 300 140 13.42 300 183 13.55 rail 10 245 224 34.71 420 168 74.06 420 228 76.06 100 245 224 34.71 420 198 76.53 420 227 74.68 1836 P. Kurschner ¨ Gramian P , but its rank r is clearly larger. In all cases, no additional steps of the At rational Krylov method were necessary once the approximation of e B was found. The results for the observability Gramians were largely similar. For the bips 3078 example, the observability Gramians appeared to be more demanding for Algorithm 3 than the reachability Gramians. 5.3 Model reduction results Now we execute BT [36]aswellas(modified)TLBT[25, 29] using the square-root method (Algorithm 1) with the low-rank factors of the reachability and observability Gramians computed in the previous section. For this purpose, we restrict ourselves to the reduction to fixed specified orders r. The approximation quality of the con- structed reduced order models is assessed via the point wise and maximal relative error norms y(t)−y (t ) r 2 E (t ) := ,t ≤ t , E := max E (t ) f T y(t) t ∈[0,t ] of the output responses y(t), y( ˜ t) of original and reduced order models. We consider the response y (t ) for the impulse input u(t ) = δ(t)1 for all examples. Moreover, δ m for each used test system, the transient response with respect to an additional input signal u(t ) is also considered. We use step like input signals u(t ) = 1 and u(t ) = 501 for the bips 3078 and, respectively, rail example, and u := [5 · 10 · m ∗ 2 T 0.198(sin(t π/100) ), 4, 2, 1, 3, 1] for vertstand. The results are given in Table 3 listing the largest relative error E in [0,t ] and the overall computation time t , i.e., the computation time for computing low-rank mor factors of both reachability and observability Gramians by Algorithm 3 plus the time for Algorithm 1 to execute the BT variants. We also indicate whether the produced reduced order models are asymptotically stable (s = 1) or unstable (s = 0).For some selected settings of u(t ) and t , the system responses and point wise relative Table 3 Model reduction results by BT, TLBT, and modified TLBT using different t and u(t ): reduced order r, largest relative error E in [0,t ], overall computation time t ,and s ∈{0, 1} indicates if T e mor reduced system is asymptotically stable or unstable BT TLBT mod. TLBT Example ut E t s E t s E t s e T mor T mor T mor bips 3078, δ 3 5.10e-04 291.3 1 1.08e-06 331.9 0 5.15e-04 399.6 1 r = 100 1 3 6.90e-06 291.3 1 6.33e-09 331.9 0 1.20e-05 399.6 1 vertstand, δ 300 8.90e-02 13.2 1 2.44e-02 30.4 1 6.48e-02 31.6 1 r = 20 u 300 8.26e-03 13.2 1 9.18e-04 30.4 1 9.95e-03 31.6 1 rail, δ 10 6.60e-01 62.2 1 5.66e-04 130.8 0 6.59e-01 137.4 1 r = 50 δ 100 6.60e-01 62.2 1 7.90e-03 135.7 1 6.60e-01 136.4 1 50 10 1.78e-03 62.2 1 6.08e-07 130.8 0 2.61e-03 137.4 1 50 100 1.78e-03 62.2 1 2.57e-05 135.7 1 2.15e-03 136.4 1 Balanced truncation model order reduction in limited time... 1837 Fig. 2 Responses of original and reduced systems obtained by different BT versions: (top left) bips 3078, t = 3, u(t ) =impulse, r = 100, (top right) vertstand, t = 300, u(t ) = u , r = 20, e e ∗ (bottom left) rail, t = 10, u(t ) =impulse, r = 50, (bottom left) rail, t = 100, u(t ) = 50, r = 50 e e errors E (t ) are plotted against the time t in Figs. 2 and 3, respectively. Figure 4 shows the behavior of E as the reduced order r increases. Fig. 3 Relative errors E (t ) obtained by different BT versions: (top left) bips 3078, t = 3, u(t ) =impulse, r = 100, (top right) vertstand, t = 300, u(t ) = u , r = 20, (bottom left) rail, e ∗ t = 10, u(t ) =impulse, r = 50, (bottom left) rail, t = 100, u(t ) = 50, r = 50 e e 1838 P. Kurschner ¨ Fig. 4 Maximum relative errors E in [0,t ] against increasing reduced orders r for different u(t ) and T e BT variants: (top left) bips 3078, t = 3, u(t ) =impulse, (top right) vertstand, t = 300, u(t ) = u , e e ∗ (bottom left) rail79k, t = 10, u(t ) =impulse, (bottom left) rail79k, t = 100, u(t ) = 50 e e Apparently, for the chosen orders r and in the time regions of interest, the largest relative errors E of the reduced order models returned by TLBT are in most experi- ments more than one order of magnitude smaller compared to standard and modified time-limited BT. The plots in Fig. 4 also indicate that much larger reduced order models are needed for BT and modified TLBT to achieved the same accuracy as unmodified TLBT. Figure 3 shows that after leaving the time region T , TLBT deliv- ers larger errors. However, Table 3 also reveals that executing TLBT and its modified version is more time consuming than standard BT. This is a direct consequence of the higher computation times for getting the low-rank factors of the (modified) time-limited Gramians which was pointed out in the previous subsection (Table 2). Using the concept of angles between subspaces or the modal assurance crite- T 2 2 2 rion (MAC) (MAC(x, y) =|y x| /(x y ),[20]) indicated that the spaces spanned by the projection matrices T , S in TLBT and T , S in unre- TLBT TLBT BT BT stricted BT, respectively, are different. For example, computing the MAC for the right n ×100 projection matrices T , T ∈ R for the bips 3078 system showed that TLBT BT only a few of the columns of both matrices are well correlated to each other (i.e., MAC(T e ,T e ) ≈ 1 for very few i, j ∈{1,..., 100}). TLBT i BT j Moreover, albeit the higher accuracy in [0,t ], in some cases TLBT produces unstable reduced order models. This is especially visible in the upper left plot of Fig. 2 showing the impulse responses of bips 3078. For times t ≥ t outside the interval [0,t ], the impulse response of reduced order model generated by TLBT exhibits an exponential growth and departs from the original response. As illustrated Balanced truncation model order reduction in limited time... 1839 Table 4 Results of BT, TLBT, and modified TLBT reduction to r = 15 of vertstand example with respect to time frame T =[t ,t ]=[50, 100] s e BT TLBT mod. TLBT input u E t s E t s E t s T mor T mor T mor δ 7.95e-04 14.2 1 3.07e-04 37.6 0 9.21e-04 32.5 1 u 6.34e-03 14.2 1 1.78e-02 37.6 0 6.85e-03 32.5 1 with the rail examples and proven in [25], using a higher end time t can already cure this. Modified TLBT does not generate unstable reduced systems, but its approx- imation quality in [0,t ] is very close to standard BT without time restrictions. Taking also into account the higher computational costs of modified TLBT given in Table 3, the introduction of the time restriction is rendered essentially redundant because no smaller errors are achieved in the targeted time region. Hence, if stability preservation in the reduced order model is crucial, we recommend to stick to standard BT. To conclude, TLBT fulfills the goal to acquire smaller errors in the desired time interval [0,t ], but at the price of somewhat larger execution times because the computation of required low-rank Gramian factors is currently more costly. Now we carry out one experiment to evaluate the approximation qualities of TLBT for nonzero t .The vertstand example is used and the reduced order model should approximate the output y(t) with respect to u(t ) = δ(t)v and u(t ) = u =[5 · 10 · 2 T 0.198(sin(t π/100) ), 4, 2, 1, 3, 1] in the time window T =[t ,t ]=[50, 100]. s e The low-rank factors of the Gramians are computed as before using Algorithm 3 with the required small extensions mentioned in Section 4.1. The results are summarized in Table 4 and Fig. 5 illustrates the relative errors plots as well as the largest relative error in T against the reduced order r. Compared to stan- dard BT, TLBT delivers, as expected, more accurate results of the impulse response, but fails for u(t ) = u . Moreover, the bottom left plot Fig. 5 shows that the decay of E with respect to the impulse response is less monotonic as in the case t = 0s.t. T s for some reduced orders r, especially larger ones, standard BT outperforms TLBT. The failure of TLBT to approximate the transient response for arbitrary inputs u(t ) was also observed in other experiments with different time intervals and examples, and occasionally even for the impulse response. In experiments with smaller systems, using exact matrix exponentials and Gramian factors did not yield any improvement of the reduction results such that the problems with TLBT are most likely not a result of inaccurate Gramian approximations. In summary, although TLBT in the current form does indeed deliver significantly more accurate reduced order models in the time intervals [0,t ], the same cannot be said when time intervals [t ,t ], t > 0are e s e s considered. Improving the performance of TLBT in this scenario is therefore subject of further research. As the final experiment we briefly test the application of TLBT for the reduction of unstable systems. For this purpose a modification of the bips 606 example is ˆ ˆ used with A := A + 0.2M leading to max Re (λ) = 0.323, λ ∈ (A). The final 1840 P. Kurschner ¨ Fig. 5 Results for the reduction of the vertstand example in [t ,t ]=[50, 100]. The top plots show s e the relative error norms against t for u(t ) = δ(t)v (left) and u(t ) = u (right) after a reduction to r = 15. The behavior of E in [t ,t ] for increasing reduced dimensions r is illustrated in the bottom plots T s e time is set to t = 5 and a reduced order model with r = 100 is generated. Due to the comparatively small system dimension (n = 606), the matrix exponentials and Gramians could be computed by direct, dense methods for these tasks. Moreover, the employed rational Krylov methods from the experiments before were not able to compute the required low-rank Gramians factors. The impulse response of exact and reduced system and the relative error are illustrated in Fig. 6. Apparently, TLBT was able to reduce this mildly unstable system to a reduced −2 order model with maximal relative error E ≈ 1.2 · 10 in [0,t ]. Hence, provided T e t is chosen in a reasonable way, e.g., with a sufficient distance to the exponential Fig. 6 Results for the reduction of the unstable variation of the bips 606 example for t = 5, r = 100, u = δ(t)v (left: system response, right: relative error) Balanced truncation model order reduction in limited time... 1841 growth of y(t), TLBT appears to be a potential candidate from model order reduc- tion of unstable systems. However, the applications to large-scale unstable systems is currently still difficult because algorithms for computing low-rank solutions of Lya- punov equations usually require that A is (anti)stable. Advances in this direction are, therefore, necessary to pursue this type of reduction further. 6 Conclusions BT model order reduction for large-scale systems restricted to finite time inter- vals [25] was investigated. The resulting Lyapunov equations that have to be solved numerically also include the matrix exponential in their inhomogeneities. We first showed that the difference of time-limited and infinite Gramians decays for increas- ing times in a similar way as the impulse response of the underlying system. Hence, for small time intervals, a reduced numerical rank of the time-limited Gramians can be observed. Future research should further investigate the influence of the chosen end time on the eigenvalue decay of the Gramians. As in frequency-limited BT [8], we proposed to handle the matrix exponen- tials and the Lyapunov equations by an efficient rational Krylov subspace method incorporating a subspace recycling idea. While this numerical approach already led to satisfactory results, there is some room for improvement, especially regarding the approximation of the action of the matrix exponential. In this context, further work could include, for instance, enhanced strategies like tangential directions [17], different inner products [23], or using different methods [1, 12] altogether. The numerical reduction experiments indicated that TLBT is able to acquire sev- eral orders of magnitude more accurate reduced order models in time intervals of the form [0,t ] at a somewhat higher, although still comparable, numerical effort. Similar techniques were applied for stability preserving modified TLBT [29]which, however, could not keep up with standard or time-limited BT in terms of efficiency and accuracy of the reduced order models. Hence, we recommend to use standard BT if the preservation of stability is an irrevocable goal. For keeping the high accuracy in the time-interval of interest, a different way to make TLBT a stability preserving method has to be found. TLBT for time regions [t ,t ] with nonzero start times t > 0 s e s provided much worse results than with t = 0, except for approximating the impulse response. In the form introduced in [25], TLBT appears to be incapable of produc- ing good reduced order models with respect to [t ,t ] and, thus, further research is s e necessary in this direction. The results in [5, 33] might be one possible ingredient for this. A brief experiment also indicated that TLBT can be employed to reduce unsta- ble systems. The efficient computation of low-rank Gramian factors in this case is currently not as advanced as in the stable situation, making this a further interesting research topic. Acknowledgments Open access funding provided by Max Planck Society. I thank the referees for their helpful comments. Moreover, I am grateful for the constructive discussions with Maria Cruz Varona, Serkan Gugercin, and Stefan Guettel. 1842 P. Kurschner ¨ Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, dis- tribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. References 1. Al-Mohy, A.H., Higham, N.J.: Computing the action of the matrix exponential, with an application to exponential integrators. SIAM J. Sci. Comput. 33(2), 488–511 (2011) 2. Antoulas, A.C., Sorensen, D.C., Zhou, Y.: On the decay rate of Hankel singular values and related issues. Syst. Cont. Lett. 46(5), 323–342 (2002) 3. Baker, J., Embree, M., Sabino, J.: Fast singular value decay for Lyapunov solutions with nonnormal coefficients. SIAM J. Matrix Anal. Appl. 36(2), 656–668 (2015) 4. Bartels, R.H., Stewart, G.W.: Solution of the matrix equation AX + XB = C: Algorithm 432. Comm. ACM 15, 820–826 (1972) 5. Beattie, C., Gugercin, S., Mehrmann, V.: Model reduction for systems with inhomogeneous initial conditions. Syst. Control Lett. 99, 99–106 (2017). https://doi.org/10.1016/j.sysconle.2016.11.007 6. Beckermann, B., Reichel, L.: Error estimates and evaluation of matrix functions via the faber transform. SIAM J. Numer. Anal. 47(5), 3849–3883 (2009). https://doi.org/10.1137/080741744 7. Benner, P.: Solving large-scale control problems. IEEE Control Syst. Mag. 14(1), 44–59 (2004) 8. Benner, P., Kurschner, ¨ P., Saak, J.: Frequency-limited balanced truncation with low-rank approxima- tions. SIAM J. Sci. Comput. 38(1), A471–A499 (2016). https://doi.org/10.1137/15M1030911 9. Benner, P., Saak, J.: Numerical solution of large and sparse continuous time algebraic matrix Ric- cati and Lyapunov equations: a state of the art survey. GAMM Mitteilungen 36(1), 32–52 (2013). https://doi.org/10.1002/gamm.201310003 10. Berljafa, M., Guttel, ¨ S.: Generalized rational krylov decompositions with an application to rational approximation. SIAM J. Matrix Anal. Appl. 36(2), 894–916 (2015) 11. Breiten, T., Beattie, C., Gugercin, S.: Near-optimal frequency-weighted interpolatory model reduc- tion. Sys. Control Lett. 78, 8–18 (2015) 12. Caliari, M., Kandolf, P., Ostermann, A., Rainer, S.: Comparison of software for computing the action of the matrix exponential. BIT 54(1), 113–128 (2014) 13. Davies, P.I., Higham, N.J.: Computing f(A)b for matrix functions f . In: Boric¸i, A., Frommer, A., Joo, ´ B., Kennedy, A., Pendleton, B. (eds.) QCD and Numerical Analysis III, Lect. Notes Comput. Sci. Eng., vol. 47, pp. 15–24. Springer, Berlin (2005). https://doi.org/10.1007/3-540-28504-0 2 14. Druskin, V., Knizhnerman, L., Zaslavsky, M.: Solution of large scale evolutionary problems using rational krylov subspaces with optimized shifts. SIAM J. Sci. Comput. 31(5), 3760–3780 (2009) 15. Druskin, V., Lieberman, C., Zaslavsky, M.: On adaptive choice of shifts in rational krylov subspace reduction of evolutionary problems. SIAM J. Sci. Comput. 32(5), 2485–2496 (2010) 16. Druskin, V., Simoncini, V.: Adaptive rational Krylov subspaces for large-scale dynamical systems. Syst. Control Lett. 60(8), 546–560 (2011) 17. Druskin, V., Simoncini, V., Zaslavsky, M.: Adaptive tangential interpolation in rational Krylov subspaces for MIMO dynamical systems. SIAM J. Matrix Anal. Appl. 35(2), 476–498 (2014). https://doi.org/10.1137/120898784 18. Du, X., Benner, P.: Finite-Frequency Model Order Reduction of Linear Systems via Parameterized Frequency-dependent Balanced Truncation. Tech. Rep. 1602.04408. ArXiv e-prints (2016) 19. Du, X., Benner, P., Yang, G., Ye, D.: Balanced truncation of linear time-invariant systems at a sin- gle frequency. Preprint MPIMD/13-02, Max Planck Institute Magdeburg. Available from http://www. mpi-magdeburg.mpg.de/preprints/ (2013) 20. Ewins, D.: Modal testing: theory, practice, and application (2nd edn) Research Study Press LTD (2000) 21. Fehr, J., Fischer, M., Haasdonk, B., Eberhard, P.: Greedy-based approximation of frequency-weighted Gramian matrices for model reduction in multibody dynamics. Z. Angew. Math. Mech. 93(8), 501– 519 (2013) 22. Freitas, F., Rommes, J., Martins, N.: Gramian-based reduction method applied to large sparse power system descriptor models. IEEE Trans. Power Syst. 23(3), 1258–1270 (2008) Balanced truncation model order reduction in limited time... 1843 23. Frommer, A., Lund, K., Szyld, D.B.: Block Krylov subspace methods for computing functions of matrices applied to multiple vectors. Tech. Rep. 17-03-21, Department of Mathematics Temple University (2017) 24. Frommer, A., Simoncini, V.: Matrix functions. In: Schilders, W., van der Vorst, H., Rommes, J. (eds.) Model Order Reduction: Theory, Research Aspects and Applications, Mathematics in Industry, vol. 13, pp. 275–303. Springer, Berlin (2008). https://doi.org/10.1007/978-3-540-78841-6 13 25. Gawronski, W., Juang, J.: Model reduction in limited time and frequency intervals. Int. J. Syst. Sci. 21(2), 349–376 (1990). https://doi.org/10.1080/00207729008910366 26. Goyal, P., Redmann, M.: Towards time-limited H -optimal model order reduction. Tech. Rep. WIAS Preprint No. 2441 (2017) 27. Grasedyck, L.: Existence of a low rank or H -matrix approximant to the solution of a Sylvester equation. Numer. Lin. Alg. Appl. 11, 371–389 (2004) 28. Großmann, K., Stadel, C., Galant, A., Muhl, ¨ A.: Berechnung von Temperaturfeldern an Werkzeug- maschinen. Zeitschrift fu,r ¨ Wirtschaftlichen Fabrikbetrieb 107(6), 452–456 (2012) 29. Gugercin, S., Antoulas, A.C.: A survey of model reduction by balanced truncation and some new results. Internat. J. Control 77(8), 748–766 (2004). https://doi.org/10.1080/00207170410001713448 30. Guttel, ¨ S.: Rational Krylov Methods for Operator Functions. Ph.D. thesis, Technische Uni- versitat ¨ Bergakademie Freiberg, Germany. Available from http://nbn-resolving.de/urn:nbn:de:bsz: 105-qucosa-27645 (2010) 31. Guttel, ¨ S.: Rational Krylov approximation of matrix functions: Numerical methods and optimal pole selection. GAMM-Mitteilungen 36(1), 8–31 (2013). https://doi.org/10.1002/gamm.201310002 32. Halevi, Y.: Frequency weighted model reduction via optimal projection. IEEE Trans. Automat. Control 37(10), 1537–1542 (1992) 33. Heinkenschloss, M., Reis, T., Antoulas, A.C.: Balanced truncation model reduction for systems with inhomogeneous initial conditions. Automatica 47(3), 559–564 (2011). https://doi.org/10.1016/j.auto matica.2010.12.002 34. Higham, N.: Functions of matrices: Theory and Computation. Society for Industrial and Applied Mathematics, Philadelphia (2008) 35. Knizhnerman, L.A.: Calculation of functions of unsymmetric matrices using Arnoldi’s method. Comput. Math. Math. Phys. 31(1), 1–9 (1992) 36. Moore, B.C.: Principal component analysis in linear systems: controllability, observability, and model reduction. IEEE Trans. Autom. Control AC–26(1), 17–32 (1981). https://doi.org/10.1109/TAC.1981. 37. Penzl, T.: Eigenvalue decay bounds for solutions of Lyapunov equations: the symmetric case. Syst. Cont. Lett. 40, 139–144 (2000). https://doi.org/10.1016/S0167-6911(00)00010-4 38. Petersson, D.: A Nonlinear Optimization Approach to H2-Optimal Modeling and Control. Ph.D. thesis, Linkoping ¨ University. Available from http://www.diva-portal.org/smash/get/diva2:647068/ FULLTEXT01.pdf (2013) 39. Petersson, D., Lofber ¨ g, J.: Model reduction using a frequency-limited H -cost. Sys. Control Lett. 67, 32–39 (2014) 40. Redmann, M., Kurschner, ¨ P.: An H -Type Error Bound for Time-Limited Balanced Truncation. Tech. Rep. 1710.07572v1. ArXiv e-prints (2017) 41. Ruhe, A.: Rational Krylov sequence methods for eigenvalue computation. Linear Algebra Appl. 58, 391–405 (1984) 42. Ruhe, A.: The rational Krylov algorithm for nonsymmetric Eigenvalue problems. III: complex shifts for real matrices. BIT 34, 165–176 (1994) 43. Saad, Y.: Numerical Solution of Large Lyapunov Equation. In: Kaashoek, M.A., van Schuppen, J.H., Ran, A.C.M. (eds.) Signal Processing, Scattering, Operator Theory and Numerical Methods, pp. 503– 511, Birkhauser (1990) 44. Saad, Y.: Analysis of some Krylov subspace approximations to the matrix exponential operator. SIAM J. Numer. Anal. 29(1), 209–228 (1992) 45. Sabino, J.: Solution of large-scale Lyapunov equations via the block modified smith method. Ph.D. thesis, Rice University, Houston, Texas. http://www.caam.rice.edu/tech reports/2006/TR06-08.pdf (2007) 46. Simoncini, V.: Analysis of the rational Krylov subspace projection method for large-scale algebraic Riccati equations. SIAM J. Matrix Anal. Appl. 37(4), 1655–1674 (2016). https://doi.org/10.1137/16M 1059382 1844 P. Kurschner ¨ 47. Sinani, K., Gugercin, S.: Iterative Rational Krylov Algorithms for Unstable Dynamical Systems and Optimality Conditions for a Finite-Time Horizon (2017). http://meetings.siam.org/sess/dsp talk.cfm? p=81168 SIAM CSE (2017) 48. Truhar, N., Veselic, ´ K.: Bounds on the trace of a solution to the Lyapunov equation with a general sta- ble matrix. Syst. Cont. Lett. 56(7–8), 493–503 (2007). https://doi.org/10.1016/j.sysconle.2007.02.003 49. Vuillemin, P.: Frequency-limited model approximation of large-scale dynamical models. Ph.D. thesis, Universite de Toulouse. https://hal.archives-ouvertes.fr/tel-01092051 (2014)

Journal

Advances in Computational MathematicsSpringer Journals

Published: Jun 5, 2018

References

You’re reading a free preview. Subscribe to read the entire article.


DeepDyve is your
personal research library

It’s your single place to instantly
discover and read the research
that matters to you.

Enjoy affordable access to
over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month

Explore the DeepDyve Library

Search

Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly

Organize

Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place.

Access

Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals.

Your journals are on DeepDyve

Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more.

All the latest content is available, no embargo periods.

See the journals in your area

DeepDyve

Freelancer

DeepDyve

Pro

Price

FREE

$49/month
$360/year

Save searches from
Google Scholar,
PubMed

Create lists to
organize your research

Export lists, citations

Read DeepDyve articles

Abstract access only

Unlimited access to over
18 million full-text articles

Print

20 pages / month

PDF Discount

20% off