On shift selection for Krylov subspace based model order reduction

Mechanical systems are often modeled with the multibody system method or the finite element method and numerically described with systems of differential equations. Increasing demands on detail and the resulting high complexity of these systems make the use of model order reduction inevitable. Frequently, moment matching based on Krylov subspaces is used for the reduction. There, the transfer functions of the full system and of the reduced system are matched at distinct frequency shifts. The selection of these shifts, however, is not trivial. In this contribution we suggest an algorithm that evaluates an increasing number of shifts iteratively until a reduced model that approximates the full model in a subspace with very low approximation error is found. Thereafter, the projection matrix that spans this subspace is decomposed with singular value decomposition and only most important directions are retained. In this way, small reduced models with good approximation properties that do not exceed a predefined error bound can be found or low-error models for a given reduced order can be generated. The evaluation of more shifts than necessary and further reduction by means of singular value decomposition is the novelty of this contribution. In this paper, this novel approach is extensively studied and, furthermore, applied to the numerical example of an industrial helicopter model.


Introduction
In engineering today, virtual prototyping and digital twins are frequently used to reduce production and development costs or to enable model-based control algorithms. Therefore, the P. Eberhard peter.eberhard@itm.uni-stuttgart.de L. Frie lennart.frie@itm.uni-stuttgart.de 1 Institute of Engineering and Computational Mechanics, University of Stuttgart, Pfaffenwaldring 9, 70569 Stuttgart, Germany dynamics of the underlying mechanical systems are often modeled with Elastic Multibody Systems (EMBS) or the Finite Element Method (FEM). This results in a set of ordinary differential equations that describe the system behavior, see [22]. The desire for increasingly accurate predictions on the one hand and the increasing complexity of the systems on the other hand leads to very high orders of the differential equations. Despite the simultaneous increase in computing power, it is usually necessary to reduce the number of equations to enable efficient system evaluations. This is the basis for applications in the area of, e.g., parameter studies, optimization, uncertainty analyses, or control, see [12,13,20]. The approximation of high order of these sets of differential equations with low-rank systems is called Model Order Reduction (MOR) and has become an active field of research over the past years.
One method that is frequently used for the reduction of mechanical systems is moment matching with Krylov subspaces. In [12] it is used to reduce a helicopter model similar to the one used in this article and in [5,19] for other systems. In moment matching, the transfer function of the full order model is matched at chosen shifts in the frequency domain. The reduced system depends strongly on the chosen shifts and, yet, there is no general automatic procedure to choose them. In the basic form of Krylov reduction, the number of shifts is proportional to the rank of the reduced system and, therefore, fewer shifts are desirable to obtain a small reduced system. On the other hand, few shifts often lead to large approximation errors in areas not covered by the shifts. This creates the issue of choosing just as many shifts as necessary and placing them in the right places. Frequently, the shifts are linearly spaced in the frequency domain or set from experts based on experience or trial and error. Apparently, this choice is often far from optimal. An approach to overcome this issue was first presented in [14] with an algorithm that selects shifts iteratively in a way that the H 2 -norm of the approximation error is minimized. This algorithm is called Iterative Rational Krylov Algorithm (IRKA). It is extended to multiple-input-multiple-output (MIMO) systems in [7] and to second order systems in [21]. However, the IRKA does not ensure convergence, and on convergence it is only locally optimal, see [10]. Another drawback of this approach is that it is not possible to weight frequency intervals in terms of their importance, which is often useful for the approximation of mechanical systems, see, e.g., [15].
In this paper, we suggest a novel, practical, and automated approach to iteratively select the shifts for moment matching with Krylov subspaces. For this purpose, many shifts are first evaluated to find a possibly large subspace in which the original model is represented with a very small error. The dimension of this subspace is thereafter reduced by decomposing the projection matrix that spans this subspace with Singular Value Decomposition (SVD), choosing only most important spatial directions to be kept. This combination of Krylov subspace based moment matching and SVD is the novelty of this contribution. The functionality and merits of the presented approach are illustrated by academic examples throughout this paper and compared to state-of-the-art methods as well. Furthermore, the applicability to a numerical model of an industrial helicopter airframe is shown.
This article is structured as follows. Sect. 2 reviews the theoretical foundations of projection based model order reduction with Krylov subspaces and of SVD. Sect. 3 shows the dependency of reduced models on the selection of shifts by using state-of-the-art methods. In Sect. 4 a novel algorithm is presented, where as small as possible low-rank model that does not exceed a predefined error bound is generated. In Sect. 5 this approach is inverted by determining a small-error model for a given reduced order. The presented method is then applied to a helicopter model in Sect. 6, which is followed by the conclusions.

Theoretical background
The following section contains some theoretical foundations on the mathematical description of dynamical systems, projection based MOR for these systems, moment matching with Krylov subspaces, and matrix approximation with SVD.

Projection based model order reduction with Krylov subspaces
Mechanical systems are often spatially discretized and then described with the FEM, see [22]. Also elastic bodies in an EMBS are described individually by a finite element (FE) model and then connected to form a multibody system. This leads to linear second order differential equations Here, M, D, K ∈ R N×N are the mass, damping, and stiffness matrices, respectively. The state vector q(t) ∈ R N describes the nodal displacements. The notationq(t) := dq(t) /dt and q(t):= d 2 q(t) /dt 2 is used for the time derivatives in this paper. The right-hand term f (t) contains all external forces that act on the system. The time dependency is omitted in the following for better readability.
By writing the force vector as product f = Bu and extracting system outputs with y = Cq, Equation (1) can be rewritten as the second order input-output system The matrices B ∈ R N×b and C ∈ R c×N are the input and output matrices of the system and the vectors u ∈ R b and y ∈ R c are the inputs and outputs. In this work we only consider systems where the inputs and outputs coincide, i.e., b = c and, more importantly, B = C.
The accurate modelling of complex mechanical systems usually requests a fine discretization and, thus, results in very high orders N of the system of differential equations. In industrial applications they not uncommonly exceed millions. MOR is used to develop a reduced dynamical system with n N differential equations that appropriately describes the system behavior, i.e., it approximates the full model well in a defined space of interest. One possibility to build such models is linear projection based MOR. Its basic idea is to approximate the state vector within a subspace V = span{V } ⊂ R N with

Plugging Equation (3) in Equation (2) leads to
where r is a residuum that arises from the approximation in Equation (3). The matrix V is called projection matrix and is a core element of the projection based MOR. Leftmultiplying by a matrix W with makes Equation (4) specific. In this work, only orthogonal projections with W = V are considered because they have some advantages as, e.g., preservation of symmetry and stability properties. This choice also satisfies Equation (5) since r is orthogonal to the subspace V , see [11]. By multiplying V , one obtains the reduced model by the so-called Galerkin projection withMq +Dq +Kq =Bu, The tilde indicates quantities of the reduced model and its quite small system matrices result fromM The question remains how to choose V . A commonly used method that often works very well for approximating large systems is moment matching. It aims at approximating the transfer function H (s) = C s 2 M + sD + K −1 B of the system, where s is the Laplace variable. Therefore, the transfer function is written as a power series around expansion points s k , also called shifts, with as explained in [9,11]. The first J k moments are then matched at different shifts s k . This is implicitly done by the use of Krylov subspaces.

Definition 1 A second order Krylov subspace is defined as
The matrices A 1 , A 2 ∈ R N×N and G ∈ R N×r are called basic blocks and starting matrix.
For the reduction of Equation (2), the matrices A 1 , A 2 , and G are selected as Choosing the projection matrix V so that ensures the transfer function of the reduced system to match the transfer function of the full system at the shifts s k until the order J k . A numerically stable algorithm that produces such subspaces is the second order Arnoldi algorithm. It uses a modified Gram-Schmidt orthogonalization to ensure stability and is presented in [16,18]. Moment matching based on Krylov reduction, hereafter referred to as Krylov reduction, in its basic form produces reduced models with an order that is a multiple of the number of inputs and the number of shifts. For undamped systems with D = 0, b real vectors are added to the projection matrix for each shift so that n = bk. For damped systems, the Krylov vectors v are complex, and two vectors Re(v) and Im(v) are added for every input at each shift. This is possible because for a complex projection matrix with its complex conjugate at the complex conjugate shift s k is given by because of symmetry of the system matrices. Combining V k and V k and the transformation of basis shows that the complex projection matrix V can be represented by a purely real matrix, also see [12] on this topic. In this way, the moments are also matched at the complex conjugate shifts s k , but for every shift with Im(s k ) = 0 two vectors are added to V . For shifts placed in frequency domain on the real axis and with s = θ + i f /2π, the size of the reduced system is n = 2bk − 1 if one shift is placed at f = 0 and n = 2bk if f = 0 is not a shift. So far it has not been discussed how to choose the shifts at which the transfer functions are matched. The easiest way is to use linear spaced shifts in a frequency interval of interest. However, this choice can be inadequate for some systems since the dynamical behavior might be differently complex in different areas and, thus, needs coarser or finer sampling, accordingly. Therefore, experienced users often choose the shifts manually by trial and error as in [12]. Another possibility is to choose the shifts H 2 optimal as suggested in [14]. The goal of H 2 optimal reduction is the minimization of the H 2 -norm error of the linear, time invariant system with the H 2 -norm This system norm gives also an assessment to the time domain with for an output y, see [3,14]. An iterative algorithm that chooses the shifts in order to minimize the H 2 -norm error is presented in [14] and extended to second order systems in [21]. It is called IRKA. The algorithm is simple and often effective as shown in [8], but local convergence is proven only for state-space-symmetric systems only in [10], which does not apply to the system from Equation (2). In contrast to the basic Krylov reduction, the size of IRKA reduced models does not necessarily need to be a multiple of the number of inputs. This is because tangential directions that transform the MIMO problem to a single-inputsingle-output (SISO) problem in every iteration are used, see [2,7,14].

Singular value decomposition
The singular value decomposition (SVD) is one of the most powerful tools in linear algebra [1] used in a broad field of applications as, e.g., statistics and machine learning. It is defined in [1,6] as follows.

Definition 2 Given an arbitrary matrix
The column information of A is contained in U , whereas P contains the row information. The singular values stretch and clinch the mapping in the corresponding directions. Note that in general the SVD is not limited to real numbers but also applicable to every complex matrix with complex conjugate instead of transposed matrices. However, we stay in the real domain in the context of this paper.
The SVD can be used to approximate matrices as described in [1,6]. The matrix A is decomposed and with written as a sum of rank-one matrices. This sum can be truncated after n ≤ v terms leaving which is the optimal solution to the problem with the induced 2-norm ||A|| 2 = σ max (A). Hence, the singular values can be interpreted as importance measures that sort the basis vectors in terms of the importance to their contribution to approximate A. Just retaining the first n columns of U builds an orthogonal basis, which is the best approximation to the column span of A in terms of the 2-induced norm, see [1,6]. For the special case n = v, where all nonzero singular values are kept, the decomposition is called economy SVD.

Error measure
Any reduced model can only approximate the full model, i.e., there always occurs an approximation error. For the evaluation of different reduced models, a measure of this error is inevitable and, therefore, introduced in the following. The approximation error in frequency domain can be quantified with the relative error where

||H (s)|| F = trace H (s)H H (s) (23) is the Frobenius norm with the conjugate transpose H H (s). For SISO systems, H (s) is a scalar and
holds.

Influence of shift selection in Krylov reduction
This section addresses the influence of the shift selection on the size of the reduced model as well as on the approximation error. The limitations of existing methods are pointed out.

Linear shift selection
In mechanical engineering problems, often distinct frequency intervals are of special interest. Thus, the shifts for Krylov reduction are placed inside this interval and an equidistant choice seems natural. A finer shift sampling leads to a smaller approximation error but also to a higher order of the reduced system. Especially, if the distance between shifts is infinitely small, i.e., s → 0, the approximation becomes an exact representation over the frequency interval. However, the projection matrix is then assembled with infinitely many column vectors and fails the objective of a reduction. Therefore, the number of shifts is usually lower than the size of the original system n. It cannot be stated in general that more linear-spaced shifts lead to smaller approximation errors, i.e., the error does not decrease monotonically over the reduction order as shown in the following example. Fig. 1a is discretized with 945 volume elements.

Example 1 The clamped beam shown in
A system of 2700 ordinary differential equations describes its deformation due to external forces. Now, let f be a force acting at the beam's tip, and let the observed output y be the deflection at the tip. Figure 1b visualizes the transfer function of this SISO system. The peaks correspond to the excited natural eigenfrequencies. Higher eigenfrequencies are not visible due to damping. The beam is modeled with Rayleigh damping. To find a suitable reduced model, different numbers of shifts are now chosen and the approximation error is reviewed. In Fig. 2 the number of shifts increases iteratively. Therefore, we begin with two shifts at the boundaries of the interval [0Hz, 5000Hz] and add shifts between old shifts in the next iterations. The evaluated shifts are kept for the following iterations. From the fourth to the fifth iteration, there is no significant improvement and the two corresponding error plots are almost coincident. Here, the error is already very small because of the simplicity of the system. In contrast, in Fig. 3    Example 1 shows that if the shifts of two reduced models are complementary, it is possible that a reduced model generated with more shifts reveals larger approximation errors. This observation demonstrates, that some shifts give more information than others and make a greater contribution to the approximation of the transfer behavior over the frequency interval. However, the choice of the ideal shifts is far from trivial, and the importance of a shift cannot be determined beforehand. Tendentially, the maximal error drops for higher numbers of shifts, which confirms that sufficiently many shifts result in a reduced model that approximates the full one well. However, the order of the reduced model grows with the number of shifts.

H 2 optimal shift selection
A different approach to the shift selection is presented in [14]. It aims at choosing the shifts in a way that the H 2 -error of the reduced system is minimized. Therefore, a reduced system is generated with some initially chosen shifts. The eigenvalues of the reduced systems are computed and the shifts are moved to the mirror images of the poles of the reduced system. The eigenvalues are then recomputed and the shifts are moved to the mirror images of the new poles again. Previously computed shifts are not used for further iterations. This is done until convergence, i.e., until the eigenfrequencies do not change anymore.

Example 2
We apply the IRKA to the clamped beam from Example 1. Four different models are generated with 3, 5, 9, and 13 shifts, respectively. Fig. 4 shows the relative error for all four reduced models. The smaller two of them do not yield a good approximation over the whole frequency interval but only for lower frequencies. Only the larger two models with 9 and 13 shifts reveal a small error over the full frequency range.
The IRKA generates models with very small errors if sufficiently many shifts are evaluated. However, a comparison of Example 1 and Example 2 shows that IRKA requires similarly many shifts to yield small errors over the whole frequency space as the linear-spaced shifts strategy. Small reduced models are only appropriate for low frequencies as the shifts are placed at the mirror images of the first eigenfrequencies and, thus, in a low frequency range. For higher reduced orders, the error does not become noticeably smaller because the shifts are placed outside the frequency range of interest.

An error controlled automated shift selection method
In this section, an iterative algorithm for automated shift selection is presented. The idea behind it is to choose more shifts than necessary first and select the most relevant directions thereafter with SVD. This is based on the considerations from Sect. 2 that infinitely many shifts would lead to a vanishing approximation error but also to no reduction and that an optimal matrix approximation can be found with SVD. Since we cannot evaluate infinitely many shifts in finite time, we have to define a termination criterion that determines whether enough shifts have been evaluated. A realization of these thoughts can be found in Algorithm 1 where an automated shift selection is conducted until enough shifts are added to  form a large projection matrix and a further reduction of the dimension is achieved via SVD of this matrix.
For a system which is to be reduced, two inputs have to be passed initially. These provide an error bound that shall not be exceeded by the reduced model and a frequency interval in which this error is considered and in which the shifts are also placed. For mechanical systems, usually only closed intervals are relevant and, therefore, the systems must only match in this frequency range. With this interval, the first shifts to evaluate are set onto the boundary values of the interval. The original transfer function has to be computed once in order to calculate the approximation error in the following. The size of the reduced model is initially set to the size of the full model and the projection matrix is empty since we have not evaluated any shifts so far. The flag nChanged serves as termination criterion if it turns to false and is, thus, initialized with true. While true, Krylov subspace vectors are added to the matrix V full with an Arnoldi algorithm. However, no Gram-Schmidt orthogonalization is conducted. Instead, the matrix is decomposed with SVD in every iteration. The advantage of the SVD is the ordering of the columns of the new projection matrix U , which contains the left singular vectors sorted by their importance. The reduced system and its error can now be computed. Note that this is only possible because we restrict ourselves to matching moments of the order J = 1. For higher moments, a Gram-Schmidt orthogonalization has to be executed to ensure numerical stability, see [17]. If the error is already smaller than the predefined bound, an integer bisection algorithm finds the number n of columns that have to be kept in the projection matrix to not exceed the error limit. The algorithm stops as soon as the same order n results for two following iterations since it is then assumed that enough shifts are evaluated. If more shifts are needed, then shifts are added in between the shifts of the previous iteration, i.e., if in iteration ξ the shifts s ξ = are evaluated in the following iteration ξ + 1.

Example 3
We again consider the beam from previous examples. The goal is now to find a reduced model where the relative error never exceeds 10 −3 . Fig. 5 shows the error for all iterations of the presented algorithm with n = 3 to n = 33 and the error for the resulting reduced models where as many columns of V as possible are truncated with n = 10. The size of the reduced models is given in the legend. It can again be seen that the error drops for a higher number of shifts. The first time that the error does not exceeded the defined bound is for n = 17. For n = 9 instead, it is still larger at some frequencies. The difference between the models for n = 17 and n = 33 is comparatively small. Both result in an error which does not exceed 10 −3 , and after truncating 7 and 23 columns, respectively, the system is reduced to n = 10. This is the smallest order that is found in both iterations and does not violate the error bound and, therefore, the algorithm stops at this point.
The reduced model that results from the algorithm presented above exhibits an error that is relatively smooth, i.e., the error is not as small at the shifts as for the untruncated model but never as large as for a similar size model with basic linear shift selection. The SVD allows to transform the vectors of projection matrix so that the new vectors span the same space but are orthogonal. Additionally, an importance criterion for the new spatial directions is available with the singular values. In this way, a model that is sufficiently accurate for the given use case and small at the same time can be found. However, it cannot be said which shifts are particularly important and which shifts do not contribute heavily to the approximation of the transfer function.

Extension to MIMO systems
For MIMO systems, the number of shifts is even more important. This is because for every shift, b vectors, or even 2b vectors in the damped case, are added to the projection matrix. Thus, by using more shifts, the size of the reduced model grows quickly. Algorithm 1 is equally applicable to MIMO systems. Only the error norm considered must also map the behavior of different inputs to different outputs to one scalar. The Frobenius norm or the maximum over all inputs to all outputs can be used in this regard.

Example 4
A slightly more complex model than the clamped beam that is well suited as a benchmark for spatially distributed inputs and outputs is a clamped plate as visualized in Fig. 6a. Here, five input and output nodes are selected randomly, and their vertical movements are considered as inputs and outputs for the transfer function visualized in Fig. 6b. Fig. 7 shows the relative error for all iterations of Algorithm 1 if again an error bound of ε = 10 −3 is set. The first model with linear spaced shifts that does not exceed the bound is the one with n = 85. The SVD and the truncation lead to a model with = 47 for which the error is below the bound in the entire frequency range. Example 4 confirms that the effect of the SVD based truncation is even more valuable for MIMO systems since the order of the reduced system can be reduced from n = 85 to n = 47 in this way.

Frequency weighting
For some systems there might exist different frequency intervals where different errors are allowed. In mechanical systems, the lower frequency range is often most important as the amplitudes of low frequency oscillations are usually higher. Nevertheless, one might wish for a reasonably good approximation for other frequency ranges as well. With a slight adaption of Algorithm 1, this additional information can be included in the reduction process. Therefore, the entire frequency range is divided into intervals, and shifts are selected for every interval in Algorithm 2 according to its importance.

Algorithm 2
Error controlled automated shift selection with frequency weighting frequencyIntervals ← [I 1 , . . . , The main difference to Algorithm 1 is that the frequency intervals and error bounds are given as vectors. The shift selection and SVD are then conducted for each interval. In the end, the nonorthogonal matrix V full is assembled with all Krylov vectors and decomposed with economy SVD. In a final step it is checked how many columns of U can be truncated without violating an error bound in any interval.

Example 5
As an example for the frequency weighting, we separate the transfer function of the clamped beam from Example 1 into two intervals. The first interval covers the first half of the space of interest with [0Hz, 2500Hz] and the second covers the other half [2500Hz, 5000Hz]. Two cases are considered for demonstration purpose. It is assumed that the approximation quality has to be very good in the lower range, i.e., the error shall not exceed 10 −5 in case A, and it shall not exceed 10 −3 in case B. The error in the second half is bounded by 10 −1 only each time. Fig. 8 shows the error for both cases. It can be seen that the error in the right interval exceeds the error allowed in the left interval in both cases. However, the limits specified for the respective interval are not exceeded. The reduced models have a size of n = 8 and n = 10. If we compare the size with the model from Fig. 5, we see that a smaller model can be found if we restrict the interval where an error smaller than 10 −3 has to be fulfilled to the first half instead of the entire interval from 0 Hz to 5000 Hz.
Even if the error remains well below the specified limits for some intervals, frequency weighting can be applied by dividing the frequency space and providing different limits. The difference between limit and real error comes from the fact that the different intervals influence each other in their approximation quality. Nevertheless, it can be ensured that different errors are not exceeded in their intervals.

Low-error models for desired reduced orders
Up to now, we have always tried to find a model that is as small as possible with a predefined error bound. In this section, we consider the reverse problem, i.e., that a reduced order is predefined and a low-error model is to be found. Algorithm 3 addressing this problem is presented in the following. Algorithm 3 is similar to Algorithm 1, but the use of error and reduced order is interchanged. In the beginning, a shift vector that contains the minimum number of shifts to achieve the desired reduced order is initialized. Therefore, the shifts are again distributed equidistantly. The error for this model is calculated and then shifts are iteratively added. As long as the error changes and as long as it is larger than a defined tolerance, Krylov vectors are added for the new shifts to a matrix V full , which is then decomposed and truncated after the number of desired vectors. The error is recomputed for the new model and the termination criterion is checked.

Example 6
If we apply this algorithm to the beam model from Example 1, we obtain the results visualized in Fig. 9. A reduced order model with n = 10 is generated. The error is shown for the reduced model obtained from the truncated projection matrix where more than six shifts are evaluated. However, all models are truncated to the same size. The model with six shifts reveals the largest error. The other models are very similar and have a smaller maximal error. However, the model generated with 81 shifts does not yield the smallest error of all models. Example 6 demonstrates that it cannot be guaranteed that a model with more shifts always yields a smaller error everywhere. Nevertheless, the error for the resulting model, which is created after many shifts, reveals a much lower error than the model that would usually be used with the minimal required number of linear spaced shifts.

Application to an industrial numerical example
In this section, the presented algorithms are now applied to a complex numerical model of a helicopter airframe. The model represents a helicopter in an assembly line station where measurements are conducted to review quality. The FE model used is shown in Fig. 10.
There, also all considered inputs and outputs and the boundary nodes where the helicopter is mounted with springs are marked. The inputs and outputs are chosen at characteristic points for the vibration behavior of the airframe. Their designation is listed in the Algorithm 3 Automated shift selection for a given reduced order shifts ← shifts in the middle between old shifts end if end while  For the numerical simulation of the measurement with this model in a reasonable time, we need model order reduction techniques. Krylov reduction is well suited to generate a reduced model that is able to represent the measurements because it aims at approximating the transfer behavior from given inputs to given outputs. However, as stated above, the selection of shifts is not trivial. In Fig. 11 the relative approximation errors for different reduced models are compared. The largest error reveals the model reduced with IRKA, although the algorithm converged to a solution after eight iterations. The basic Krylov reduction with linearly spaced shifts leads to models where the order is a multiple of the number of inputs. Therefore, the model created with linear shift selection is slightly larger than the IRKA model. A model with three linear spaced shifts and n = 182 is the smallest model that can be achieved with linear shift selection, which does not exceed the error bound of ε = 10 −3 . By truncating as many columns as possible with Algorithm 1, a reduced model with n = 109 that does not exceed the same error bound can be found. For comparison, also a model with linear spaced shifts and n = 130 is plotted. This model is still larger than the SVD truncated one, but exhibits much larger errors. Thus, the presented method improves existing methods by either reducing the order for models with similar approximation errors or by reducing the error for models with similar reduction orders.
The results for the inverse approach, where n is constant and a low error is looked for, are presented in Fig. 12. Here, all models have a reduced order of n = 60 and the termination criterion is chosen to be ε = 10 −1 . It can be seen that the maximal error drops for more shifts. Even for 33 shifts, it does not get very low though. This indicates that a higher dimension than 60 is needed to approximate the range from 0 Hz to 50 Hz well. If the error tolerance is reduced, then further iterations are conducted. However, the full projection matrix then gets very large. In this case, for 65 shifts and 26 inputs, it would be V full = R 1434996×1690 for the following iteration that exceeds the memory of the 128GB RAM computer used for this work. The high memory consumption can be a limitation of the presented algorithm.
The reduced models can then be used to simulate a measurement. Therefore, the model is first decoupled with the modal matrix , which consists of the eigenvectors of the reduced system byM The decoupled system can be integrated efficiently. The results for time integration with the model generated with Algorithm 1 and n = 109 are presented in Fig. 13. There, a measurement is simulated where the airframe is excited with an impulse in lateral direction at the main gear box, and the lateral motion of different output nodes is observed. This simulation is only possible because of the use of MOR.

Numerical effort for different models
When evaluating the numerical effort, one can distinguish between offline costs for operations that only have to be done once and online costs that possibly repeat multiple times. In this framework, all operations to create a reduced model account to offline costs, whereas the evaluation of models is considered as online cost because often many evaluations are necessary. The separation of the two costs is partially blurred since, e.g., the evaluation of the transfer function is also needed to create the reduced models. In Fig. 14 the offline costs for the models presented in Fig. 11 are contrasted. The first model is created by selecting equidistant shifts. To compare this approach to the SVD based approach, shifts are iteratively increased starting with two. After every iteration, the error of the model is checked, and the algorithm terminates as soon as the error bound of ε = 10 −3 is not exceeded anymore. More specific, models are build with n s k = 2, 3, 4, . . . shifts until the smallest reduced model with linear shift selection that suits the given error bound is found. This procedure results in a model with n = 182 DOF. Of course, the computation time would be much less if we knew the number of necessary shifts beforehand. However, there is no approach to guess the shift number in advance and, thus, an iterative increase is necessary. The IRKA algorithm is by far the most time-consuming algorithm of the three. With more than 57 hours of CPU time, the generation of the reduced model takes more than 27 times longer than the iterative linear shift selection. An Intel(R) Xeon(R) CPU E5-1650 v3 @ 3.50GHz processor is used for this work. The additional expense would be justifiable if the algorithm led to models with significantly smaller errors or smaller models for similar sized errors. Since this is not the case here, the IRKA is not effective for the presented problem. The fastest algorithm of the three is the greedy shift selection with SVD based truncation that was presented in Algorithm 1. For an error bound of ε = 10 −3 , it takes 1.35 CPU hours and leads to the model shown in Fig. 11 with n = 109.
The high offline costs for the preparation of adequate models are legitimized if one compares the costs for evaluations of the full model with evaluations of the reduced model. The computation of the transfer function from all inputs to all outputs in the frequency range [0Hz, 50Hz] takes more than 13 hours of CPU time. By contrast, the evaluation of the transfer function of a reduced model with n = 200, which is rather large in the scope of this paper, takes about 40 s. This corresponds to a speedup by a factor greater than 1000. A time evaluation of the full model is not even possible with the used computer. The computation costs for the evaluation of different size reduced models are compared in Fig. 15. There, a time integration for a measurement scenario with impulse excitation at the main gear box and measurements at all other points is executed for 1 s. The time consumption for the evaluation of the transfer function, for the time integration, and for the calculation of the reduced system matrices is shown. Apparently, the costs rise for all three operations with the order n of the reduced model. Therefore, a model as small as possible is desirable and higher preprocessing costs might be acceptable, especially if many system evaluations are necessary. In this case, a model of the size n = 109, which was above generated with SVD based truncation and showed small approximation errors, takes two fifths of the time of the time integration for a model with n = 182. Because of that, the higher costs for the presented algorithm amortize after few evaluations. An application scenario where many evaluations are necessary would be, e.g., the determination of a parameter by comparison of simulated measurement at different parameter points and real measurement. In this context, parametric model order reduction can be applied so that the reduction algorithm has to be applied for only some parameter samples in the parameter space to obtain a reduced model that is valid for the entire parameter space, see [4].

Conclusion
In this paper, a new approach was presented to choose the shifts in moment matching based on Krylov subspaces. This approach chooses the shifts automatically so that a given error bound is not exceeded or so that a given reduction order is obtained. The advantage of this procedure is that no user input regarding the position of the shifts is necessary and that small models can be found with good approximation quality. To achieve this, in this contribution Krylov subspaces were approximated with SVD, and only the most important directions were kept for the reduced model. Two algorithms were presented where either a low-rank model for a given error bound or a small-error model for given reduced order is created. Compared to the state-of-the-art, the presented approach is able to achieve good results for academic examples as well as for an industrial example of a helicopter cell, both in terms of computation time and approximation error.