Stability preservation in Galerkin-type projection-based model order reduction

We consider linear dynamical systems consisting of ordinary differential equations with high dimensionality. The aim of model order reduction is to construct an approximating system of a much lower dimension. Therein, the reduced system may be unstable, even though the original system is asymptotically stable. We focus on projection-based model order reduction of Galerkin-type. A transformation of the original system guarantees an asymptotically stable reduced system. This transformation requires the numerical solution of a high-dimensional Lyapunov equation. We specify an approximation of the solution, which allows for an efficient iterative treatment of the Lyapunov equation under a certain assumption. Furthermore, we generalise this strategy to preserve the asymptotic stability of stationary solutions in model order reduction of nonlinear dynamical systems. Numerical results for high-dimensional examples confirm the computational feasibility of the stability-preserving approach.


Introduction
The mathematical modelling of problems from science and engineering often yields dynamical systems. The increasing complexity of industrial applications causes high-dimensional systems by electronic design automation, for example. Thus a numerical simulation of the model may become too costly. Methods of model order reduction (MOR) are required to decrease the dimensionality of the dynamical systems, see [1,4,24].
We consider linear implicit systems of ordinary differential equations (ODEs), which are asymptotically stable. Projection-based MOR yields linear ODEs of a lower dimensionality. However, the reduced system may be unstable and thus useless. Firstly, some solutions become unbounded in the time domain. Secondly, error bounds, which follow from the transfer functions in the frequency domain, are not valid any more. Hence stability-preserving MOR methods are necessary to guarantee adequate reduced systems.
The method of balanced truncation, see [11], always yields stable reduced systems, while the computational effort is often relatively large. Krylov subspace techniques, see [9], are cheaper, whereas stability can easily be lost. A stabilitypreservation of a Krylov subspace approach is given by special assumptions and methods in [13]. The stability property can also be satisfied by a post-processing using the poles of the transfer function, see [2].
We examine projection-based MOR of Galerkin-type, where each scheme is defined by a single orthogonal projection matrix. Prominent methods are the onesided Arnoldi algorithm and the proper orthogonal decomposition (POD), for example. Prajna [21] introduced an approach to guarantee the stability in such an MOR of a nonlinear system of ODEs by a basis transformation in both the state space and the image space. This technique was also applied to a stochastic Galerkin projection in [23]. Castañé Selga et al. [7] investigated a stabilisation of linear systems of ODEs by a transformation in the image space only. The main effort consists in solving a single Lyapunov equation, whose efficient numerical solution is critical. The high dimensionality excludes direct methods and thus approximate methods have to be used. The alternating direction implicit (ADI) algorithm, see [15,20], requires an input matrix in the form of a symmetric low-rank factorisation.
We perform a specific ansatz for an approximate solution of the Lyapunov equation. Our strategy yields an alternative Lyapunov equation, where a symmetric factorisation can be chosen as input matrix. Its rank depends on the signs of eigenvalues in a symmetric part. Consequently, the ADI algorithm is applicable with a small number of iteration steps in the case of low ranks. Moreover, our method guarantees a non-singular and often well-conditioned mass matrix in the reduced system. The technique can also be employed to preserve the asymptotic stability of stationary solutions in an MOR of nonlinear ODEs.
The paper is organised as follows. We describe the class of MOR methods as well as the stability-preserving transformation in Section 2. Our numerical method is derived and analysed in Section 3. We demonstrate the applicability of this approach to nonlinear dynamical systems in Section 4. Finally, Section 5 presents numerical results of illustrative examples.

Projection-based reduction and stability
The projection-based MOR and the stability problem are defined in this section.

Linear dynamical systems and stability
Let a linear dynamical system be given in the form of implicit ODEs with state variables x : [0, ∞) → Ê n , inputs u : [0, ∞) → Ê n in and outputs y : [0, ∞) → Ê nout . The system includes constant matrices A, E ∈ Ê n×n , B ∈ Ê n×n in and C ∈ Ê nout×n . We assume that the mass matrix E is non-singular.
In the frequency domain, a transfer function completely describes the inputoutput behaviour of the system (1), see [1]. This transfer function H : \Σ → nout×n in reads as The mapping (2) is a rational function with a finite set of poles Σ ⊂ . The magnitude of a transfer function can be characterised by norms in Hardy spaces. The H 2 -norm is defined by, see [26, p. 92], with i = √ −1 and the Frobenius (matrix) norm · F provided that the integral exists.

R. Pulch
The stability issues of the system (1) are independent of the definition of inputs and outputs. To discuss the stability, we recall some general properties of matrices.
Definition 1 Let A ∈ Ê n×n and λ 1 , . . . , λ n ∈ be its eigenvalues. The spectral abscissa of the matrix A is the real number The definition of dissipativity is in agreement to [10, p. 62]. A dissipative matrix is also stable. Vice versa, a stable matrix is not dissipative in general. The linear dynamical system (1) is asymptotically stable, if and only if each eigenvalue λ ∈ satisfying det (λE − A) = 0 (4) has a strictly negative real part, see [5, p. 376]. Equivalently, E −1 A is a stable matrix. The asymptotic stability guarantees the existence of the integral in (3). If the eigenvalues of the problem (4) all have a non-positive real part and a real part zero appears, then Lyapunov stability may still be satisfied. However, we consider this case also as a loss of stability, because the asymptotic stability is not valid any more. Likewise, the linear dynamical system (1) is called dissipative, if and only if E −1 A is a dissipative matrix.

Projection-based model order reduction
We assume that the linear dynamical system (1) exhibits a huge dimensionality n. Thus the involved matrices A and E typically are sparse. The purpose of MOR is to decrease the complexity. An alternative linear dynamical system has to be constructed with state variablesx : [0, ∞) → Ê r and the matrices A,Ē ∈ Ê r×r ,B ∈ Ê r×n in ,C ∈ Ê nout×r , where the dimension r is much smaller than n. Again we assume thatĒ is non-singular. Initial valuesx(0) =x 0 are approximated from the initial values x(0) = x 0 . Nevertheless, the output of (5) should be a good approximation to the output of (1), i.e.,ȳ(t) ≈ y(t) for all relevant times. The system (5) is called the reduced-order model (ROM) of the full-order model (FOM) given by (1). The linear dynamical system (5) has its own transfer functionH : \Σ → nout×n in of the form (2). If both the original system (1) and the reduced system (5) are asymptotically stable, then error bounds are available in the case of x 0 = 0 andx 0 = 0. It holds that, see [3, p. 496 the H 2 -norm (3), the maximum (vector) norm · ∞ and the Euclidean (vector) norm · 2 .
In projection-based MOR, see [1], two projection matrices V, W ∈ Ê n×r of full rank are specified. We obtain the matrices of the ROM (5) bȳ The orthogonality V T V = I r and sometimes the biorthogonality W T V = I r are supposed with the identity matrix I r ∈ Ê r×r . Often the projection matrices result from the determination of subspaces, i.e., V = span(V ) ⊂ Ê n and W = span(W ) ⊂ Ê n .
On the one hand, the original state variables are approximated within the space V by x ≈ Vx. On the other hand, the residual is kept small by the requirement s(t) ⊥ W and thus W T s(t) = 0 for all t.
It holds that W = V in a Galerkin-type projection, where just an appropriate projection matrix V has to be identified. Important examples are the one-sided Arnoldi method, see [9], and proper orthogonal decomposition (POD), see [1]. However, the investigation of optimal choices for the matrix V is not within the scope of this paper. We consider an arbitrary matrix V satisfying V T V = I r .
In each MOR approach, the reduced system (5) is often useless if it is not at least Lyapunov stable. Many Krylov subspace methods or moment matching techniques do not guarantee a stable ROM.

Stability preservation by transformations
We apply the following well-known property for Galerkin-type projection-based MOR, see [7,21].
Theorem 1 In the linear dynamical system (1), let E be symmetric positive definite and A be dissipative. It follows that the reduced system (5) is asymptotically stable for the matrices (8) with a full-rank matrix V and W = V . Proof: Let B ≡ 0 in (1) without loss of generality. Since E is symmetric positive definite and V has full rank, the matrixĒ = V T EV is also symmetric positive definite. LetĒ =LL T be any symmetric decomposition, for example, the Cholesky factorisation. Thus the reduced system reads asLL withz =L Tx . It follows that V ′ = VL −T is a full-rank matrix again. We investigate the symmetric part of the transformed system matrix in (11) The dissipativity of A implies that A + A T is negative definite by Definition 3. Due to the full rank of V ′ , the symmetric part (12) is negative definite again. The matrix V ′T AV ′ becomes dissipative and thus stable.
An important special case of Theorem 1 is given for a system (1) with an identity matrix E = I n , where explicit ODEs arise.
We assume just a non-singular matrix E now. If E is symmetric and positive definite, then we suppose a non-dissipative matrix A. Hence Theorem 1 is not applicable. The idea is to transform the system (1) into an equivalent system with a symmetric positive definite mass matrix and a dissipative matrix on the right-hand side. Let M ∈ Ê n×n be symmetric and positive definite. We multiply the system (1) by the non-singular matrix E T M and obtain the equivalent system This operation can be seen as a transformation in the range of the linear mappings. The mass matrix E T ME is always symmetric and positive definite. The matrix M has to be chosen such that the matrix E T MA is dissipative. This task is achieved by solving the generalised Lyapunov equations, cf. [19, p. 449], for any symmetric positive definite matrix F ∈ Ê n×n . Since E −1 A is a stable matrix, the Lyapunov equations (14) exhibit a unique solution M, which is also symmetric positive definite. Moreover, the matrixÂ = E T MA is dissipative in the system (13), becauseÂ +Â T = −F is negative definite. We summarise this important result in a theorem.
Theorem 2 Let the dynamical system (1) be asymptotically stable. If M is the solution of the Lyapunov equations (14) including a symmetric positive definite matrix F , then each Galerkin-type projection-based MOR of the linear dynamical system (13) yields an asymptotically stable reduced system (5).
The proof follows from Theorem 1 and the discussion above.

Numerical solution of Lyapunov equations
There are two classes of numerical linear algebra methods for the solution of the (generalised) Lyapunov equation (14): i) Direct methods: Either the solution M or its Cholesky factor L without solving for M first can be computed, see [12]. The computational effort reads as O(n 3 ). The memory requirement is about n 2 2 machine numbers.
The approximate methods yield a low-rank factor Z ∈ Ê n×q with q ≪ n On the one hand, the dimension n is high in MOR and thus the direct methods are often computationally infeasible or their effort is much higher than solving the FOM (1). On the other hand, approximate techniques typically require much less computational effort. In addition, the decomposition M ≈ ZZ T with a low-rank factor Z allows for cheap matrix-matrix multiplications with M in the projections. However, we encounter two difficulties within this approach: 1. The approximationM = ZZ T is singular by construction and thus does not describe a basis transformation.
R. Pulch (14) require a positive definite matrix F due to Theorem 1. Yet the ADI iteration is efficient only in the case of

The Lyapunov equations
We discuss the first item further. The mass matrix reads as in the reduced system (5). Since V has full rank, EV ∈ Ê n×r is also a full-rank matrix. The subspace denotes the range of the linear mapping induced by a matrix. If it holds that R(Z) ⊥ R(EV ), then we obtain Z T EV = 0 and thusĒ = 0. Although this event may happen due to q, r ≪ n, it is rather unlikely. Nevertheless, the matrix (15) can easily become ill-conditioned or singular with respect to working precision.
We observe that the stabilisation is not straightforward. An alternative technique is suggested to omit the two difficulties above in the next section.

Numerical method
A numerical approach is derived and investigated, which guarantees the preservation of stability in a Galerkin-type projection-based MOR, while Lyapunov equations can be solved approximately.

Setup of projection matrices
Let V ∈ Ê n×r with V T V = I r be a projection matrix constructed for a reduction of the system (1). We apply the Galerkin-type MOR with V to the transformed system (13). The matrices of the associated reduced system (5) can be written in the form (8) with the projection matrix Hence this reduction is of the form (8) and thus represents a special case of a (non-Galerkin type) MOR for the original system (1). However, it holds that W T V = I r in general, i.e., biorthogonality is not given.
We briefly discuss the computation work to obtain the matrices (8). The mass matrix E is often more sparse than the matrix A. Thus the matrix-matrix multiplication EV is relatively cheap. The calculation of (17) requires W = M(EV ), where M is dense. We will use an approximation of M in Section 3.2, which allows for a cheap computation of this matrix-matrix product. Given V and W , the matrices (8) can be computed as in any projection-based MOR.
It is important to note that the two Galerkin-type MORs with projection matrix V for the system (1) and the transformed system (13) are not equivalent. The reason is that the residual (10) is mapped to The property s(t) ⊥ V is not equivalent to the property E T Ms(t) ⊥ V for each t.

Alternative ansatz and Lyapunov equations
In view of Theorem 2, the construction of the projection matrices involves the efficient solution of the Lyapunov equation (14) together with the identification of an appropriate matrix F . Concerning the system (1), we assume that the matrix E −1 A is stable but non-dissipative. Definition 3 implies that the symmetric part has k ≥ 1 non-negative eigenvalues. Since the matrix E −1 A has exclusively eigenvalues with negative real part, we expect a small number k. Let µ 1 ≥ µ 2 ≥ · · · ≥ µ n be the eigenvalues of (18) and u 1 , u 2 , . . . , u n the associated eigenvectors, which form an orthonormal basis. It holds that µ max = µ 1 is the maximum eigenvalue. We define U = (u 1 , . . . , u k ) ∈ Ê n×k . The following lemma motivates a choice of the matrix F in (14).
where E is non-singular and E −1 A is a nondissipative matrix. It follows that the symmetric matrix with any real number δ > 0 is positive definite. Proof: The eigenvectors u 1 , . . . , u n of the symmetric matrix (18) represent an orthonormal basis. It follows that for j = 1, . . . , k and Hence the eigenvectors of F are u 1 , . . . , u n again. Due to µ 1 ≥ µ j for all j and µ j < 0 for j > k and µ 1 ≥ 0, all eigenvalues of F are positive.
The matrix (19) represents a rank-k update of a matrix, where the eigenvalue problem was investigated in [8], for example.
Concerning the solution of the Lyapunov equation (14), we perform the ansatz The simplification M = I n + ∆M arises in the case of E = I n . Furthermore, an approximation of (20) reads asM including a low-rank factorisation with a factor Z ∈ Ê n×q of rank q < n. Hence the approximation consists in ∆M ≈ ZZ T .
We require a positive definite matrix F . Inserting (20) and (22) in (14), we obtain the alternative Lyapunov equation Due to Lemma 1, we choose which guarantees that the matrix (22) is positive definite. The Lyapunov equations (23) can be written as with the factorŨ = √ µ max + δ U ∈ Ê n×k . An efficient approximate solution of the Lyapunov equations (25) with factorisation ∆M ≈ ZZ T and Z ∈ Ê n×q (q ≥ k) is feasible in the case of low numbers k. The determination of the matrixŨ requires the computation of the eigenvectors associated to the k dominant eigenvalues for a real symmetric matrix. Appropriate iterative methods exist for the numerical computation of just the dominant eigenvalues and their eigenvectors.
The ROM (5) includes the mass matrixĒ. We obtain a bound on its condition number.

Theorem 3
The matrixĒ = V T E TM EV including the approximation (21) and V T V = I r is symmetric positive definite. The condition number with respect to the spectral (matrix) norm · 2 exhibits the upper bound Proof: If holds that Obviously, the matrixĒ is symmetric. Let y ∈ Ê r be an eigenvector ofĒ with y 2 = 1 and η be the associated eigenvalue. It follows that η = y TĒ y = 1 + Z T EV y 2 2 = y 2 2 + Z T EV y 2 2 ≤ 1 + E 2 2 Z 2 2 due to V 2 = 1 and Z T 2 = Z 2 . Thus the eigenvalues are bounded uniformly by 1 ≤ η ≤ 1 + E 2 2 Z 2 2 . All eigenvalues are positive. Since the matrix is symmetric, the condition number satisfies which completes the proof.

Hence the upper bound (26) depends only on the magnitude of the factors E, Z
and not on the position of the subspace R(Z) in Ê n .

Symmetric decomposition
In [21], a symmetric decomposition of the solution of Lyapunov equations is used to perform a basis transformation. In Section 3.2, we introduced an approach, which does not require a symmetric decomposition of the approximation (21). Yet a numerical technique is feasible including such a symmetric decomposition with a reasonable computational effort in the case of E = I n for (1). We outline this scheme, even though it is not used for the numerical computations in Section 5.
It holds thatM x ∈ R(Z) for x ∈ R(Z) andM x = x for x ∈ R(Z) ⊥ concerning subspaces (16). Due to q ≪ n, the linear mapping described by the matrixM represents the identity on a relatively large subspace. We require a partition Ê n = R(Z) ⊕ R(Z) ⊥ . Therefore, we apply the matrix square root for the symmetric decompositionM =M 1 2M 1 2 of (21). Note that the Cholesky algorithm is inappropriate, because it calculates the complete matrix ZZ T in (21) followed by a computational effort of O(n 3 ).

R. Pulch
Theorem 4 Consider the QR-factorisation where Q ∈ Ê n×n is an orthogonal matrix, R ∈ Ê n×q and R ′ ∈ Ê q×q is an upper triangular matrix. The eigendecomposition includes an orthogonal matrix S ∈ Ê q×q and a diagonal matrix D ∈ Ê q×q with positive diagonal elements. The matrix square root of (21) exhibits the formulã with the matrices including identity matrices I j ∈ Ê j×j for j = q, n − q. Proof: Since we always assume rank(Z) = q, it follows that rank(R) = rank(R ′ ) = q.
Hence the diagonal elements of D are strictly positive. We calculatẽ We have achieved an orthogonal eigendecomposition ofM. Taking the square roots of the elements in the diagonal matrix yields the desired matrix square root (29).
The matrix square root is a symmetric matrix again. However, the basis change induced byM 1 2 is not an orthogonal transformation. We directly obtain the inverse matrix of (29) byM where just reciprocal values are required within the diagonal matrix.
Householder transformations perform the QR-factorisation of a matrix, see [27, p. 223]. Thus the matrix Q is given by successive rank-one updates. The complexity becomes 2nq 2 − 2 3 q 3 operations for (27), which can be characterised as O(nq 2 ) due to q ≪ n. Likewise, the eigendecomposition (28) is cheap, because matrices of dimension q are involved.
Basis transformations of projection matrices V require matrix-matrix multiplications, which can be performed by just a few matrix-vector productsM 4. v (4) =Sv (3) : As in step 2.
In conclusion, the effort for each matrix-vector product is negligible in comparison to the combination of QR-factorisation (27) and eigendecomposition (28).

Application to nonlinear dynamical systems
We show the applicability of the stability-preserving technique from Section 3 to nonlinear systems of ODEs.

Projection-based model order reduction
Let a nonlinear system of ODEs be given in the form with state variables x : [0, ∞) → Ê n , inputs u : [0, ∞) → Ê n in , outputs y : [0, ∞) → Ê nout , matrices B ∈ Ê n×n in and C ∈ Ê nout×n , a non-singular mass matrix E ∈ Ê n×n and a nonlinear smooth function f : Ê n → Ê n . Again initial values x(0) = x 0 are predetermined. Concerning general nonlinear dynamical system, the input-output behaviour cannot be described by a transfer function in the frequency domain.
A projection-based MOR involves matrices V, W ∈ Ê n×r of full rank again. The ROM reads asĒẋ (t) =f (x(t)) +Bu(t) y(t) =Cx(t). (31) The matrices are defined by (8). The nonlinear function is approximated bȳ Again the choice W = V yields a method of Galerkin-type.

Stability of stationary solutions
We assume an autonomous system of ODEs (30) with a constant input u(t) ≡ u 0 . Without loss of generality, we drop the term Bu in (30), because Bu 0 can be included in the function f . A stationary solution or equilibrium x * ∈ Ê n of the nonlinear dynamical system is characterised by f (x * ) = 0. We assume isolated stationary solutions. The stationary solution is asymptotically stable, if and only if the matrix is a stable matrix with respect to Definition 2, see [25, p. 22]. We assume an equilibrium x * = 0 ∈ Ê n without loss of generality. It holds that f (0) = 0.
In a projection-based MOR, the reduced system (31) features the stationary solutionx * = 0 ∈ Ê r now. The stability of this equilibrium is determined by the due to (32) and (33). Yet the stationary solution may be unstable, even if the zero-state of the FOM (30) is asymptotically stable.

Preservation of stability
We mimic the method from Section 2.3. Let M ∈ Ê n×n be symmetric and positive definite. A multiplication of the system (30) by E T M yields (assuming The state x * = 0 also represents a stationary solution of (34). We arrange the Lyapunov equation including a symmetric positive definite matrix F ∈ Ê n×n . The following theorem summarises the stability-preserving approach.
Theorem 5 Let the nonlinear dynamical system (30) with Bu = 0 have the asymptotically stable stationary solution x * = 0. Let V ∈ Ê n×r be any projection matrix with full rank. The Lyapunov equation (35) yields a solution M for some symmetric positive definite matrix F . It follows that the reduced system (31) given by (8) and (32) with the matrix W from (17) features a stationary solutionx * = 0, which is asymptotically stable.

Proof:
The projections (8) and (32) with the matrix (17) generate the reduced system The matrixĒ = V T E T MEV is symmetric and positive definite. LetĒ =LL T . Usingz =L Tx , we obtain the equivalent dynamical systeṁ The stability of the stationary solutionz * =x * = 0 is determined by the matrix . Since M is the solution of the Lyapunov equation (35), it follows that E T M ∂f ∂x x=0 is dissipative. As in the proof of Theorem 1, we conclude thatĀ ′′ is dissipative and thus a stable matrix. Now the ansatz (20) and the approximation (21) can be applied to solve the Lyapunov equation (35) approximately. The technique proceeds as in Section 3.2. In contrast to [21], a transformation is not required in the state space and thus a symmetric decomposition of the solution to the Lyapunov equation (35) can be omitted.

Illustrative examples
We apply the approach from Section 2 and Section 3 to two examples now. All numerical computations were performed by the software package MATLAB [17].

Stochastic model of mass-spring-damper system
Lohmann and Eid [16] introduced a mass-spring-damper configuration, which can be modelled by a first-order system (1) with n ′ = 8 equations. The system is single-input-single-output (SISO). In [22], this example was extended to a stochastic model, where physical parameters are replaced by random variables with uniform distributions. A polynomial chaos expansion, see [29], is used with m orthogonal basis polynomials in the random space. A stochastic Galerkin method yields a linear dynamical system (1) of dimension n = mn ′ . The system is single-input-multiple-output (SIMO). More details can be found in [22].
We consider this example again, where the uniformly distributed random parameters exhibit a variation of 10% around their mean values. All basis polynomials up to degree three are included (m = 680). In the system (1), the mass matrix E is symmetric positive definite and the matrix E −1 A is stable. However, A is not dissipative and thus Theorem 1 cannot be applied. Table 1 shows further properties of the linear dynamical system. We focus on the constant basis polynomial, which is associated to an approximation of the expected value of the original single output. Thus our stochastic Galerkin system becomes SISO, whose Bode plot is depicted in Figure 1. We want to employ a transformation to guarantee the stability in ROMs. Yet the symmetric part (18) of the matrix E −1 A exhibits k = 2720 non-negative eigenvalues, which violates the assumption of k ≪ n. Thus we do not use the approximation (21). Alternatively, we directly solve the generalised Lyapunov equation (14) by the MATLAB function lyap, where the identity matrix F = I n is chosen. Given any projection matrix V , the stability-preserving reduction is done via (8),(17).
The one-sided Arnoldi algorithm with the real expansion point s 0 = 1 yields the projection matrices V ∈ Ê n×r for each integer r now. We examine the cases r = 1, 2, . . . , 60. On the one hand, the ROM is determined by the conventional form (8) with W = V . Figure 2 (left) illustrates the spectral abscissa from Definition 1 for the matricesĒ −1Ā . It follows that just 18 of the 60 reduced systems (5) are stable. On the other hand, the stability-preserving reduction (8) with (17) is arranged. Figure 2 (right) shows the associated spectral abscissas. As expected, all reduced systems (5) are stable in this method.
We also compare the approximation quality between the conventional reduced systems and the stabilised reduced systems, because the two approaches are not equivalent even if the systems are stable. Therefore we compute the error bounds (6) for a unit norm (7) of the input ( u L 2 [0,∞) = 1), which are the H 2 -norms (3) of the differences between the transfer functions of original system and reduced system. Figure 3 depicts approximations of these upper bounds computed using a grid on the imaginary axis. We observe that the magnitudes of the error indicators are often lower and thus better in the stabilised systems.

Anemometer benchmark
The anemometer system represents a benchmark from the MOR Wiki [18]. A semi-discretisation of a convection-diffusion partial differential equation yields a linear dynamical system of the form (1) with SISO. Table 2 depicts its properties. Since the mass matrix E is a non-singular diagonal matrix, we just scale this system to obtain the case of E = I n used in the following. Now the matrix A is stable and non-dissipative. Figure 4 illustrates the Bode plot of the system. We perform an MOR by POD as described in [1, p. 277]. The input of the system (1) is chosen as the harmonic oscillation We select the initial values x(0) = 0. The MATLAB function ode45 implements a Runge-Kutta method with time step size selection based on local error control. This Runge-Kutta method performs 1602 steps and generates s = 6409 snapshots (including the internal stages) in the time interval [0, 10 −3 ]. The resulting output is depicted in Figure 5. We compute a singular value decomposition of the snap- shot matrix X ∈ Ê n×s , where just the largest singular values are determined. The r dominant singular vectors build the projection matrix V ∈ Ê n×r .
We investigate the cases r = 1, 2, . . . , 40. The conventional POD method employs the matrix W = V in (8). Figure 6 (left) shows the spectral abscissa from Definition 1 for the matricesĀ (Ē = I r ) in the ROMs (5). We observe that this MOR generates an unstable system only in the three cases r = 2, 6, 18.
Alternatively, the stability-preserving technique is used as presented in Section 3.2. We compute the largest eigenvalues of the symmetric part A + A T (E = I n ). It follows that k = 39 eigenvalues are non-negative. Hence the assumption k ≪ n is satisfied. We arrange the matrix (24) with δ = 1 and solve the Lyapunov equation (25) iteratively. We use the ADI algorithm in the function lyapchol of the sss toolbox, see [6]. Ten iteration steps yield an approximation (21) with a factor Z of rank q = 390. The projection matrices V and W from (17) produce the ROMs (5), where it holds thatĒ = I r . illustrates the spectral abscissa of the matricesĒ −1Ā . All reduced systems are asymptotically stable now.
We compare the approximation quality of the conventional reduction and the stabilisation technique. Numerical solutions of initial value problems for the FOM and the ROMs are computed by the trapezoidal rule using 1000 time steps of constant size in [0, 10 −3 ]. The maximum absolute errors are depicted for the output in Figure 7. We remind that the magnitude of the output is about 10 −5 , see Figure 5. Often the stabilised systems cause a slightly lower error in comparison to the conventional systems. Surprisingly, the instabilities do not affect the approximation quality of the ROMs in the cases r = 2, 6, 18.

Conclusions
Stability of Galerkin-type projected dynamical systems can be guaranteed by a transformation, which requires the solution of a Lyapunov equation. We derived an approximation of the solution, whose unknown part satisfies an alternative Lyapunov equation including an input matrix given by a low-rank factorisation. This approximation ensures a non-singular transformation matrix. Furthermore, the reduced mass matrices are well-conditioned in general. Numerical computations confirmed that our stability-preserving technique is feasible and does not compromise the accuracy of the projection-based MOR. Moreover, the ADI method produces an efficient iterative solution of the Lyapunov equation in the case of input matrices with factors of sufficiently small rank. However, this restriction on the rank is not always fulfilled. The rank depends on the number of non-negative eigenvalues in a symmetric matrix.