Model reduction by iterative error system approximation

ABSTRACT The analysis of a posteriori error estimates used in reduced basis methods leads to a model reduction scheme for linear time-invariant systems involving the iterative approximation of the associated error systems. The scheme can be used to improve reduced-order models (ROMs) with initial poor approximation quality at a computational cost proportional to that for computing the original ROM. We also show that the iterative approximation scheme is applicable to parametric systems and demonstrate its performance using illustrative examples.


Introduction
Model order reduction (MOR) of dynamical systems has become a central topic in the computational sciences and engineering. In particular for linear time-invariant (LTI) systems, many MOR methodologies have been proposed in recent years. For a survey we refer to the books [1][2][3] and the recent survey paper [4]. Here, we will discuss the situation that a reduced-order LTI model has been computed by any projection-based method, and it is desirable to improve its approximation quality. This may be due to the fact that design specifications are not met with the original reduced-order model (ROM), or an a posteriori error analysis shows insufficient accuracy of the approximation. The goal is to employ the already obtained information from the ROM and an a posteriori error analysis in order to improve the approximation quality of the ROM in a systematic manner. Compared with the reduced basis (RB) method, which updates the ROM iteratively, the approach proposed here is more general in that any appropriate model reduction technique can be applied in the iteration steps, and the way the ROM is improved differs from the RB method in several aspects. We are only aware of the approach in [5,6], which yields a procedure similar to the one we suggest in this paper, but is derived from a different perspective. (The relations and differences of both approaches are discussed later in this paper.) The present paper has been inspired by the RB approach to MOR of parametrized systems [7]. In particular a careful analysis of the methodology presented in [8] leads to the developments reported later. The essence of our approach is an iterative method which at each step is performing MOR to the error system; as a result this error is reduced at each step of the iteration in a targeted fashion, thus yielding good reduced-order systems. Due to the different derivation of the method as compared to [5,6], it also becomes clear quickly that an extension to linear parametric problems is straightforward.
The paper is structured as follows. An analysis of the a posteriori error bounds in [8] leads to Lemma 2.1 which reveals the new structure of the error system. Consequences follow in Section 2.1 and the proposed new procedure is described in Section 3. The paper concludes with numerical examples illustrating the performance of the method and a generalization to parametric systems.

Some preliminaries
Consider the single-input single-output (SISO) full-order model (FOM) of order n (i.e. the state xðtÞ is in R n at any time t, A; E; 2 R nÂn , B; C T 2 R n ) with transfer function HðsÞ ¼ CðsE À AÞ À1 B. We construct the (SISO) ROM E d dtx ¼Âx þBu;ŷ ¼Ĉx (2) of order k by means of a Petrov-Galerkin projection defined by V; W 2 C nÂk , i.e.
The transfer function of the ROM, thus, isĤðsÞ ¼ĈðsÊ ÀÂÞ À1B . Assume for the time being that E ¼ I, and W Ã V ¼ I (as in [8]). Given the FOM and ROM defined by Equations (1) and (2), respectively, the error and residual in time domain are defined as follows: Hence, the temporal evolution of the error is described by Throughout this paper, we assume that A and all reduced order matricesÂ (or the associated matrix pairs ðA; EÞ, ðÂ;ÊÞ) have all their eigenvalues in the open left half plane. This implies that the associated linear systems are (asymptotically) stable, as well as that the constant γ 1 in the following proposition is finite.
Remark 1.1 Denoting the Laplace transforms ofx; u byX; U, the residual in frequency domain is where ΦðsÞ ¼ sI À A. This expression (with UðsÞ ¼ 1) will appear prominently in the sequel.
1.2. Improved error estimation using the dual system Following [8], we define the dual system on the interval ½0; T, with fixed final time T. The dual state will be denoted by x du : Given the dual projection defined by V du , In this framework, the projection defined by V, W in Section 1.1 is sometimes referred to as the primal projection. The resulting error and residual of the dual FOM and ROM are: e du ðtÞ ¼ x du ðtÞ À W duxdu and (9) Thus, the dual error satisfies the differential equation: d dt e du ðtÞ ¼ A du e du ðtÞ þ r du ðtÞ; where e du ðTÞ ¼ ðI À W du V du Ã ÞC Ã : Proposition 1.2 (Dual a posteriori error estimate [8].) Let γ du 1 ! sup t2½0;T k expðÀA du tÞ k . Then the following bound holds: The next result provides an exact expression for the output error y Àŷ, where the two outputs are defined by Equations (1) and (2), respectively. This expression involves both the primal and the dual reduced systems. In the sequel, the inner product of the complex vectors x, y is denoted by hx; yi ¼ x Ã y. Lemma 1.1 (Output error equality [8].) The reduced output error at time T satisfies At this point (following [8]), the primal and dual systems are combined to define an improved reduced system output at time T, namely: Our goal in the sequel is to analyse this expression so as to obtain an explicit reduced-order system and hence an explicit error system for this improved reduced system output.

Analysis of the improved reduced system
In order to determine the linear system behind Equation (13), we start with some consequences of the above definitions: CVe W Ã AVðTÀτÞ W Ã BuðT À τÞdτ; In order to investigate the term T 4 , we note that W duxdu ðtÞ ¼ W du e ÀÂ du ðTÀtÞxdu ðTÞ ¼ W du e ðV du Þ Ã A Ã W du ðTÀtÞ ðV du Þ Ã C Ã Therefore, with rðtÞ ¼ ðAV À VÂÞxðtÞ þ ðI À VW Ã ÞBuðtÞ, we get CV du e ðW du Þ Ã AV du ðTÀtÞ ðW du Þ Ã ðAV À VÂÞxðtÞdt: Thus, T 2 þ T 4 represents the analytical solution to the system: Hence putting the equations together yields the following augmented reduced system, where the reduced state of the primal systemx is redefined asx 1 : where Hence, Rewriting these equations with ¼x 1 Thus defining and re-inserting EÞI, we get the following result.
Lemma 2. 1 The system with outputỹ as defined in (13) has the following generalized state form: From the above lemma, it follows for the error of the transfer functionH associated to the improved reduced system with outputỹ: HðsÞ ÀHðsÞ ¼ HðsÞ À ½CV 1 |{z} Here, ΦðsÞ ¼ sE À A.
The transfer function of the error system is Time domain representation of the error system x a ðtÞ ¼ A a x a ðtÞ þ B a uðtÞ; yðtÞ ¼ C a x a ðtÞ ) yðsÞ ¼ HðsÞuðsÞ; in other words the augmented system is a non-minimal realization of the original system ðE; A; B; CÞ. Consider the augmented error system: E a d dt e a ðtÞ ¼ A a e a ðtÞ þ r a ðtÞ; y err ðtÞ ¼ C a e a ðtÞ; where r a ðtÞ is the associated augmented residual defined by

Some comments and consequences
(a) The interpretation of the ROM obtained here can be done without reference to the dual. This system namely has triangular structure (the first projection affects the second, but not vice-versa).
In addition, while the improved system output in Equation (13) is defined only for time T, here, this restriction is lifted andỹðtÞ in Lemma 2.1 is defined for all time.
(b) The error system (Equation (14)) factors in a product of two residues, the first coming from the original projection and the second from the second (also referred to as dual) projection. Specifically, the transfer function H err ðsÞ can be rewritten as As we have assumed SISO systems, it is straightforward to get where R pr ðsÞ is the frequency-domain expression of the residual rðtÞ of the ROM (2), and R du ðsÞ is the frequency-domain expression of the residual r du ðtÞ of the reduced dual system (Equation (8)). Essentially, the error bound ΔðsÞ decays quadratically as it is the product of the two residuals jjR du ðsÞjj 2 and jjR pr ðsÞjj 2 . H err ðsÞ j j is the true error, often decaying faster than quadratic, and sometimes, it is much smaller than the product of jjR du ðsÞjj 2 and jjR pr ðsÞjj 2 , see Figure 7 for an experimental illustration, where H ð2Þ err is the error system (Equation (17)) of a parametric model. (c) The primal and dual projections can be of different dimension k, ,, respectively; the dimension of the ROM is then k þ ,. Notice also that the ROM above cannot be obtained by means of an (explicit) Petrov-Galerkin projection applied to (E,A,B,C). Instead the (block diagonal) projection defined by V , W , must be applied to a non-minimal realization of the (d) The error system (14) can be written as Thus, one suggestion for making use of the above formula is Step 1. Choose the first (primal) projection V 1 , W 1 .
Step 2. Compute the second (dual) projection V 2 , W 2 by means of weighted reduction, where the to-be-reduced system is CΦ À1 ðsÞ and the (right) weight is H ð1Þ w . (e) The output of the error system (Equation (15)) describes the time-domain output error of the improved ROM defined in Lemma 2.1. In Section 4, we plot the time-domain output errors of the successively constructed ROMs for a parametric system in Figure 6, where the monotonic decay of the errors can be observed.
(f) In this framework, it readily follows that three or more stages can be considered. In the case of three stages, we project the system The transfer function of the associated error system is Thus, the to-be-approximated system here is still CΦ À1 ðsÞ, while the weighting function is H ð21Þ w . Hence, we may conclude that this method leads to model reduction by successive approximation of the ensuing error systems.
(g) Some remarks on the related literature are in order.
(1) Some of the above considerations have been used in [9], e.g. some of the discussions in the comments (b) and (c). (2) In [5,6], Panzer, Wolf, and Lohmann have obtained an expression similar to Equation (14) and its generalization to more than two stages. The motivation that led to these results, however, is different from the reduced-basis motivation used here. The papers above contain an additional result which is not used here: namely the fact is used that if V is constructed by means of a rational Krylov method based on E, A, B, then ΠAV is rank one, and in particular it can be factored as ΠAV ¼ ðΠBÞðCÞ, whereC is a row vector. In these expressions Π ¼ I À EVðW Ã EVÞ À1 W Ã , and W is arbitrary.
3. The new procedure: weighted reduction of successive error systems In the following, we give a detailed description of the iterative procedure to improve a ROM using successive approximation of the error system. We use the following notation throughout: (A) Data: given is the system (C,E,A,B) of McMillan degree n, with transfer function (B) 1st step: find V 1 ; W 1 2 R nÂk 1 , and construct the corresponding ROM: It follows that the error system can be written as Introducing the projection we can eliminate ΦðsÞ as follows: Hence the error system above can be expressed as follows, where H ð1Þ w ðsÞ is the weight: (C) 2nd step: we now seek V 2 ; W 2 2 R nÂk 2 , where the resulting reduced system with It readily follows that the resulting error system is Consequently, the second step amounts to the approximation of CΦ À1 ðsÞ, subject to the weight (Note that CΦ À1 H ð1Þ w is nothing but the first error system H ð1Þ err ðsÞ (Equation 19). Therefore, approximation of CΦ À1 ðsÞ, subject to the weight H ð1Þ w ðsÞ, implicates that V 2 ; W 2 are obtained from applying MOR to the error system H ð1Þ err ðsÞ.) (D) q-th step: we seek V q ; W q 2 R nÂk q , where the resulting reduced system with It readily follows that the resulting error system is H ðqÞ err ðsÞ ¼ HðsÞ À H ðqÞ red ðsÞ it is easy to see that the matrices V i , W i are obtained by performing MOR to the error system H ðiÀ1Þ err ðsÞ, whereĤ ðiÀ1Þ err ðsÞ is the resulting reduced error system. After V i , W i are computed, they are combined with V j ; W j ; j < i, to form V 1i ; W 1i , so as to construct the ith updated ROM H ðiÞ red ðsÞ of the original system.
For i > 1, if the balanced truncation (BT) method is used to compute V i ; W i , i ¼ 1; . . . ; q, it consists of applying BT to the error system H ðiÀ1Þ err ðsÞ (23), i.e. the weighted system C; E; A; B ðiÞ w . If an interpolatory method is used, such as moment-matching (MM), selection of good interpolation points is important. At each step i > 1 of updating the ROM, the interpolation point can be selected as the one which corresponds to the largest magnitude of H ðiÀ1Þ err ðsÞ. The idea of selecting the interpolation point is analogous to the greedy algorithm of the RB method [7], where at each iteration step, the parameter corresponding to the biggest error is selected to compute a new basis vector, and usually a single vector is added to enrich the RB. For MM, several vectors could be selected at once.  The new procedure explored in the current section can, more or less straightforwardly, be applied to these parametric systems by replacing ΦðsÞ with ΦðμÞ, C with CðμÞ, and B with BðμÞ. For the second-order system in Equation (25), in addition to the block matrices of A introduced in Section 2.1, the block matrices of M, K should also be introduced to compute the updated ROMs as explained in the comment (c) in Section 2.1. The block structures of M.K,A are the same as the block structures of E,A for the first-order system. In the next section, we also demonstrate the new procedure applied to an example of a second-order parametric system.

Examples
We will now illustrate the above considerations by means of three examples. For the weighted model reduction needed in the successive approximation procedure, we can use weighted balanced truncation (W-BT). For completeness, we sketch this method here, more details can be found in [1].
Let (C,E,A,B,D) be a realization of the original system and (C w , E w , A w , B w ,D w ) a realization of the weighting function. 1 We define the composite system with transfer func-tionĤðsÞ ¼ HðsÞH w ðsÞ: Let P be the controllability Gramian of ðÊ;Â ;BÞ, and Q the observability Gramian of (C,A). Finally, let P n be the leading n Â n minor of P. With P n ¼ U Ã U and Q ¼ L Ã L, we compute the SVD UL Ã ¼ WΣV Ã . The projection matrices are then Example 1. The first system has order 16 and is defined as follows: A ¼ blkdiag À:1 40 À40 À:1 ! ; À:01 25 À25 À:01

!
; À:02 10 À10 À:02 ! ; B ¼ onesð16; 1Þ; C ¼ ½2; 1; À1; 3; 1; À1; À1; À2; À2; 5; 3; 1; À1; À2; À4; 1; The W-BT procedure described above is applied in three steps. The reduction order at each step is 2, i.e. the three reduced systems have order 2, 4, 6. In Table 1, the H 2 -and the H 1 -norms of the three error systems are shown, together with the norms of the sixth order BT (last row of the table). Here, Σ is the original system, Σ 2 , Σ 4 , Σ 6 are the systems obtained by successive W-BT, and Σ bal 6 is the sixth order reduced system obtained by BT (listed for comparison purposes). Recall that BT of order 6 involves a full, as opposed to block triangular, projection and hence the reduced system has a better fit than the one obtained by a three step-reduction where each step is of order 2.
4.2125e+00 7.4014e+00 k Σ À Σ 6 k 9.5388e-01 1.3787e+00 k Σ À Σ bal 6 k 9.5371e-01 1.3790e+00 Example 2. The second system is the well-known CD player model with McMillan degree 120 [10]. The reduction is also performed in three steps (orders of the resulting ROMs: 6; 12; 18). We apply two reduction methods for each successive step. The first is W-BT as described earlier and the second is Rational Krylov, where the interpolation points are chosen from the peaks of each error system. The necessity for the second method arose from the fact that W-BT missed the second peak of the amplitude Bode plot (between 10 4 and 10 5 Hz); this is due to the low magnitude of this part of the frequency response compared with the magnitude of the dominant peak. Capturing the secondary peak is thus achieved at the expense of obtaining higher error norms as shown in Table 2. There, Σ is the original system; Σ 6 , Σ 12 , Σ 18 are the systems obtained by successive W-BT or rational Krylov; and Σ bal 18 is the 18 th -order reduced system obtained by BT. Again, the monotonicity of the decay of these norms holds. Finally, for comparison purposes, we add the norm of the 18th-order reduced system obtained by BT.
The corresponding Bode plots of the FOM and ROM transfer functions as well as the resulting error systems are shown in Figures 2 and 3.
Example 3. We study a parametric model of a plate (floor inside a building near a highway). The parametric model is of second-order form, The parameter p is the damping coefficient. The size of the original system is n ¼ 22.
We use the parametric MOR algorithm proposed in [11] to compute the ROM. From the algorithm, a projection matrix V is computed, and the ROM is obtained by Galerkin projection as below, We have computed H ð1Þ err ðs; pÞ, H ð2Þ err ðs; pÞ, and H ð3Þ err ðs; pÞ, which are the transfer functions of the error systems of the first, second and the third ROMs of sizes r ¼ 4; 7; 10, respectively. The ith ROM is obtained by interpolating the i À 1st error system H ðiÀ1Þ err , where the interpolation point of ðs; pÞ corresponds to the peak of H ðiÀ1Þ err , for i ¼ 2; 3. In Figure 4 on the left, we compare H ð1Þ err ðs; pÞ with H ð2Þ err ðs; pÞ, and in Figure  From Figure 4(a), we see that the peak of H ð1Þ err is close to f ¼ 140Hz. After interpolation of H ð1Þ err at the point close to f ¼ 140Hz, the error of the second ROM, i.e. the error system H ð2Þ err , has the lowest magnitude at that interpolation point. Similarly in the right plot of Figure 4, the peak of H ð2Þ err is at a point close to f ¼ 40Hz. After interpolating H ð2Þ err at that point, the error of the third  ROM, i.e. the error system H ð3Þ err , has the lowest magnitude at the same point. It is entirely in agreement with the theoretical analysis in Remark 3.1. The maximal H 1 -norm of the error systems over the parameter sample space p is given in Table 3.
In Figure 5, we plot the magnitudes of the transfer functions of the original system Σ, and those of the first ROM Σ r¼4 , and the second ROM Σ r¼7 at 500 random samples of frequency and the parameter p. The second ROM is already indistinguishable from the original model.
In addition, we plot the error behaviour of the ROMs in the time domain in Figure 6. The plots clearly describe the behaviour of the output error according to the time evolution as well as the parameter variation.
In Figure 7, the norms of the two residuals jjR du jj 2 and jjR pr jj 2 , as well as the magnitude of H ð2Þ err are plotted, where H ð2Þ err is indeed smaller than the product of jjR du jj 2 and jjR pr jj 2 , as implied by the error bound discussed in the comment (b) in Section 2.1.

Conclusions
We have discussed an iterative procedure to improve a ROM of a LTI or parametric system computed by an arbitrary method. The procedure is based on successively reducing the error systems associated with the current reduced systems. The formulas used are derived by analysing known a posteriori error estimates used in RB methods. The iterative process uses either weighted BT or rational Krylov methods to reduce the error systems, represented as weighted transfer functions. Numerical examples illustrate the performance of the suggested procedure. Note 1. If the weight is the product of two transfer functions defined by means of (C i ;E i ;A i ;B i ;D i ); i = 1; 2, then

Disclosure statement
No potential conflict of interest was reported by the authors.

Funding
This work was supported by NSF Grants CCF-1017401 and CCF-1320866, as well as the DFG Grants AN 693/1-1 and BE 2174/18-1.