Frequency- and Time-Limited Balanced Truncation for Large-Scale Second-Order Systems

Considering the use of dynamical systems in practical applications, often only limited regions in the time or frequency domain are of interest. Therefor, it usually pays off to compute local approximations of the used dynamical systems in the frequency and time domain. In this paper, we consider a structure-preserving extension of the frequency- and time-limited balanced truncation methods to second-order dynamical systems. We give a full overview about the first-order limited balanced truncation methods and extend those methods to second-order systems by using the different second-order balanced truncation formulas from the literature. Also, we present numerical methods for solving the arising large-scale sparse matrix equations and give numerical modifications to deal with the problematic case of second-order systems. The results are then illustrated on three numerical examples.


Introduction
The modeling of, e.g., mechanical and electrical systems often leads to linear dynamical systems containing second-order time derivatives. In this paper, we consider linear second-order inputoutput systems of the form with M, E, K ∈ R n×n , B u ∈ R n×m and C p , C v ∈ R p×n , and u(t) ∈ R m , the inputs, x(t) ∈ R n , the states, and y(t) ∈ R p , the outputs of the system. In the frequency domain, the input-to-output relation is directly given as y(s) = H(s)u(s), whereby the so-called transfer function is given by with s ∈ C. In applications, the number of differential equations, n, describing the system, can become very large. This complicates using the model for simulations and controller design due with E, A ∈ R N ×N , B ∈ R N ×m , C ∈ R p×N , and the corresponding transfer function with s ∈ C. For simplicity, we are assuming that E is invertible and the system is c-stable, i.e., all eigenvalues of λE − A lie in the open left complex half-plane. The extension of the balanced truncation method to the descriptor system case (E non-invertible) can be found in [11,45]. The system Gramians of (3) are defined as with P ∞ the infinite controllability Gramian and E T Q ∞ E the infinite observability Gramian. Note that in the infinite case, the frequency and time representations of the Gramians are equal. It can be shown that those Gramians (5) are the unique, symmetric positive semi-definite solutions of the following Lyapunov equations The Hankel singular values are then defined as the positive square-roots of the eigenvalues of P ∞ E T Q ∞ E, which are a measure of how much influence the corresponding states have on the input-output behavior of the system. The main idea of balanced truncation is to balance the system such that with the Hankel singular values σ 1 ≥ σ 2 ≥ . . . ≥ σ N > 0 and then to truncate states corresponding to small Hankel singular values [34]. The complete balanced truncation square-root method is summarized in Algorithm 1. The balanced truncation method provides an a posteriori error bound in the H ∞ norm where H is the transfer function of the original model (4) and H the transfer function of the reduced-order model. The bound (7) depends only on the truncated Hankel singular values, which allows an adaptive choice of the reduction order. Also, this method preserves the stability of the original model, i.e., if H was a c-stable model then also H will be c-stable. The application of the balanced truncation method to large-scale sparse systems is possible by approximating the Cholesky factors of the Gramians via low-rank factors P ∞ ≈ Z R∞ Z T R∞ , Q ∞ ≈ Z L∞ Z T L∞ , with Z R∞ ∈ R N ×k R , Z L∞ ∈ R N ×k L and k R , k L N ; see, e.g., [9]. The approximation of the Gramians is reasonable due to a fast singular value decay arising by the lowrank right-hand sides [1]. For the computation of those factors, appropriate low-rank techniques are well developed [10].

Frequency-limited approach
A suitable method to localize the approximation behavior of the balanced truncation method in the frequency domain is the frequency-limited balanced truncation [20]. The idea is based on the

Algorithm 1: Balanced Truncation Square-Root Method
Input: System matrices A, B, C, E from (3). Output: Matrices of the reduced-order system A, B, C, E. 1 Compute Cholesky factorizations of the Gramians by solving the Lyapunov equations (6) such that with Σ 1 = diag(σ 1 , . . . , σ r ) containing the r largest Hankel singular values. 3 Construct the projection matrices 1 . 4 Compute the reduced-order model by frequency representation of the system Gramians (5), such that the frequency-limited Gramians of (3) are given by where Ω = [−ω 2 , −ω 1 ] ∪ [ω 1 , ω 2 ] ⊂ R is the frequency range of interest. It can be shown that the left-hand sides in (8) are also given as the unique, symmetric positive semi-definite solutions of the two Lyapunov equations AP Ω E T + EP Ω A T + B Ω B T + BB T Ω = 0, with new right hand-side matrices B Ω = EF Ω B, C Ω = CF Ω E containing the matrix functions with ln(.) the principle branch of the matrix logarithm. Note that in case of Ω = [−ω, ω], the function evaluation (10) simplifies to Also, the frequency-limited Gramians can be extended to an arbitrary number of frequency intervals, i.e., for  (9) such that P Ω = R Ω R T Ω , Q Ω = L Ω L T Ω . 2 Follow the steps 2-4 in Algorithm 1.
with 0 < ω 1 < . . . < ω , leads to the following modification of (10) See [6] for a more detailed discussion of the theory addressed above. The extension of this method to the large-scale system case can also be found in [6] and an extension to descriptor systems in [26]. The resulting frequency-limited balanced truncation method is summarized in Algorithm 2.

Time-limited approach
The counterpart of the frequency-limited balanced truncation from the previous section in the time domain is the time-limited balanced truncation [20]. This method aims for the approximation of the system on a time interval T = [t 0 , t f ], where 0 ≤ t 0 < t f , based on the limitation of the time domain representation of the Gramians (5). The time-limited Gramians of (3) are then given by and it can be shown, that the left-hand sides in (11) are the unique, positive semi-definite solutions of the two following Lyapunov equations where the new right hand-side matrices At 0/f contain the matrix exponential. The right hand-sides of (12) simplify in case of t 0 = 0 since B 0 = B and C 0 = C. A more detailed discussion of the time-limited theory, especially for the large-scale system case, can be found in [28]. Also, the extension of the theory to the case of descriptor systems is given in [22]. It can be noted that considering more than one time interval at once [t 0,1 , t f,1 ] ∪ · · · ∪ [t 0, , t f, ] is not practical and usually one cannot guarantee a good approximation behavior in the single intervals. Instead it is common to take the smallest and largest time points in the intervals to construct a new overarching time interval [t 0,min , t f,max ], where t 0,min = min{t 0,1 , . . . , t 0, } and t 0,max = max{t f,1 , . . . , t f, } such that The resulting time-limited balanced truncation method is summarized in Algorithm 3.

Second-order case
After recapitulating the basic ideas of the classical as well as the frequency-and time-limited balanced truncation methods for first-order systems, in this section we will extend those methods to second-order systems (1).

Second-order balanced truncation methods
Over time, there have been many attempts for the generalization of the classical balanced truncation method to the second-order system case [16,33,37]. All of them have in common the idea of linearization, i.e., the second-order system (1) is rewritten as a first-order system. The usual linearization of choice for (1) is its so-called first companion form where T is the new combined state vector. The matrix J ∈ R n×n is an arbitrary invertible matrix but usually chosen as J = I n or J = −K, which can lead to symmetric A and E matrices in case of mechanical systems. For system (13), the first-order Gramians are used, as given by (5) or (6), and then partitioned according to the block structure in (13) such that where P p , Q p are the the position Gramians of (1) and P v , Q v the velocity Gramians. Due to P ∞ = P T ∞ ≥ 0 and Q ∞ = Q T ∞ ≥ 0, also the position and velocity Gramians are symmetric positive semi-definite and can be written in terms of their Cholesky factorizations Based on those, the different second-order balanced truncation methods are defined by balancing certain combinations of the four position and velocity Gramians. For most of the methods, the resulting balanced truncation is computed as second-order projection method where the different choices for W and T can be found in Table 1. There, the different transformation formulas are summarized and denoted by the type as used in the corresponding references. The subscript 1 matrices denote the part of the singular value decompositions corresponding to the largest characteristic singular values. In contrast to the balancing methods that describe the reduced-order model by (15), the secondorder balanced truncation (so) from [16] computes the reduced-order model by

SVD(s) Transformation Reference
3 Compute the singular value decompositions and transformation matrices as in Table 1. 4 Compute the reduced-order model by either (15) for the methods p, pm, pv, vp, vpm, v and fv or by (16) for so.
where S = W p JT v and the transformation matrices W p , W v , T p , T v are given in the last line of Table 1. This type of balancing can be seen as a projection method for the first-order realization (13) with a recovering of the second-order structure. The general second-order balanced truncation square-root method is summarized in Algorithm 4.

Remark 1.
In contrast to the first-order balanced truncation described in Section 2.1.1, none of the second-order balanced truncation methods provides an error bound in the H ∞ norm or can preserve the stability of the original system in the general case. A collection of examples for the stability issue is given in [37]. In case of symmetric second-order systems, i.e., M = M T , E = E T , K = K T , C p = B T u and C v = 0, it can be shown that the position-velocity balancing (pv) as well as the free-velocity balancing (fv) are stability preserving. Note that the position-velocity balancing also belongs to the class of balanced truncation approaches, which define system Gramians by using the underlying transfer function structure (2). Those balancing approaches have been generalized in [14] for systems with integro-differential equations.

Second-order frequency-limited approach
The generalization of the frequency-limited balanced truncation method for second-order systems has been discussed in [23] for the position (p) and position-velocity (pv) balancing from [37]. Here we will summarize their results and give a more general extension for the frequency-limited secondorder balanced truncation method. The basic idea for the approach comes from the observation that the block partitioning of the Gramians (14) can be written as Therefor, the extension of the existing second-order balanced truncation methods to the frequencylimited approach can be done by replacing the infinite first-order Gramians P ∞ and Q ∞ in (17) by the first-order frequency-limited Gramians P Ω and Q Ω from (8) corresponding to the first-order realization (13). The frequency-limited second-order Gramians are then given by where P Ω,p and P Ω,v are the frequency-limited position and velocity controllability Gramians, and J T Q Ω,p J and M T Q Ω,v M are the frequency-limited position and velocity observability Gramians. Note that P Ω and Q Ω are given by (9) using the first-order realization (13). As for the infinite Gramians, one observes that the frequency-limited position and velocity Gramians are symmetric positive semi-definite. According to [20,23,37], we can now define the corresponding frequency-limited characteristic values as follows. Following the observations in the first-order frequency-limited case as well as the second-order balanced truncation method, those characteristic singular can be interpreted as a measure for the influence of the corresponding states to the input-output behavior of the system in the frequency range of interest. Anyway, there is no energy interpretation as for the first-order balanced truncation method.
With (18) and the Definition 1, the resulting second-order frequency-limited balanced truncation square-root method is written in Algorithm 5.

Remark 2.
The second-order frequency-limited balanced truncation method is in general not stability preserving. Also, the approach from [23] does not necessarily lead to a one-sided projection as suggested by the authors and also might not produce a stable second-order system in the end. Even so, we will discuss approaches that can have stability-preserving properties in Section 3.4.

Algorithm 5: Second-Order Frequency-Limited Balanced Truncation Square-Root Method
Input: System matrices M , E, K, B u , C p , C v from (1), frequency range of interest Ω. Output: Matrices of the reduced-order system M , E, K, B u C p , C v . 1 Compute Cholesky factorizations of the first-order frequency-limited Gramians by solving (9), where the linearization (13) is used, such that P Ω = R Ω R T Ω , Q Ω = L Ω L T Ω . 2 Follow the steps 2-4 in Algorithm 4.
Compute Cholesky factorizations of the first-order time-limited Gramians by solving (12), where the linearization (13) is used, such that Follow the steps 2-4 in Algorithm 4.

Second-order time-limited approach
The extension of the time-limited balanced truncation to the second-order system case was first discussed in [24]. As in the previous section, we are generalizing the ideas from [24] to all secondorder balanced truncation methods. In any case, the same idea as for the frequency-limited case is applied here. That means, we replace the infinite first-order Gramians in (17) by the first-order time-limited Gramians from (11) to get where again the first-order realization (13) was used. Following the naming scheme of [37], P T,p and P T,v are the time-limited position and velocity controllability Gramians, and J T Q T,p J and M T Q T,v M the time-limited position and velocity observability Gramians. Note that P T and Q T are given by (12) with the first-order realization (13). As for the infinite Gramians, one observes that the time-limited position and velocity Gramians are symmetric positive semi-definite. According to the frequency-limited characteristic singular values, we are giving the following definition for the time-limited version.
on the first-order case and does not guarantee the preservation of stability for second-order systems in general. Approaches that can be more beneficial in terms of preserving stability are discussed in Section 3.4.

Numerical methods
In this section, we will discuss points concerning the numerical implementation of the proposed second-order frequency-and time-limited balanced truncation methods.

Matrix equation solvers for large-scale systems
A substantial part of the numerical effort in the computations of the second-order frequencyand time-limited balanced truncations goes into the solution of the arising matrix equations (9) and (12). In general it has been shown for the first-order case, that the singular values of the frequency-and time-limited Gramians are decaying possibly faster than of the infinite Gramians; see, e.g., [6] for the frequency-limited case. That leads to the natural approximation of the Gramians by low-rank factors, e.g., N . Those low-rank factors then replace the Cholesky factors in the balanced truncation algorithms 1-6.
In the following three sections, we will shortly review existing approaches for these problems and give comments on existing implementations.

Quadrature-based approaches
A natural approach based on the frequency and time domain integral representations of the limited Gramians (8) and (11) is the use of numerical integration formulas. As used for example in [23,26], the low-rank factors of the Gramians can be computed by rewriting the full Gramians by quadrature formulas, e.g., where γ k are the weights and ω k the evaluation points of a fitting quadrature rule, which can be again rewritten for the low-rank factors by Note that this approach becomes unhandy considering the timelimited case, since there, for each step of the quadrature rule, an approximation of the matrix exponential has to be computed. A different approach was suggested in [6], which writes the right-hand side of the frequencylimited Lyapunov equations (9) as integral expressions, such that the right-hand side is first approximated and afterwards the large-scale matrix equation is solved, using one of the approaches in Section 3.1.2 or 3.1.3. In general it is possible to approximate the right-hand sides of (9) and (12) with matrix functions by using the general quadrature approach from [25]. We are not aware of a stable, available implementation of quadrature-based matrix equation solvers for the frequencyand time-limited Lyapunov equations and, therefor, use the following approaches rather than the quadrature-based methods.

Low-rank ADI method
The low-rank alternating direction implicit (LR-ADI) [8,31] method is a well established procedure for the solution of large-scale sparse Lyapunov equations. Originally developed for the Lyapunov equations corresponding to the infinite Gramians (6), the LR-ADI produces low-rank approximations of the form [4][5][6] for more details on this method.
The right-hand sides of the limited Lyapunov equations (9), (12) can be rewritten as , which shows that the righthand side matrices are indefinite. The LR-ADI method can be extended to this case by using an LDL T -factorization for the right-hand side as well as for the solution [29]. Note that for applying this method for the solution of the large-scale matrix equations, an approximation of the matrix functions in the right-hand sides is needed beforehand. It was noted in [6], that the information used for the approximation of the matrix functions cannot be used in the LR-ADI method. A stable version of the LR-ADI method in the low-rank and LDL T formats is implemented in [39]. We will use this implementation in case the methods, described in the following section, are failing to converge for the solution of the matrix equation but give approximations to the function right hand-sides.

Projection methods
An approach that can be used to approximate the matrix functions in the right-hand sides of the limited Lyapunov equations, as well as to solve the large-scale matrix equations at the same time, is given by projection-based methods. Here, low-dimensional subspaces V k = range(V k ) are used to obtain the low-rank solutions as, e.g., whereP Ω is the solution of the projected Lyapunov equation are the projected matrices of the frequencylimited controllability Lyapunov equation (9). The equation (20) is now small and dense and can be solved using established dense solvers. As one can observe, this method gives also the opportunity to approximate the matrix function right-hand side by the low-dimensional subspace V k , for which one can also use dense computation methods [25].
Usually, the low-dimensional subspace V k is constructed as standard [27], extended [44] or rational Krylov subspace [17], all of which can be easily computed for large-scale sparse systems. The implementation of the limited balanced truncation methods for second-order systems [13], we provide, is also based on rational Krylov subspaces. We refer the reader to [6, Algorithm 4.1] for the underlying idea of the implementation.
A drawback of the projection-based approach, especially for second-order systems, is that the projected system matrices T k are not necessarily c-stable, even if the original first-order realization of the second-order system was. In fact, the quality and performance of the projection-based solvers strongly depend on the chosen first-order realization. Therefor, we are going to use the so-called strictly dissipative realization of second-order systems [36] in our computations. Assuming M, E, K to be symmetric positive definite, the second-order system (1) can be described by a first-order realization using the following matrices with the parameter 0 < γ < λ min (E(M + 1 4 EK −1 E) −1 ). The advantage of this realization is that E is symmetric positive definite and A + A T symmetric negative definite. Following that, projection methods can preserve the stability in the projected matrices T k if the computations are made on the corresponding standard state-space realization, obtained by a symmetric statespace transformation using the Cholesky factors E = LL T , i.e., the algorithms work implicitly on a realization of the formq Remark 4. Note that the realization (21) is computationally more involved than the classical first companion form (13) or its second companion form, since it is not possible to make use of occurring zeros in the block structure.
Also, by changing the first-order realization to (21), the computed Gramians change compared to the definition of the second-order balancing methods. Therefor, let P and Q be Gramians computed for the strictly dissipative first-order realization (21) and P and Q be the Gramians from the first companion form realization (13). Then it holds P = P and Q = T T QT, with the transformation matrix That means we can use the strictly dissipative realization (21) for the solution of the matrix equations and for the balancing procedure just perform the easy back transformation of the observability factor.

Stabilization and acceleration by α-shifts
So far, it was always assumed that the second-order system (1) is c-stable. But in practice, the eigenvalues of λ 2 M + λE + K can be very close to the imaginary axis or even on the axis, e.g., in the case of marginal stability. This makes the usage of the model reduction methods and matrix equation solvers very difficult. A strategy to overcome those problems has been proposed in, e.g., [19]. There, a shift in the frequency domain was used to move the spectrum of the pencil λE − A, which had eigenvalues at zero, away from the imaginary axis to compute the system Gramians. This approach cannot be used the same way for the first-order realizations (13) or (21) of second-order systems since it destroys the block structure one can exploit in the numerical implementations of the solvers or rather the block structure that is used for the second-order balancing approaches. Therefore, we will transfer the concept of α-shifts to the case of secondorder systems. Let α ∈ R >0 be a real, strictly positive shift and consider the second-order differential equations in the frequency-domain with E = E + 2αM and K = K + αE + α 2 M . Also, the second equation (22b) can be rewritten as where C p = C p + αC v . Now, the new system described by (M, E, K, B u , C p , C v ) is used for the computation of the reduced-order projection matrices W, T ∈ R n×r . Then, the projected system ( M , E, K, B u , C p , C v ) yields the following relations where E = W T ET , K = W T KT and C p = C p T are the transformed non-shifted matrices. Now, we consider the transformed system again in the frequency domain with the Laplace variable ρ and using the back-substitution ρ = s − α, such that The back-substitution gives the resulting reduced-order model ( M , E, K, B u , C p , C v ). The α-shift strategy can be interpreted as a structured perturbation in the frequency domain during the computations. Experiments have shown that such an approach works fine for α small enough. It has to be noted that there are no theoretical results on the influence of the chosen α concerning the quality of the reduced-order model or properties like stability preservation and error bounds.

Remark 5.
The α-shift approach can also be used either to improve the conditioning of the used matrix equation solvers by improving the condition number of the shifted linear systems solving with (σ 2 M + σ E + K), or to improve the convergence of those solvers by pushing the eigenvalues of λ 2 M + λ E + K further away from the imaginary axis.

Two-step hybrid methods
The idea of two-step (or hybrid) model reduction methods has been used for quite some time in different applications [18,30,46]. In general, two-step methods are based on the division of the model reduction process into two phases. First a pre-reduction, which can be easily computed and gives a very accurate approximation for the system's behavior. The model resulting from the pre-reduction is usually of medium-scale dimensions, on which the second reduction step by a more sophisticated model reduction method is applied. This procedure has the advantage that there is no necessity of applying difficult approximation methods for the large-scale matrix equations arising in the balancing related approaches. Instead, the exact methods can be used on the, usually dense, pre-reduced system. In order to have a structure-preserving pre-reduction method, we suggest the use of interpolation by rational Krylov subspaces [2,41,42]. This has been shown to be equivalent to the use of shiftbased approximation methods for the large-scale matrix equations in Section 3.1; see [46]. The second-order rational Krylov subspaces are generated as with s k ∈ C, k = 1, . . . , , chosen interpolation points. Let V and U be Hermitian bases of the same size such that V ⊂ range(V ) and U ⊂ range(U ), respectively, the pre-reduced model is then generated by For preservation of stability and the realness of the system matrices, we choose the interpolation points to appear in complex conjugate pairs s k and s k , and replace one of the projection matrices by U = V . The choice of points s k is crucial for the quality of the pre-reduced model. While there are strategies for an adaptive or optimal choice of s k , we suggest a simple oversampling on the imaginary axis, which is usually enough as a global pre-reduced model.

Remark 6.
For the frequency-limited case, a natural choice for the interpolation points would be to take jΩ instead of aiming for a global approximation. In this case, the resulting frequencylimited balanced truncation will very likely not give the same results as the large-scale approach. This observation comes from the fact, that the frequency-limited balanced truncation still takes information about the complete system structure into account and the pre-reduced system can be completely different from the original one, if only a local pre-reduction is performed.
Due to the required accuracy of the pre-reduced model, the dimension of it can be still very large. Therefore, we suggest an efficient iterative solver for the Lyapunov equations appearing in the second reduction step. In general, we consider the following stable Lyapunov equations where Q ∈ R m×m and R ∈ R p×p are symmetric and possibly indefinite. The solution of (23) can then be factored in the same way as the right-hand sides, i.e., where Y 1 and Y 2 are also symmetric matrices. For efficiently computing the solutions of (23), we extend the dual sign function iteration method from [3] for the LDL T -factorization of the solutions. As a result, we get a sign function iteration, that solves both Lyapunov equations with symmetric indefinite right hand-sides (23) at the same time; see Algorithm 7.
The implementation of Algorithm 7 as well as dense versions of the second-order frequency-and time-limited balanced truncation methods can be found in [12].

Remark 7. In
Step 4 of Algorithm 7, the memory requirements and operations are doubling in every step due to the extension of the solution factors. It is suggested to do LDL T column and row compressions at that point to keep the size of the factors small.

Modified Gramian approach
A drawback of the frequency-and time-limited balanced truncation methods is the loss of stability preservation. For the first-order system case, there are different modifications of the methods to regain the preservation of stability, e.g., the replacement of one of the limited Gramians by the infinite Gramian [22,26].
A different technique, proposed in [21], is the modified Gramian approach. Therefor, the indefinite right-hand sides (19) are replaced by definite ones. Using eigenvalue decompositions, the right-hand sides can be rewritten as
Compute the scaling factor for convergence acceleration Compute the next iterates of the solution factors Compute the next iteration matrix Set k = k + 1.

end 8 Construct the solution factors
as the solutions of the following Lyapunov equations Using those modified Gramians for the limited balanced truncation methods also preserves the stability in the reduced-order models in the first-order case. There also exists an H inf error bound for the modified frequency-limited balanced truncation for first-order systems [6]. Note that the limited Gramians can also be easily computed using the projection-based matrix equation solvers with only minor changes in the algorithm [6,28].

Remark 8.
Neither the replacement of limited Gramians by the infinite ones nor the modified Gramian approaches are guaranteed to preserve the stability in the reduced-order model when it comes to the second-order case. The stability preserving methods in [23,24] are just based on the assumption, that the same procedure as in the first-order case also works for second-order systems. This is not the case, since already the classical second-order balanced truncation methods are in general not stability preserving [37]. Figure 1: Setup of the single chain oscillator.

Remark 9.
Also, it has been mentioned and shown by numerical examples in [6,28] that the modified Gramian approach usually does not pay off since the quality of the reduced-order models is often the same as for the global approaches, i.e., the local approximation property of the limited balanced truncation methods gets lost.

Numerical examples
In the following, some mechanical systems of second-order form from the literature have been chosen as benchmark examples. The experiments reported here have been executed on machines with 2 Intel(R) Xeon(R) Silver 4110 CPU processors running at 2.10GHz and equipped with either 192 GB or 384 GB total main memory. The computers are running on CentOS Linux release 7.5.1804 (Core) and using MATLAB 9.4.0.813654 (R2018a). For the computations, the following software has been used: • MORLAB version 5.0 [12], for all evaluations in the frequency and time domain, the generation of the pictures and the dense implementations of the limited model reduction methods used in the two-step approach, • the limited balanced truncation for large-scale sparse second-order systems code package [13], for the computations of the full-order limited Gramians and the implementation of the balancing formulas from Table 1, • the M-M.E.S.S. library version 2.0 [39], for computing the full Gramians with already approximated right hand-sides.
In general, we used the projection-based methods from [13] to approximate the right hand-sides and the Gramians. But in case that the Gramians did not converge, we used the computed approximation of the right hand-sides from the projection methods in the ADI method from [39] to compute a solution to the matrix equation. where [t min , t max ] is the time frame as shown in the plots or rather the local time range [t 0 , t f ] chosen for the time-limited methods.
As criterion for the computed approximation order, the characteristic values from Definition 1 and 2 have been used. Therefore, we truncated all states corresponding to the singular values that in sum were smaller than the largest singular values multiplied with the tolerance 10 −4 , i.e.,

Single chain oscillator
As first example, we consider the single chain oscillator benchmark from [32], where we removed the holonomic constraint to get a mechanical system without algebraic parts. Figure 1 shows the basic setup of the system, where the parameters are chosen as in [32], i.e. in our experiments we have

Frequency domain
The frequency range of interest in this example is chosen, just for demonstration reasons, to be between 1 and 100 Hz. The computations have been done with no α-shift (α = 0). In Figure 2, the resulting reduced-order models (ROMs) can be seen in terms of their transfer functions (a), the point-wise absolute error (b) and point-wise relative error (c). The frequency range of interest is marked as the area between the dashed vertical lines. Table 2 gives an overview for all applied second-order frequency-limited. It can be noted that all computed ROMs are of order 2, stable and have absolute and relative errors in the same order of magnitude. Also we note that as wanted, the errors in the frequency range of interest are significantly smaller than in the overall considered frequency region. For the two-step approach, we used, on the one hand, a logarithmically  Since no significant differences between the full-order Gramian and two-step approaches could be seen, we refer the reader also to Figure 2 and Table 2 for the results.

Time domain
In the time domain, we apply two different input signals to test our ROMs u step (t) = δ(t − 5) and u sin (t) = sin(t)δ(t − 5), for t ∈ [0, 100] and δ(t) the Heaviside function. As time range of interest, [0, 20] has been chosen. While Figure 3 shows the results for the time-limited balanced truncation methods in terms of absolute and relative errors for the two applied input signals (24), in Table 3, the ROM sizes, absolute and relative errors are given. One can observe that all ROMs are of order 4, stable and have locally significantly smaller errors than globally.
Again, the result of the two-step approaches are only marginal distinguishable from the results of the full-order Gramians, where we used the global sampling between 10 −4 and 10 4 Hz to preapproximate the system's behavior. Therefore, those results are also not shown here.

Crankshaft
The crankshaft is a model from the University Stuttgart, describing the crankshaft of a fourcylinder engine [35], which is shown in Figure 4. After discretization by the finite element method, the constraint model is of dimension n = 42 126 with m = p = 35 inputs and outputs. Due to the rigid elements, coupling the interface nodes, the system has several eigenvalues at zero. Therefore, we apply the shift α = 0.01, as suggested in Section 3.2, to make the system asymptotically stable during the computations of the matrix equations and low-rank projection matrices.

Frequency domain
In the frequency domain, we are interested in the actual working range of the crankshaft between 4 and 20 kHz. Figure 5 shows the results for using the full-order frequency-limited Gramians. The frequency range of interest lies again between the two vertical dashed lines. We can see that all ROMs approximate the frequency region of interest better than the global region. Also Table 4 shows the desired approximation behavior in terms of the errors. In this example, some of the computed ROMs are unstable as denoted by x-marks in Table 4. It should be noted that even for the same order some methods might produce unstable models while others do not.
In this example, we also applied the two-step approach with 200 frequency sample points in the region of interest to generate the intermediate model of order 447. Those results can be seen in Figure 6. Table 5 shows that the ROMs produced by the two-step approach are slightly larger in dimension and also partially in errors, while the same methods (pm, vp, vpm, so) as for the full-order Gramian approach produce unstable models.

Time domain
In the time domain, we consider just the first 0.01 s of using the crankshaft, while the full simulation runs over a time range of [0, 0.05] s. As test input signals, we apply where 1 35 denotes the ones vector of length 35. The results for the time-limited balanced truncation with the full-order Gramians can be seen in Figure 7 and Table 6. Only one unstable model (vpm) was computed, which still gives suitable approximation results, and all ROMs have small enough errors in the time domain. Even so, we recognize that the local approximation error is only in some cases a bit smaller than the global one. For the two-step approach, we computed 200 logarithmically equidistant distributed samples in the frequency domain between 10 −2 and 10 6 Hz. The intermediate model had the order 876. Since the resulting ROMs are of the same order as the ones computed via the full-order Gramians, featuring the same stability properties, and are only slightly worse in terms of the time domain errors than in Table 6, we skip the additional presentation of those results here.

Artificial fishtail
The artificial fishtail is a mechanical system, describing the movement of a fishtail-shaped structure by using the fluid elastomer actuation principle. Figure 8 shows a transparent sketch of the fishtail model consisting of a carbon beam in the center and a silicon hull around. A more detailed description of the model as well as a comparison of structure-preserving second-order model reduction techniques for this example can be found in [40]. After spatial discretization by the finite element method, the resulting second-order system has n = 779 232 states describing the model. By the actuation principle, we have m = 1 input and a sensor is measuring the displacement of the fishtail's tip in all spatial dimensions, i.e., we have p = 3 position outputs and no velocity outputs. The discretized data is available as open benchmark at [43]. The computations were done without an α-shift (α = 0).

Frequency domain
In the frequency domain, the range of interest for the fishtail model lies between 0 and 20 Hz, since higher frequencies are physically not realizable. Figure 9 shows the results for the frequency-limited balanced truncation methods, based on the full-order Gramians. Except for the fv balancing there is no visible difference between the ROMs and the full-order model. The error plots show that the approximation reached a sufficiently small error in the region of interest. Table 7 shows the corresponding maximum absolute and relative error in the local and global frequency regions. It is remarkable that the methods were able to approximate the original model, having around 780 000 states, by stable order 1 systems in the region of interest. While the absolute errors are comparable between local and global region, the relative errors show again the strength of the frequency-limited method. Time t (sec)

Time domain
In the time domain, the fishtail is simulated from 0 to 2 s. For our time-limited methods we consider the time range up to 0.5 s and as inputs, the following two signals are considered u step (t) = 5000δ(t − 0.1) and u sin (t) = 2500(sin(10π(t − 1.35)) + 1)δ(t − 0.1). Figure 10 and Table 8 show the results. Except for the models generated by pm, vpm and fv, the computed ROMs have acceptable small errors in the time domain. Also, only the vpm ROM is unstable. The errors in the local region are sometimes a bit smaller than the global one as we were aiming for by the method. The two-step approach here used 200 logarithmically equidistant sample points in the frequency range from 10 −4 to 10 4 Hz, which gave an intermediate model of order 100. The results of the ROMs computed by the two-step approach differ a bit from the ones generated by the full-order Gramians. Those results can be seen in Table 9. There, shown errors are partially smaller or larger than in Table 8 and also we note that for the two-step approach, the vpm model is also unstable but still gives usable results for both applied input signals.