Model reduction for second-order systems with inhomogeneous initial conditions

In this paper, we consider the problem of finding surrogate models for large-scale second-order linear time-invariant systems with inhomogeneous initial conditions. For this class of systems, the superposition principle allows us to decompose the system behavior into three independent components. The first behavior corresponds to the transfer between the input and output having zero initial conditions. In contrast, the other two correspond to the transfer between the initial position and the initial velocity conditions having zero input, respectively. Based on this superposition of systems, our goal is to propose model reduction schemes allowing to preserve the second-order structure in the surrogate models. To this aim, we introduce tailored second-order Gramians for each system component and compute them numerically, solving Lyapunov equations. As a consequence, two methodologies are proposed. The first one consists in reducing each of the components independently using a suitable balanced truncation procedure. The sum of these reduced systems provides an approximation of the original system. This methodology allows flexibility on the order of the reduced-order model. The second proposed methodology consists in extracting the dominant subspaces from the sum of Gramians to build the projection matrices leading to a surrogate model. Additionally, we discuss error bounds for the overall output approximation. Finally, the proposed methods are illustrated by means of benchmark problems.


Introduction
Second-order dynamical systems arise in many engineering applications, e.g., electrical circuits, structural dynamics, and vibration analysis. In many setups, these systems are modeled by partial differential equations having secondorder time-derivatives. In order to compute the numerical simulations, spatial discretizations are needed, leading to high fidelity models. However, those high fidelity models may present a high number of degrees of freedom, which are not suitable for numerical computations. Consequently, model order reduction techniques are used to construct reduced-order models, that approximate the behavior of the original system.
Most reduction techniques assume that the considered systems have zero initial conditions. Consequently, these methods fail in approximating systems if they have inhomogeneous initial conditions. Additionally, the corresponding error estimators of these methods are not applicable in this case. This work is dedicated to finding surrogate models for second-order systems with inhomogeneous initial conditions while preserving the system structure.
In the literature, there exist several reduction methods dedicated to systems with homogeneous initial conditions. Examples are singular value based approaches such as balanced truncation [9,20,31] and Hankel norm approximations [15]. Additionally, there exist Krylov based methods such as the iterative rational Krylov algorithm (IRKA) [9,13,16], as well as, data driven methods such as the Loewner framework [19]. In this work, we consider secondorder continuous-time dynamical systems governed by the system of differential equations where M, D, K ∈ R n×n , B ∈ R n×m , C ∈ R p×n , x(t) ∈ R n , u(t) ∈ R m and y(t) ∈ R p . We assume that the position and velocity initial condition are not known a priori. However, they are assumed to lie in two known subspaces X 0 := span {X 0 } and V 0 := span {V 0 }, respectively, with X 0 ∈ R n×nx 0 and V 0 ∈ R n×nv 0 . Hence, the initial conditions can be expressed as Our main goal in this work is to find low dimensional surrogate models for the system (1) with inhomogeneous conditions (2) preserving the second-order system structure. By structure-preserving, we mean to determine Petrov-Galerkin projection matrices W, V ∈ R n×r leading to a reduced second-order system Mẍ r (t) +Dẋ r (t) +Kx r (t) =Bu(t), x r (0) =X 0ẑ0 ,ẋ r (0) =V 0ŵ0 , withM = W T MV,D = W T DV,K = W T KV,B = W T B,Ĉ = CV,X 0 = W T X 0 ,V 0 = W T V 0 , and x r (t) ∈ R r . In the literature, there exist several methods enabling model order reduction preserving the second-order structure [11,14]. These techniques range from balanced truncation as well as balancing related model order reduction [12,23,29] to moment matching approximations based on the Krylov subspace method [4,26]. The recent work [25] provided an extensive comparison among common secondorder model reduction methods applied to a large-scale mechanical fishtail model. Additionally, [5] proposed interpolation based methods for systems possessing very general dynamical structures. More recently, the authors in [7] propose a new philosophy enabling to find the dominant reachability and observability subspaces enabling very accurate reduced-order models preserving the structure. Moreover, an extension of the Loewner framework was proposed in [8] for the class of Rayleigh damped systems and in [28] for general structured systems.
Up to our knowledge, there is no dedicated work on system theoretical model reduction of second-order systems with inhomogeneous initial conditions. For the class of firstorder systems with inhomogeneous conditions, we briefly review four proposed approaches from the literature. In [2], the authors proposed to shift the state by the initial condition x 0 , e.g. the new state is given asx(t) := x(t) − x 0 . That way, the initial condition is included in the input and output equation and therefore considered in the reduction process. This method, however, is not straightforward applicable to second-order systems if we have a velocity initial condition and want to preserve the second-order structure.
In [17] the input Bu(t) is extended by the initial condition space X 0 . More detailed, a new input matrixB := [B X 0 ] and a new input [u(t) z 0 ] T are defined such that the initial condition is taken into account applying reduction methods. As in the previous method, this approach is not feasible if we consider velocity initial conditions in the second-order case.
In [3], the authors' strategy is to decompose the system into a zero initial condition system and a system with initial conditions but no input. The sum of the two corresponding outputs provides the original output. This superposition is used to reduce these two systems, separately. Extensions of the proposed methodology for the class of bilinear systems is proposed in [10] and [22] based on different splittings.
A recent approach [27] proposes a new balanced truncation procedure based on the shift transformation on the state. This transformation is depending on design parameters allowing some flexibility and enabling the generalization of the methodologies proposed in [17] and [3]. Additionally, those parameters can be optimized, leading to accurate reduced-order models.
In this paper, the superposition ideas in [3] are extended to the class of second-order systems. For the later class, we show that, due to the superposition principle, the original system can be decomposed into three subsystems. The first subsystem corresponds to the map between the input u(t) and the output while the initial conditions are set to zero. Additionally, the second subsystem is the output resulting from the position initial condition x(0) and the third one corresponds to the output obtained using the velocity initial conditionẋ(0). Hence, we analyze the three corresponding subsystems separately. For the systems which result from the initial conditions, we develop balancing based reduction techniques based on tailored Gramians that are introduced in this paper.
Here, two model reduction schemes are proposed. The first one consists in reducing each of the components independently using a suitable balanced truncation procedure. Hence, the sum of these reduced systems provides an approximation of the original system. As a consequence, an advantage of this approach is that the reduced dimensions and therefore the accuracies can be chosen flexibly. The second proposed methodology consists in extracting the dominant subspaces from the sum of Gramians to build the projection matrices leading to one surrogate model.
The rest of the paper is organized as follows. In Section 2, we present balanced truncation for first and secondorder systems. Afterwards, in Section 3, we deduce a superposition of the second-order system (1). In Section 4, tailored Gramians for inhomogeneous second-order systems are derived. Based on these Gramians, two model reduction schemes are proposed in Section 5. Finally, Section 6 provides the resulting error estimation and in Section 7, the methodologies are illustrated in two numerical examples.

Balanced truncation
In this section, we briefly present balanced truncation method for first-order and second-order systems having zero initial conditions.

First-order systems
We consider the first-order dynamical system with zero initial conditions We assume that the system is asymptotically stable, i.e. all eigenvalues λ of the matrix pencil A − λE fulfill Re(λ) < 0.
The goal of balanced truncation is to find a reduced-order model, that approximates the input-output behavior of (4). The corresponding controllability and the transformed observability Gramian are defined as The Gramian P and the transformed Gramian Q can be computed by solving the Lyapunov equations Small singular values of P and E T QE correspond to states, that are difficult to reach and to observe. In order to truncate small singular values of P and E T QE, simultaneously, we transform the system such that the transformed Grami-ansP,Q fulfillP =Q = Σ = diag (σ 1 , . . . , σ n ), where σ 1 , . . . , σ n are called Hankel singular values. This transformation is called balancing. Afterwards, we truncate the n − r smallest Hankel singular values σ r+1 , . . . , σ n , r ≪ n. Therefore, we consider the low-rank factors RR T = P and SS T = Q and compute the singular value decomposition S T ER = UΣX T . The resulting projection matrices that do both, balance and truncate, are where Σ r := diag (σ 1 , . . . , σ r ) and U r and X r include the r leading columns of U and X. The balanced and truncated system is then given bẏ since W T EV = I r . For more details about standard balanced truncation, see [6,20].

Second-order systems
Balanced truncation for second-order systems (1) with zero initial conditions is presented in [12]. By applying the Laplace transform to the second-order system system (1) with zero initial conditions provides the following transfer function: First, we define the input to state and the state to output mapping that result from the transfer function.
Definition 2.2. The input to state mapping R SO and the state to output mapping S SO of the second-order system (1) with zero initial conditions are The corresponding second-order controllability Gramian P SO and observability Gramian Q SO are defined by In order to reduce the second-order system (1) with zero initial conditions, we transform it to a first-order system (4) by setting The first-order system is then equivalent to the second-order system and the corresponding transfer function is given by Note that there exist several first-order representation that are equivalent to system (1). We determine the secondorder controllability Gramian of system (1) with zero initial conditions by considering the Gramian of the first-order system (4) as described in the following proposition.
Proposition 2.1. The second-order controllability Gramian P SO of System (1) with zero initial conditions is equal to the upper left block P 1 of the first-order controllability Gramian Proof. Applying the Schur complement provides that P 1 is given by The matrix P SO is called position controllability Gramian. We observe that P SO encodes the important subspaces of the map between the input to state of the homogeneous secondorder system (1). Hence, P SO spans the controllability space and is used to apply balanced truncation in the second-order case.
The same argumentation is used to extract the state to output mapping space from the first-order observability As in the first-order case we use the low-rank factors R 1 and S 3 with P 1 = R 1 R T 1 and Q 3 = S 3 S T 3 and compute the singular value decomposition S T 3 R 1 = UΣX T . The resulting balancing and truncating projection matrices are where Σ r is the diagonal matrix containing the r largest singular values of Σ. Moreover, U r and X r include the r leading columns of U and X. Projecting by W and V provides the reduced system (3) with zero initial conditions.

Superposition principle for second-order systems
This section aims at decomposing the original system behavior of the second-order system (1) in simpler subsystems. This system decomposition will be the inspiration of the proposed model reduction schemes. By applying the Laplace transform to equation (1a) we obtain where X is the Laplace transform of x and U the Laplace transform of u. Hence, it holds that Applying the Laplace transform to equation (1b) and defining Y as the Laplace transform of y provides We observe that the output is a superposition of the input to output mapping, the position initial condition to output mapping and the velocity initial condition to output mapping. As a consequence, the global input-output behavior is given by Up to now, we saw that three independent transfer functions characterize the inhomogeneous behavior of the second-order realization (1). The transfer function H SO (s) = CΛ(s)B corresponds to the input and output map without initial conditions. Hence, it is associated with the following realization x SO (0) = 0,ẋ SO (0) = 0.
The transfer function H x0 (s) = CΛ(s)(D + sM)X 0 corresponds to the transfer between the initial position condition and the output. Hence, the following realization is associated to it: x x0 (0) = X 0 z 0 ,ẋ x0 (0) = 0.
Finally, we write the realization for H v0 (s) = CΛ(s)MV 0 . This transfer function corresponds to the transfer between the initial velocity condition and the output. The following realization is associated to it: x v0 (0) = 0,ẋ v0 (0) = V 0 w 0 .
To summarize, we have seen that the output of the inhomogeneous second-order system in (1) can be decomposed as

Gramians of inhomogeneous second-order systems
In order to derive the proposed model reduction schemes, we analyze separately the three subsystems and we introduce tailored Gramians for each one of them. Notice that subsystem (7) corresponds to a second-order realization with homogeneous initial conditions. Hence, the controllability and observability Gramians P SO and Q SO presented in Definition 2.2 can be used to characterize its dominant subspaces.
However, subsystems (8) and (9) have a different structure, and hence, tailored Gramians are required. In Section 4.1 and 4.2, we propose tailored Gramians for these subsystems. Afterwards, in Section 5, we propose two different MOR schemes based on these Gramians.

Gramians of H x 0
Considering the transfer function H x0 (s) of system (8) more detailed shows that the input to state mapping differs from the structure in Definition 2.2. The state to output mapping, however, is the same. Hence, we define the input to state mapping and the corresponding second-order Gramian.
Definition 4.1. The input to state mapping R x0 and the corresponding controllability Gramian P x0 of the secondorder system (8) are Proof. Applying the Schur complement to Γ(iw) provides that its upper left block is −((iw) 2 M + iwD+ K) −1 (iwM + D) and hence it holds that Proposition 4.1 shows that the second-order controllability Gramian P x0 of System (8) is given by the upper left part P 1 of the controllability Gramian P of the first-order system (4) with B := X 0 0 . Moreover, the second-order observability Gramian Q x0 is equal to Q SO since the state to output mapping S x0 (s) := C s 2 M + sD + K −1 that is used to derive the observability Gramian is equal to the state to output mapping S SO (s) from Definition 2.2 for the homogeneous system case.

Gramians of H v 0
In order to apply balanced truncation to system (9), we define the corresponding input to state mapping. As in the previous section, the state to output mapping is the same as for system (7).
Definition 4.2. The input to state mapping R v0 and the corresponding controllability Gramian P v0 of the secondorder system (9) are We note that the input to state mapping R v0 (s) and hence the second-order controllability Gramian P v0 are of the same structure as in the homogeneous case, presented in Definition 2.2.
Applying the Schur complement to Γ(iw) provides that its upper right part is −((iw) 2 M − iwD + K) −1 and hence Again the second-order observability Gramian Q v0 is equal to the one of the homogeneous second-order system Q SO since their state to output mappings coincide.

Model reduction schemes
In this section, we present two model reduction schemes for the class of systems in (1). The procedures use the tailored Gramians presented in Section 4 and construct second-order reduced-order models via balanced truncation as presented in Subsection 2.2.

Method 1: Reducing each subsystem
The first method we propose utilizes the superposition properties to reduced the subsystems presented in Section 3 separately based on the Gramians presented in Definition 2.2 and in Section 4.
For the homogeneous subsystem (7), we aim to apply the reduction procedure from Subsection 2.2 using the Gramians from Definition 2.2 to derive a reduced-order system with the transfer function where W SO and V SO are the corresponding projection matrices given in (6). For the subsystem (8) describing the system behavior that results from the state initial condition, we build the corresponding balanced truncation projection matrices W x0 and V x0 as in Equation (6) based on the Gramians presented in Subsection 4.1. We reduce system (8) accordingly and obtain the reduced position initial condition transfer function Applying second-order balanced truncation to system (9) using the second-order Gramians from Subsection 4.2 provides the projection matrices W v0 and V v0 from Equation (6). Reducing the system accordingly generates the corre-sponding reduced transfer function that describes the velocity initial condition to output behavior. Summarizing, we apply balanced truncation to the three systems to generate the corresponding reduced transfer functionsĤ SO ,Ĥ x0 andĤ v0 associated with the outputŝ y SO (t),ŷ x0 (t), andŷ v0 (t), respectively, such that the overall behaviorŷ =ŷ SO (t) +ŷ x0 (t) +ŷ v0 (t) approximates the original output y(t). The detailed reduction procedure for each subsystem is given in Algorithm 1.

Method 2: Combined Gramians
We have discussed the approach where we use separated projections for each subsystem. However, for some applications it might be advantageous to have only one projection that reduces the original system including the initial conditions at once. That means that we need to determine a projection based on a controllability space which corresponds to the input and the initial conditions. This controllability space is spanned by the columns of the sum of the controllability Gramians introduced in the previous sections where R T SO , R T x0 , R T v0 and S SO are the corresponding lowrank factors of P T SO , P T x0 , P T v0 and Q SO . Applying balanced truncation for second-order systems based on the low-rank factors of the combined Gramains P c and Q c from equation (10) results in a reduced-order system that takes into account the input to state and the initial conditions to state mappings.
Another approach that results in the same controllability Gramian P c and therefore the same reduced-order system would be a modification of the method presented in [17] for second-order systems. Therefore, we consider the homogeneous first order system (4) with the input matrix Proposition 5.1. The position controllability Gramian of the first order system (4) with B as in equation (11) is equal to the Gramian P c in equation (10).
Proof. The position controllability Gramian of system (1) is described by the upper left part P 1 of Applying the Schur complement to Γ(iw) provides that the upper left block is −((iw) 2 M + iwD + K) −1 (iwM + D) and the upper right block is −((iw) 2 M+iwD+K) −1 . It follows that The final method is presented in Algorithm 2 and generates the reduced transfer function The advantage of this approach is the fact that we obtain only one second-order reduced-order model approximating the behavior of the original system. The disadvantage of this method is the inflexibility of the different controllability space dimensions. It follows, that the combined Gramian leads possibly to reduced dimensions significantly larger than the dimensions of the separately reduced systems to reach the same approximation quality.

Error bounds
In this section we develop a posteriori error bounds for the methods of this article. Therefore, we use the fact that Algorithm 2 BT method for inhomogeneous second-order systems by combined Gramians.
Require: The original matrices M, K, D, B, C, X 0 , V 0 and the order r. Ensure: The reduced matricesM,K,D,B,Ĉ,X 0 ,V 0 .
1: Build the input matrix 2: Compute low factors of Gramians P ≈ RR T from (10) and Q ≈ SS T from Definition 2.2. 3: Perform the SVD of S T R, and decompose as Firstly, we use the above inequality to find a posteriori error bounds for the reduction scheme presented in Subsection 5.1. As a consequence, a possible error bound for the reduced subsystems approximation is Using the H 2 norm has the advantage of less computational costs. However, one needs to have u ∈ H ∞ , which applies some restrictions to the family of inputs u because u ∈ H ∞ is a stronger condition than u ∈ H 2 . We where P andP are the controllability Gramians of the original first-order system and the reduced first-order system. The cross GramianP solves the Sylvester equation The controllability Gramians P of the full system needs to be computed anyway to apply balanced truncation and the reduced GramiansP andP are cheap to compute. The error estimation for the combined Gramian and the resulting reduced-order system is equal to the error estimate above, where the projection matrices that lead to the reduced system are the same for the three subsystems. On the other hand one can evaluate for the matrix B as in (11) and u c := . This estimation can be computed using the equation (12).

Numerical results
In this section, we illustrate the procedure presented in this article using two different examples. The first example is a vibrational model of a building and the second one a mass spring damper system. We evaluate and compare for each example three reduced systems. We obtain the first one by applying balanced truncation to the full system (1) that does not consider the initial conditions. Hence, the projection matrices result from the evaluation of the Gramians of the homogeneous system (7). The second reduced system is obtained by reducing the three subsystems, separately, as presented in this article. The third method uses the combined Gramian presented in Section 5.2 to obtain the reduced-order model. For all reduced systems, we evaluate the output behavior and the corresponding output error.
The computations were done on a computer with 4 Intel ® Core™i5-4690 CPUs running at 3.5 GHz. The experiments use MATLAB ® R2017a. In the second example, the Lyapunov equations are solved using the methods from the M-M.E.S.S. toolbox [24].
We will refer to the original system (1) as FOM, in the following, and to the reduced system generated by standard balanced truncation that considers homogeneous systems, i.e, by applying second-order balanced truncation as described in Subsection 2.2 by ROM HOM. The reduced system approximation that is obtained by applying method 1 introduced in Subsection 5.1 is referred to as ROM SPL and the reduced system that is generated by applying method 2 introduced in Subsection 5.2 as ROM COM.

Building example
In this section we consider the building example from page 17 of the technical report [1]. The data are available in [21]. The dimension of the matrices are n = 24, m = p = 1. For the projection matrix W SO that results from the balanced truncation procedure for the homogeneous secondorder system (7) we consider the singular value decomposition Assume that rank(W) = ℓ. The position and velocity initial condition are the (ℓ + 1)-st column of U: In this example, the separately reduced systems and the combined reduced system are truncated with a reduced dimension r = 10. Figure 2a shows the output behavior of the original system and the reduced ones for an input u(t) = 0.2 · e −t . We observe that the original output behavior that is depicted in green is well approximated by the separately reduced system ROM SPL that is depicted by the blue, dashed line. The reduced system ROM COM using the combined Gramian (depicted by the orange colored, dashed line) provides a proper approximation of the original output as well. Additionally, we see that the reduced output of the reduced system ROM HOM, which is depicted in red, fails in approximating the original system's transient behavior. Figure 2b depicts the errors and their ℓ 2 -norms. The light blue line with markers depicts the error of the separately reduced system ROM SPL and the dashed, brown colored line the error of the reduced system ROM COM using the combined Gramian. The reduced system ROM HOM leads to the error depicted by the dashed, orange colored line. We observe, that the separately reduced system and the reduced system that uses the combined Gramain lead to errors that are  significantly smaller than the error corresponding to the reduced system ROM HOM. Additionally, we evaluate the actual ℓ 2 -norm error. Therefore, we plot the integral t 0 y(t) −ŷ(t) 2 dt (13) that converges to the ℓ 2 -norm of the error. The dark blue, dashed line with markers is the integral (13) converging to the actual ℓ 2 -norm error of the separately reduced system ROM SPL. The error bound from Section 6 provides a value of 3.2740 · 10 −5 (depicted by the black line). We see that this error estimator provides a proper upper bound of the actual ℓ 2 -norm error. The green line with markers provides the integral (13) corresponding to the combined Gramian reduced system ROM COM and its error estimation 1.5469 · 10 −4 is depicted by the dashed, black line. The red line shows the integral (13) of the reduced system ROM HOM. It confirms again, that this method fails for this example.

Mass spring damper example
The mass spring damper model we consider in this section is presented in [30]. More detailed background can be found in [18]. We choose the model of dimensions n = 2000, m = p = 1. The input is the external forcing on the n-th mass and the output observes the n-th mass.
The initial conditions are set to be the last and the first unit vector In this example, we truncate the systems with a tolerance of 10 −4 , i.e. all Hankel singular values smaller than 10 −4 ·σ 1 are truncated. That way, we obtain reduced systems of dimensions 147, 180, 98 of the three systems resulting from the superposition method and the dimension 157 for the system reduced using the combined Gramians. Figure 3a shows the output behavior of the systems for the input u(t) = 0.2 · e −t . The output behavior of the original system is depicted in green. The blue, dashed line displays the output composed by the separately reduced systems ROM SPL and the orange colored, dashed line the reduced system ROM COM using the combined Gramian. The reduced output resulting from the reduced system ROM HOM is depicted in red. We observe that all outputs approximate the original system behavior. although ROM HOM shows oscillations of slightly higher magnitude than the FOM for some time.
The output errors and their ℓ 2 -norms are illustrated in Figure 3b. The light blue line with markers, the brown colored, dashed line and the orange colored, dashed line show the error of the separately reduced outputs, the output corresponding to the combined Gramian and the output resulting form the reduced system ROM HOM, respectively. We observe again that the separately reduced system ROM SPL and the reduced system ROM COM using the combined Gramian lead to smaller errors. Additionally, we evaluate the actual ℓ 2 -norm error and plot the integral (13) that converges to the ℓ 2 -norm of the error. The dark blue, dashed line with markers shows the integral (13) for the separately reduced system ROM SPL and the green one the integral for the reduced system ROM COM using the combined Gramian. The error estimator from Section 6 provides ℓ 2 error estimation values of 7.5490 · 10 −3 and 3.1922 · 10 −2 for this example. It is depicted in Figure 3b by the black and black, dashed lines. We observe that the error estimations are conservative. The integral (13) of the reduced system ROM HOM is depicted in red. It converges to a ℓ 2 error that is larger than for the first two reduction methods.

Conclusion
We have proposed two approaches for constructing a reduction of second-order linear time-invariant systems with inhomogeneous initial conditions. First, we have used a superposition of the output into the input to output mapping, the state initial condition to output mapping and the velocity initial condition to output mapping. The three subsystems have been reduced, separately, such that the original system can be approximated well. Afterward, a combined Gramian has been used to derive projection matrices that reduce the system, including the initial conditions, all at once. For those reduction processes we have suggested new Gramians for inhomogeneous second-order systems.