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

In this article we investigate model order reduction of large-scale systems using time-limited balanced truncation, which restricts the well known balanced truncation framework to prescribed finite time intervals. The main emphasis is on the efficient numerical realization of this model reduction approach in case of large system dimensions. We discuss numerical methods to deal with the resulting matrix exponential functions and Lyapunov equations which are solved for low-rank approximations. 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.


Model reduction of linear systems
with A ∈ R n×n , B ∈ R n×m , and C ∈ R p×n . Typically, the vector functions x(t) ∈ R n , u(t) ∈ R m , and y(t) ∈ R p 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 towards finding a reduced order model with A ∈ R r×r , B ∈ R r×m , C ∈ R p×r ,x(t) ∈ R 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), 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Ã,B,C of the reduced system. 1 Compute Z P , Z Q (e.g., with the methods described in [9,46]), such that Z P Z T P ≈ P ∞ , Z Q Z T Q ≈ Q ∞ in (4). 2 Compute thin singular value decomposition with Σ 1 = diag (σ 1 , . . . , σ r ) containing the largest r (approximate) HSV. 3 Construct T := Z P Y 1 Σ in other words, y(t) −ỹ(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., viaH where H(s) = C(sI − A) −1 B,H(s) =C(sI r −Ã) −1B 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 observability Gramians of (1) which are the the symmetric, positive semidefinite solutions P ∞ , Q ∞ of the continuous-time, algebraic Lyapunov equations The Hankel singular values (HSV) of (1) are the eigenvalues of the product P ∞ Q ∞ and are system invariants under state space transformations. The magnitude of the HSV enables to identify components that are weakly controllable and observable. In BT this is achieved by first transforming (1) into a balanced realization such that P ∞ = Q ∞ = Σ ∞ = diag (σ 1 , . . . , σ n ). Dropping all states corresponding to small σ j gives the reduced order model. Solving (4) for the Gramians P ∞ , Q ∞ is the computationally most challenging part of balanced truncation. For large-scale systems one therefore uses low-rank approximations of the Gramians instead, e.g., Z P Z T P ≈ P ∞ with low-rank solution factors Z P ∈ R n×k P , rank(Z P ) = k P n, and likewise for Q ∞ . This strategy is backed up by the often 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 P , Z Q of P ∞ , Q ∞ can be done efficiently by state-of-the-art numerical algorithms for solving large Lyapunov equations [9,46]. BT using low-rank factors Z P , Z Q of the Gramians (4) is illustrated in Algorithm 1. The H ∞ -norm H H∞ = sup ω∈R ( H(ıω) 2 ). of a stable system (1) is the L 2 induced norm of the convolution operator. It connects to the time-domain behavior via y L2 ≤ H H∞ u L2 . With exact Gramian factors, i.e., Z P Z T P = P ∞ , Z Q Z T Q = Q ∞ , BT is known to always generate a asymptotically stable reduced system for which the error bound holds.

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 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 e < ∞. Hence, we consider time-and frequency restricted versions of (3) of the formỹ 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 comparable 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 which will also be the main situation in this work, but the more general restriction will also be briefly discussed. Regarding the frequency restricted setting (7b), the typical regions are 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 2 -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.

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, executing TLBT for large systems heavily relies on how well the time-limited Gramians 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 e induces significant 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 algebraic 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.

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 R n×m (C n×m ). 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 imaginary 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 T and A H = A T we indicate the transpose and complex conjugate transpose of a real and complex matrix, respectively. If A ∈ C n×n is a nonsingular, its inverse is A

Gramians and Balanced Truncation for Finite Time Horizons
Since A is assumed to be Hurwitz, the infinite Gramians P ∞ , Q ∞ (4) can be represented in integral form as Restricting the integration limits in the integrals (9) to a time interval [t s , t e ] immediately yields the time-limited Gramians P T , Q T .
Definition 2.1 (Time-limited Gramians [25]). The time-limited reachability and observability Gramians of (1) with respect to the time-interval T = [t s , t e ], 0 ≤ t e < t e < ∞ are defined by Theorem 2.1 (Lyapunov equations for the time-limited Gramians [25]). The time-limited Gramians P T and Q T according to Definition 2.1 are equivalently given in the following ways.
a) The finite time Gramians P T , Q T from (10) are given by where P ∞ and Q ∞ are the infinite reachability and observability Gramians (4).
b) The time-limited Gramians P T , Q T satisfy the time-limited reachability and observability Lyapunov equations Note that the time-limited Gramians (10) also exist for unstable A and if Λ(A) ∩ Λ(−A) = ∅ they still solve the Lyapunov equations (12), see [40]. Hence, TLBT might be one possible approach to reduce unstable systems. Except for one short numerical experiment, we will not pursue this topic any further. In analogy to the infinite time horizon case, the eigenvalues of the product P T Q T are called time-limited Hankel singular values. A basic calculation shows that the time-limited Hankel singular 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 inhomogeneities, while the left hand sides are the same unchanged (adjoint) Lyapunov operators. This raises the question how much different the timelimited 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 factorizations Of course, before approximately solving the time-limited Lyapunov equations, the actions of the matrix exponentials e Ati to B (as well as e A T ti to C T ) have to be dealt with numerically. Numerical approaches for handling the matrix exponential and computing low-rank solution factors Z P T , Z Q 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 belowand using the low-rank solution factors Z P T , Z Q T in Algorithm 2: Required changes in Algorithm 1 for TLBT. (12).; the remaining steps. Depending on the time region of 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 2 -norm [40].

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 timelimited Gramians. Before algorithms for dealing with the matrix exponentials 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 be dealt with similarly by replacing A, B with A T , C T , respectively. Moreover, only the situation (8a) will be discussed, i.e., t s = 0 and 0 < t e < ∞.

The Difference of the Infinite and Time-limited Gramians
In order to approximate P T , Q T by a low-rank factorizations, it is desirable that their 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 λ 1 ≥ . . . ≥ λ n . The inhomogeneities of the time-limited Lyapunov equations (12) have 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 Due to the stability of A, the difference E(t e ) = P ∞ −P T will decay for increasing values of t e . Tight bounds for the eigenvalue behavior of P T and E(t e ) with respect to the parameter t e 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 e ) depends on t e . and is a Hermitian positive definite Cauchy matrix.

Proof. Apply the eigendecomposition of A and P
Consider the impulse response of (1), 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 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 difference between P T and P ∞ might be only observed when t e is chosen small enough in the sense that η(t e ) 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].

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,13,14,15,23,24,30,31,34,35,44]. Here, we are mostly interested in projection methods using (block) rational Krylov subspaces of the form where s k ∈ C + ∪ ıR ∪ {∞} are called shifts. They represent the poles of a rational approximation r k = ψ k−1 /φ k−1 of e z with polynomials ψ k−1 , φ k−1 of degree at most k − 1. Let Q k ∈ R n×km have orthonormal columns and range (Q k ) = RK k . Then a Galerkin approximation [44] of e Ate B takes the form Note that for m = 1, the rational function r k underlying the rational Krylov approximation Further information on the approximation properties can be found in, e.g., [6,10,14,15,16,30]. In the remainder we assume The orthonormal basis matrix Q k and the restriction H k can be efficiently computed by a (block) rational Arnoldi process [41]. Since H k is low-dimensional, e H k te can be computed by standard dense methods for the matrix exponential [34]. The choice of shifts s k (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 techniques [14,15,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 k+1 is selected via where z j are the eigenvalues of H k (Ritz values of A); s j , j = 1, . . . , k are the previously used shifts; and ∂S k marks a set of discrete points from the boundary of S k approximating Λ(A). We follow [16] and use S k as the convex hull of Λ(H 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 j in the denominator of r k 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).

Computing the low-rank Gramian factors
Using (14), (15), a Galerkin approximation of the time-limited Gramian is (cf. [43]). It holds B k =B 0,k = e H k ts B k for the special case t s = 0 considered here, indicating already how t s = 0 can be included as well. Hence, after B te,k of sufficient 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 T ,k leads to a sufficiently 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 k Z T k ≈ P T is illustrated in Algorithm 3. In the presence of a complex shift s k+1 , it is implicitly assumed that s k+1 is the subsequent shift. For this situation, the orthogonal expansion in Steps 5-9 of the already computed basis matrix Q k by real basis vectors goes back to [42]. The projected matrix H k 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 k Z T k and allow a rank truncation. Note that typically, once the approximation of B te is found, the generated subspace is already good enough to acquire a low-rank Gramian approximation without many additional iterations of the rational Arnoldi process. Often, the criteria in Steps 12,15 are satisfied in the same iteration step.

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 s , i.e., T = [t s , t e ], can also be done using Algorithm 3 with minor adjustments. Having computed a rational Krylov basis for approximating e Ate B for some time value t e , the same basis typically also provides good approximations for any other time values [31]. Consequently, Algorithm 3 has to be simply changed by addingB ts,k = e H k ts B k , B ts,k = Q kBts,k and appropriately adjusting the steps regarding the Gramian approximation. Also the methods relying on Taylor approximations [1] can be easily modified to handle multiple time values. However, even if a nonzero t s 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 s , t e ]. Already the original TLBT paper [25] states that TLBT in [t s , t e ] is expected to only give good approximations of the impulse response (u(t) = vδ(t), v ∈ R m ) y δ (t) = C e At 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 k ) = range (T ), e Ats B ≈ B k,ts ∈ range (Q k ), e Ate B ≈ B k,te ∈ range (Q k ), and likewise for C e Ats , C e Ats and range (S). Hence, the impulse response is accurately approximated at the relevant times.
For the response of an arbitrary input u(t), y u (t) = C t 0 e A(t−τ ) Bu(τ )dτ , such argumentation clearly does not automatically hold. The value of y u (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 u (t) will be poorly approximated for t ≤ t s 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 e − t s ]. However, this time translation will also introduce an inhomogeneous initial value x 0 , which is an additional difficulty for model order reduction. Some strategies to cope with nonzero initial values are given in [5,33] and we plan to investigate the incorporation to time-limited BT in the future.

Generalized state-space and index-one descriptor systems
Consider generalized state-space systems with M ∈ R n×n nonsingular. Simple manipulations reveal that, similar to unrestricted [7] and frequency-limited BT [8], the time-limited Gramians of (18) are P T , M T Q T M , where P T , Q T solve the time-limited generalized Lyapunov equations for t i = t s , t e . Note that B ti := e AM −1 ti B. Low-rank factors of P T , Q T can be computed by methods for generalized Lyapunov equations. In particular, the rational Krylov subspace methods which we employ for approximating B ti will implicitly work on has to be used in Step 2 of Algorithm 1. In some of the numerical experiments we will encounter the situation (18) becomes a semi-explicit index-one descriptor system. Eliminating the algebraic constraints leads to a general state-space system defined by M 1 ,Â : and an additional feed-through term −C 2 A −1 4 B 2 , see [22]. Time-limited and unrestricted BT can be applied right away to this system via (19)

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 situation, a stability preserving modification of TLBT is proposed in [29] using the Lyapunov equations where Q B ∈ R n×r B , Q C ∈ R n×r C contain the eigenvectors corresponding to the r B ≤ 2m, r C ≤ 2p nonzero eigenvalues λ B i , λ C i of the right hand sides of (12a),(19a) and, respectively, (12b),(19b). The rational Krylov subspace method in Algorithm 3 for M = I, t s = 0 can be easily modified for (20) by replacing Step 13 with the steps shown in Algorithm 4. Generalization for M = I and 0 < t s < t e are straightforward. Note that because 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 T , Q T . As observed 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.

Algorithm 4: Changes in Algorithm 3 for modified time-limited Gramians
13a Compute partial eigendecomposition

Numerical Experiments
Here, we illustrate numerically the results of Section 3 regarding the numerical rank of the frequencylimited Gramians as well as the numerical method for computing low-rank factors of P T , Q T . Afterwards, the quality of the approximations obtained 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 R 8.0.0.783 (R2012b) on a Intel R Xeon R 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 u(t) = δ(t)v, v ∈ R m , x 0 = 0 can be expressed as y δ (t) = C e At Bv, it is computed by the same integrator applied to the uncontrolled system (u(t) = 0) with initial condition x 0 = Bv. For the impulse, or step responses, the vector distributing the control to the columns of B is set to v = 1 m . 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].   Hankel singular values is visible. In the top right plot the point wise norm of the impulse response y δ (t) shows that after t e = 10, y δ (t) has already almost reached its stationary phase. This confirms the expectation that significant differences between infinite and time-limited Gramians occur only for times t e 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 e . The numerical ranks of P T , Q T clearly move towards the numerical ranks of P ∞ , Q ∞ as t e increases. In contrast, the numerical ranks of P mod T , Q mod T (Section 4.3) are always close to the ones of P ∞ , Q ∞ and appear to be largely unaffected by different values of t e . This is a very similar behavior as observed for the modified frequency-limited Gramians in [8].

Computing low-rank factors of the time-limited Gramians
We proceed by testing the computation of low-rank factors of the infinite and (modified) time-limited reachability Gramians by the rational Krylov subspace method in Algorithm 3. The stopping criteria for matrix function and Gramian approximations use the thresholds τ f = τ P = 10 −8 . To save some computational cost, the projected 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 e , the produced subspace dimensions d, the ranks r of the low-rank approximations, and the total computing times t rk for the Gramians P ∞ , P T , and P mod T . Apparently, larger rational Krylov subspaces and approximately twice as long computation Table 3: Model reduction results by BT, TLBT, and modified TLBT using different t e and u(t): reduced order r, largest relative error E T in [0, t e ], overall computation time t mor , and s ∈ {0, 1} indicates if reduced system is asymptotically stable or unstable. 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 T are in all cases smaller compared to P ∞ . The modified time-limited Gramians P mod T do not show this behavior as they have ranks similar to the infinite Gramians. For the rail example, increasing the time horizon t e does not change the required subspace dimension d for the approximation of the time-limited Gramian P T , but its rank r is clearly larger. In all cases, no additional steps of the rational Krylov method were necessary once the approximation of e Ate 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.

Model reduction results
Now we execute BT [36] as well as (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 constructed reduced order models is assessed via the point wise and maximal relative error norms of the output responses y(t),ỹ(t) of original and reduced order models. We consider the response y δ (t) for the impulse input u(t) = δ(t)1 m for all examples. Moreover, 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 m and u(t) = 501 m for the bips 3078 and, respectively, rail example, and u * := [5 · 10 4 · 0.198(sin(tπ/100) 2 ), 4, 2, 1, 3, 1] T for vertstand.
The results are given in Table 3 listing the largest relative error E T in [0, t e ] and the overall computation time t mor , i.e., the computation time for computing low-rank 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 e , the system responses and point wise relative errors E(t) are plotted against the time t in Figure 2 and Figure 3, respectively. Figure 4 shows the behavior of E T as the reduced order r increases.
Apparently, for the chosen orders r and in the time regions of interest, the largest relative errors E T of the reduced order models returned by TLBT are in most experiments more than one order of magnitude smaller compared to standard and modified time-limited BT. The plots in Figure 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 delivers 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 criterion (MAC) [20] (M AC(x, y) = |y T x| 2 /( x 2 y 2 )) indicated that the spaces spanned by the projection matrices T TLBT , S TLBT in TLBT and T BT , S BT in unrestricted BT, respectively, are different. For example, computing the MAC for the right projection matrices T TLBT , T BT ∈ R n f ×100 for the bips 3078 system showed that only a few of the columns of both matrices are well correlated to each other (i.e., M AC(T TLBT e i , T BT e j ) ≈ 1 for very few i, j ∈ {1, . . . , 100}).
Moreover, albeit the higher accuracy in [0, t e ], in some cases TLBT produces unstable reduced order models. This is especially visible in the upper left plot of Figure 2 showing the impulse responses of bips 3078. After leaving [0, t e ], the impulse response of reduced order model generated by TLBT exhibits an exponential growth and departs from the original response. As illustrated with the rail examples and proven in [25], using a higher end time t e can already cure this. Modified TLBT does not generate unstable reduced systems, but its approximation quality in [0, t e ] 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 e ], 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 s . 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 4 · 0.198(sin(tπ/100) 2 ), 4, 2, 1, 3, 1] T in the time window T = [t s , t e ] = [50, 100]. 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 Figure 5 illustrates the relative errors plots as well as    the largest relative error in T against the reduced order r. Compared to standard BT, TLBT delivers, as expected, more accurate results of the impulse response, but fails for u(t) = u * . Moreover, the bottom left plot Figure 5 shows that the decay of E T with respect to the impulse response is less monotonic as in the case t s = 0 s.t. 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 e ], the same cannot be said when time intervals [t s , t e ], t s > 0 are considered. Improving the performance of TLBT in this scenario is therefore subject of further research.
As 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 + 0.2M leading to max Re (λ) = 0.323, λ ∈ Λ(Â). The final time is set to t e = 5 and a reduced order model with r = 100 is generated. Due to the comparatively small system dimension (n f = 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 lowrank Gramians factors. The impulse response of exact and reduced system and the relative error are illustrated in Figure 6. Apparently, TLBT was able to reduce this mildly unstable system to a reduced order model with maximal relative error E T ≈ 1.2 · 10 −2 in [0, t e ]. Hence, provided t e is chosen in a reasonable way, e.g., with a sufficient distance to the exponential growth of y(t), TLBT appears to be a potential candidate from model order reduction of unstable systems. However, the applications to large-scale unstable systems is currently still difficult because algorithms for computing low-rank solutions of Lyapunov equations usually require that A is (anti)stable. Advances in this direction are, therefore, necessary to pursue this type of reduction further.

Conclusions
BT model order reduction for large-scale systems restricted to finite time intervals [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 increasing 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 exponentials 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 several orders of magnitude more accurate reduced order models in time intervals of the form [0, t e ] 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 s , t e ] with nonzero start times t s > 0 provided much worse results than with t s = 0, except for approximating the impulse response. In the form introduced in [25], TLBT appears to be incapable of producing good reduced order models with respect to [t s , t e ] and, thus, further research is necessary in this direction. The results in [5,33] might be one possible ingredient for this. A short experiment also indicated that TLBT can be employed to reduce unstable 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.