# Derivative-Extended Time Domain Reduction for Coupled Systems Using Chebyshev Expansion

Authors: Xiaolong Wang; Yaolin Jiang; Jun Liu
Journal: Discrete Dynamics in Nature and Society
 (2017)

---

## Abstract

The time domain model reduction based on general orthogonal polynomials has been presented for linear systems. In this paper, we extend this approach by taking the derivative information of the system into account in the context of model reduction of coupled systems. We expand the derivative terms over the Chebyshev polynomial basis and show that Chebyshev coefficients of the expansion possess a specific structure, making it possible to preserve much more time domain information by employing projection methods. Besides, with the well-defined projection matrices, the resulting reduced model shares the same topological structure with the original coupled system. Two numerical examples are simulated to showcase the accuracy of incorporating the derivative information into model reduction.


---

## Body

## 1. Introduction

The accurate description and increasing complexity of modeling frequently result in large scale dynamical systems, which are modeled mathematically by high order differential equations. In such large scale settings, the direct numerical simulation of these systems would be extremely time-consuming [1–3]. Model reduction reduces the large system to a moderate size system, which has much less states and retains essential properties of the original system, so as to enable a fast numerical simulation [4, 5]. Model reduction has become a powerful tool for the modeling, prediction, control, and optimization of a wide range of complex systems.In the literature, the theory and computational tools of model reduction have been studied extensively, and most of them fall into two categories: moment matching methods and balanced truncation methods. By taking Laplace transformation, moment matching methods aim to provide a Páde-type approximation to the transfer function of original systems in the frequency domain. As the approximation can be achieved implicitly by Krylov-based projection methods, this kind of methods is computationally efficient and applies to extremely large scale systems [6–9]. However, reduced models resulting from this method may be unstable. The basic idea of balanced truncation methods relies on first transforming the original model to a balancing realization which has diagonal and equal controllability and observability Gramians and then directly truncating the states corresponding to the small Hankel singular values. On the one hand, balanced truncation methods not only guarantee the stability of original models, but also provide a global error bound, which allows for an adaptive choice of the reduced order for a prescribed error tolerance. On the other hand, the computation of balancing transformation needs to solve a couple of high order Lyapunov equations, making the procedure of balanced truncation computationally intensive in large scale settings [10, 11].Recently, the direct model reduction in the time domain draws much more attentions, especially the algorithms based on orthogonal polynomials. In [12], model reduction is performed by projecting the impulse response of the original system onto a smaller dimensional subspace spanned by Chebyshev polynomials. Following the same spirit, the approach given in [13] carries out the reduction based on Laguerre polynomials and shows that a closed-form expression for the expansion coefficients can be derived more simply. Later, Laguerre function approximation is also introduced into the time domain reduction and it has been proven that in this case there exists an equivalent relationship to moment matching methods in the frequency domain [14–17]. It is not until recently that the theoretical framework of model reduction in the time domain based on general orthogonal polynomials has been established for linear systems [18]. Also, it applies to the reduction of nonlinear systems [19]. However, notice that the approach on the time domain reduction has been initially built but far more below perfection. For example, the current approach typically employs one-sided projection methods, and there still exist some degrees of freedom in the construction of reduced models. As for the structure-preserving reduction in the time domain, it has not been considered rigorously. The execution of the approach also needs to be improved. All of these issues provide the main motivation of our works.In this paper, we present a two-sided projection method for coupled systems in the time domain. To extend the existing works, we include derivatives of the time response of the system in the definition of projection matrix. Our approach towards approximating the derivatives of time responses uses Chebyshev polynomials, and we show that the Chebyshev coefficients of time responses possess a specific structure, similar to that of moments in the frequency domain but with different starting vectors. This benefits a lot the achievement of two-sided projection methods in the time domain. The topological structure of coupled systems is preserved from the perspective of subsystems by using the proposed approach. The reduced models produced by the proposed approach preserve much more time domain information, thereby providing a superior approximation. For the existing techniques on model reduction of coupled systems, we refer the reader to [20–23] and the references therein.This paper is organized as follows. Section2 explains our model and briefly reviews the existing Chebyshev approach. Section 3 examines the specific structure of Chebyshev coefficients and establishes a two-sided reduction technique in the time domain for coupled systems. The main property of the approach is proved rigorously. Numerical examples verify the theoretical analysis in Section 4. Finally, Section 5 concludes this paper.


## 2. Model Reduction of Coupled Systems

We focus on the coupled system composed ofk coupled linear time-invariable (LTI) subsystems(1)Mix˙it=Nixit+biuit,yit=ciTxit,where Mi,Ni∈Rni×ni,  bi,ci∈Rni, and  xi(t),  ui(t), and yi(t) are the state, input, and output of the subsystems, respectively, i=1,2,…,l. These subsystems are interconnected with each other by the following algebraic relationships:(2)uit=ki1y1t+⋯+kilylt+giut,yt=r1y1t+⋯+rlylt,where ri,gi,kij∈R(j=1,2,…,l) are constants. Besides, u(t) and y(t) are the input and output of the whole coupled system, respectively. For simplicity, we assume the zero initial condition xi(0)=0 for i=1,2,…,l. An example of a coupled system connected by the inputs and outputs is shown in Figure 1.Figure 1
A diagram of a system coupled through inputs/outputs.The whole coupled system (1) and (2) is an LTI system. In fact, by defining the state x(t)=x1T(t)x2T(t)⋯xlT(t)T along with the following matrices(3)Mc=diag⁡M1,M2,…,Ml,Nc=diag⁡N1,N2,…,Nl,Bc=diag⁡b1,b2,…,bl,Cc=diag⁡c1T,c2T,…,clT,R=r1r2⋯rl,G=g1g2⋯glT,K=kij,we get the following closed-form realization of the coupled system: (4)Mx˙t=Nxt+but,yt=cTxt,where coefficient matrices are expressed as(5)M=Mc,N=Nc+BcKCc,b=BcG,cT=RCc.The whole system is of order n=n1+n2+⋯+nl.For the above linear system, with the properly chosenV,W∈Rn×r, a reduced model resulting from projection methods can be defined as(6)M~x~˙t=N~x~t+b~ut,y~t=c~Tx~t,where M~=WTMV,N~=WTNV,b~=WTb, and c~=VTc. The current time domain reduction techniques typically use one-sided projection methods, that is, adopting W=V in the above definition. If the different V and W are adopted, the methods are called two-sided projection methods. In [18], the time domain reduction is presented based on the general orthogonal polynomials for linear systems. The basic idea is to approximate the impulse response of the system with orthogonal polynomials and then project the original system onto a smaller dimensional subspace spanned by the coefficient vectors of the approximation.Chebyshev polynomials are defined by the sequence [24] (7)Tit=cos⁡iarccos⁡t,-1≤t≤1,i=0,1,…, and are orthonormal regarding the weight function w(t)=1/1-t2. With the impulse function as the input function and Chebyshev polynomials for the approach proposed in [18], the impulse time response of the system is represented by the mth Chebyshev polynomial approximation(8)xt≈a0T0t+a1T1t+⋯+amTmt,where ai∈Rn are Chebyshev coefficient vectors for i=0,1,…,m and Chebyshev polynomials are rescaled as the definition interval is different than [-1,1]. By integrating the state equation of the system and plugging the approximation into the derived equality, it turns out that a0,a1,…,am satisfy the linear equation(9)Ma0=Na0-14Na1+∑i=2mNai-1i+1i2-1+b,Na0=0.5Na2+Ma1,Nai-1=Nai+1+2iMai,2≤i≤m-1,Nam-1=2mMam.Note that the recurrence relationship of Chebyshev polynomial plays a vital role in the extraction of (9). Further, the projection matrix V is defined as span⁡{V}=span⁡{a0,a1,…,am} along with W=V, leading to reduced model (6) [12, 18]. In this paper we use the notation span⁡{V} to denote the column space spanned by the columns of the matrix V, and so on.Unfortunately, the topological structure of couple systems is no longer recognizable after the above reduction procedure. More importantly, the choice ofW=V appears not to be an advisable choice to achieve superior accurate reduced models in contrast to two-sided projection methods [5, 22].


## 3. Derivative-Extended Time Domain Reduction

We study two-sided projection methods in the time domain by incorporating the derivatives of the time response into model reduction. The efficient execution of the approach benefits from the nice structure of Chebyshev coefficients.


### 3.1. Chebyshev Coefficients

Let us pay attention back to linear equation (9). Since its coefficient matrix is a block-wise Hessenberg matrix, if am is available and N is nonsingular, the rest of coefficients ai can be expressed as(10)am-1=2mN-1Mam,am-2=am+4m-1mN-1M2am,am-3=4m-1N-1Mam+8m-2m-1mN-1M3am,⋯.So instead of solving the linear equation directly, one can plug (10) into the first equation of (9) to get am first and then calculate other ai sequently with much less effort. In fact, our examination indicates that, for the time domain reduction, only the vector am is essential for our purpose.It follows from (10) by induction that, for 1≤i≤m,  am-i is a linear combination of N-1Mam,(N-1M)3am,…,(N-1M)iam when i is odd, and when i is even that is a linear combination of am,(N-1M)2am,(N-1M)4am,…,(N-1M)iam. As a result, (10) can be expressed as a compact form(11)a0⋯am-1am=N-1Mmam⋯N-1MamamΓ,where Γ is a lower triangular matrix of (m+1)×(m+1). As the diagonal elements of the matrix Γ are nonzero, equality (11) implies that there exists an invertible linear transformation between the two groups of vectors, and we conclude that (12)span⁡a0,a1,…,am=Km+1N-1M;am, where Km+1(N-1M;am) is the Krylov subspace defined by N-1M and am. It means that the time domain reduction based on Chebyshev polynomials projects the original system onto a low order Krylov subspace. This is similar to that of moment matching methods in the frequency domain, just with a different starting vector.Themth Chebyshev approximation to the system output y(t) can be obtained directly from (8). There holds y(t)≈cTa0T0(t)+cTa1T1(t)+⋯+cTamTm(t), and Chebyshev coefficients of y(t) are referred to as Yi=cTai for i=0,1,…,m. It follows from (11) that Yi satisfy(13)Y0⋯Ym-1Ym=cTN-1Mmam⋯cTN-1MamcTamΓ,which implies that preserving Yi in the reduced model amounts to the preservation of the terms cT(N-1M)iam.We extend previous works in order to consider the derivatives of the time impulse. For a given input functionu(t), assume that it has the mth Chebyshev approximation u(t)≈∑i=0muiTi(t) with scalar coefficients ui. Plugging the approximation of x(t) and u(t) into the state equation results in (14)Mx˙t≈N∑i=0maiTit+b∑i=0muiTit. If M is nonsingular, multiplying the above equation on the left by cTM-1 leads to the Chebyshev expansion (15)y˙t≈∑i=0mcTM-1Nai+cTM-1buiTit. That is, Chebyshev coefficients of the derivative are Yi(1)=cTM-1Nai+cTM-1bui for i=0,1,…,m. With the aid of (11), it has the matrix form(16)Y01⋯Ym-11Ym1=cTN-1Mm-1am⋯cTamcTN-1M-1amΓ+cTM-1bu0⋯um-1um.We proceed to the higher order derivatives. Let thejth derivative of y(t) and u(t) be y(i)(t) and u(j)(t), respectively, and let u(j)(t) have the approximation u(j)(t)≈∑i=0mui(j)Ti(t) with scalar coefficients ui(j) for i,j=0,1,…,m. Generally, there holds(17)Mxjt=Nxj-1t+buj-1t,yjt=cTM-1Nxj-1t+cTM-1buj-1t.Likewise, by induction it follows from (16) and (17) that the Chebyshev coefficients Yi(j) of the high order derivatives y(j)(t) satisfy(18)Y0j⋯Ym-1jYmj=cTM-1Nm-jam⋯cTM-1N-j+1amcTM-1N-jamΓ+cTM-1Nj-1M-1bu0⋯um-1um⋯+cTM-1bu0j-1⋯um-1j-1umj-1.Note that such an explicit expression for Chebyshev coefficients looks rather cumbersome, but coefficient matrices of the system are involved in a cyclic manner, similar to that of moments in the frequency domain [25].With the above detection on the Chebyshev coefficients, in the next subsection we introduce the derivative information into the construction of reduced models, leading to a two-sided reduction method in the time domain.


### 3.2. Main Property

Now we are in a position to present the derivative-extended time domain reduction in the framework of two-sided projection methods. For reduced model (6), let the mth Chebyshev approximation of the state and the output be (19)x~t≈∑i=0ma~iTit,y~t≈∑i=0mY~iTit. Here, we have Y~i=c~Ta~i for i=0,1,…,m. The Chebyshev approximation of the high order derivatives is denoted by(20)y~jt≈Y~0jT0t+⋯+Y~m-1jTm-1t+Y~mjTmt for j=1,2,…,m. Obviously, Chebyshev coefficients Y~i and Y~i(j) also share the same characteristic structure as shown in (13) and (18), just by replacing M,N,b,c,am with M~,N~,b~,c~,a~m.Given the inputu(t) as the impulse function, we aim to generate a reduced model such that it can preserve a desired number of Chebyshev coefficients of the impulse response and its derivatives of original systems. More precisely, we desire(21)Yi=Y~i,Yij=Y~ijfor i=0,1,…,m,j=1,2,…,m. Fortunately, it follows from (13) and (18) that the achievement of (21) can be guaranteed by the following equalities:(22)cTN-1Miam=c~TN~-1M~ia~m,cTM-1Njam=c~TM~-1N~ja~m,cTM-1Nj-1M-1b=c~TM~-1N~j-1M~-1b~for i=0,1,…,m,j=1,2,…,m.Further, one can observe that the terms in (22) possess a similar structure to moments defined in the frequency domain. For more details on moments we refer the reader to [4]. Besides, all terms in (22) tie in together. Thanks to the observation, which provides us with the insight into the definition of projection matrices V and W, we define V and W, respectively, as an orthonormal basis matrix of the following Krylov subspaces:(23)span⁡V=Km+1N-1M;am,span⁡W=Km+1M-TNT;M-Tc.Theorem 1.
LetV and W be defined as (23). If the matrices N and M are invertible, then the equalities in (22) hold; that is, the reduced model preserves the first m+1 Chebyshev coefficients of the impulse response and its jth derivatives of the original system for j=1,2,…,m, as shown in (21).Proof.
By the definition ofV, there exists a vector a^m∈Rm+1 such that am can be expressed as am=Va^m. We first prove the equality(24)N-1Miam=VN~-1M~ia^m,fori=0,1,…,m.For ease of presentation, let Xi=(N-1M)iam, and then Xi solves the linear equation (25)I-MN⋱⋱-MNX0X1⋮Xm=am0⋮0. Since Xi lies in the subspace span⁡{V}, there exists a vector X~i∈Rm+1 such that Xi=VX~i. Plugging Xi into the above linear equation and multiplying each row on the left by WT except the first row by VT lead to (26)I-M~N~⋱⋱-M~N~X~0X~1⋮X~m=a^m0⋮0.Solving the above equation iteratively, we get X~i=(N~-1M~)ia^m, which implies equality (24).
Next, we provea~m=a^m. Multiplying the first equation in (9) on the left by WT yields(27)WTMa0=WTNa0-14WTNa1+∑i=2mWTNai-1i+1i2-1+WTb.Combining (24) with (11) leads to(28)a0⋯am-1am=VN~-1M~ma^m⋯N~-1M~a^ma^mΓ.Plugging (28) into (27) leads to a linear equation on a^m, which happens to be the one for a~m corresponding to the reduced model. Owing to the unique solution to the linear equation, we get a~m=a^m.
Consequently, there exists(29)N-1Miam=VN~-1M~ia~mfori=0,1,…,m.The preservation of cT(N-1M)iam is obtained by multiplying the equation on the left by cT. Besides, with the aid of Krylov-based moment matching techniques [4], there hold(30)M-TNTjM-Tc=WWTM-TNTjM-Tc,cTM-1NM-1jW=c~TM~-1N~M~-1j.Combining (29) with (30) demonstrates that (31)c~TM~-1N~ja~m=cTM-1NM-1j-1WWTNVa~m=cTM-1Njam,c~TM~-1N~j-1M~-1b~=cTM-1NM-1j-1WWTb=cTM-1Nj-1M-1b, which concludes the proof.Remark 2.
Notice that there arem2+m conditions imposed in (21), while there are 3m+1 conditions in (22), much less than the former, which ensures the preservation of all desired Chebyshev coefficients by the two-sided projection. Moreover, the left projection matrix can be exploited for other purposes, say preserving the stability, if one just considers the preservation of Chebyshev coefficients of the impulse response.Remark 3.
Clearly, the third equality in (22) is concerned with Markov parameters, which are defined as Taylor coefficients of the transfer function in the frequency domain when it is expanded around the infinity [4]. Thus a series of Markov parameters are also guaranteed at the same time by Theorem 1, and the proposed method can be regarded as a kind of time-frequency hybrid reduction in some sense.


### 3.3. Structure-Preserving Implementation

Generally, the topological structure of the system is closely related to its practical physical structure and should be taken into account in model reduction. To this end, we partition the projection matricesV and W given in (23) according to the structure of coupled systems as follows: (32)V=V1V2⋮Vl,W=W1W2⋮Wl, where the submatrices Vi,Wi∈Rni×(m+1) for i=1,2,…,l. The resulting projection matrices are defined as(33)V^=diag⁡V1,V2,…,Vl,W^=diag⁡W1,W2,…,Wl.Starting from (5), such block-wise diagonal projection matrices immediately lead to a reduced model, composed of(34)M~ix~˙it=N~ix~it+b~iuit,y~it=c~iTx~it,where M~i=WiTMiVi,N~i=WiTNiVi,b~i=WiTbi, and c~T=cTVi. Together with the coupled relationships(35)uit=ki1y~1t+⋯+kily~lt+giut,y~t=r1y~1t+⋯+rly~lt,we obtain the whole structure-preserving reduced model.The above procedure is equivalent to projecting each subsystem onto the subspaces spanned byWi and Vi from the perspective of subsystems. So the reduction can be implemented based on each subsystem instead of reducing the entire system directly. In fact, it follows from (5) that (36)M-Tc=Mc-TCcTRT,M-TNTM-Tc=Mc-TNcTMc-TCcTRT+Mc-TCcTKTBcTMc-TCcTRT, and then we get(37)span⁡M-Tc=span⁡Mc-TCcTRT⊆span⁡Mc-TCcT,span⁡M-TNTM-Tc⊆span⁡Mc-TNcTMc-TCcT∪span⁡Mc-TCcT.Recycling (37) in the definition of W leads to (38)span⁡W=Km+1M-TNT;M-Tc⊆Km+1Mc-TNcT;Mc-TCcT. By (3), if we partition the matrix CcT=diag⁡{c1,c2,…,cl}=C1TC2T⋯ClTT with the matrix Ci∈Rni×l, Krylov vectors of the subspace Km+1(Mc-TNcT;Mc-TCcT) profit from the block-diagonal structure of Mc,Nc and have the following explicit expression: (39)Mc-TNcTjMc-TCcT=M1-TN1TjM1-TC1M2-TN2TjM2-TC2⋯Ml-TNlTjMl-TCl for j=0,1,…,m. Further, as all columns of the matrix Ci are zeros, except one column being ci, the diagonal elements Wi defined in (33) can be chosen as follows:(40)span⁡Wi=Km+1Mi-TNiT;Mi-TCi=Km+1Mi-TNiT;Mi-Tci.We point out that Wi arises only relying on the ith subsystem, while full system (23) is still required and then (33) is used to obtain Vi for the subsystems. The nice structure similar to (40) is not available for V.In the following, we sketch the procedure of the proposed derivative-extended reduction method for coupled systems.Algorithm 4 (derivative-extended reduction method for coupled systems).
Input. The inputs are coefficient matrices of (1) and (2): Mi,Ni,bi,ci;kij,ri,gi.
Output. The outputs are coefficient matrices of (34): M~i,N~i,b~i,c~i.(1)
Calculateam by solving linear equation (9).(2)
Construct the matrixVi by (23) and (33).(3)
Construct the column orthonormal matrixWi by (40).(4)
Construct subsystem (34) for i=1,2,…,l.(5)
Combine all subsystems with coupled relationships (35).Remark 5.
In Algorithm4, one can calculate the vector am by solving the first linear equation of (9) with much less effort, instead of solving the whole linear system directly. Once the vector am is found, the matrices Vi and Wi with orthogonal columns can be obtained by using the well-known Arnoldi algorithm [8].Remark 6.
Algorithm4 cannot ensure the preservation of the stability of original systems in theory. Note that even though all subsystems are stable, the whole couple system may be still unstable. The stability of reduced models can be enforced by constructing a new left projection matrix at the expense of preserving Chebyshev coefficients of the derivatives. We refer the reader to [23], where the stability of coupled systems was discussed in detail based on the closed form.In the end, we point out that while we base our discussion mainly on the single-input single-output (SISO) systems, however, all the results presented in this paper can be easily extended to the multi-input multioutput (MIMO) case. The extension is nonessential in principle and here we omit the details.


## 3.1. Chebyshev Coefficients

Let us pay attention back to linear equation (9). Since its coefficient matrix is a block-wise Hessenberg matrix, if am is available and N is nonsingular, the rest of coefficients ai can be expressed as(10)am-1=2mN-1Mam,am-2=am+4m-1mN-1M2am,am-3=4m-1N-1Mam+8m-2m-1mN-1M3am,⋯.So instead of solving the linear equation directly, one can plug (10) into the first equation of (9) to get am first and then calculate other ai sequently with much less effort. In fact, our examination indicates that, for the time domain reduction, only the vector am is essential for our purpose.It follows from (10) by induction that, for 1≤i≤m,  am-i is a linear combination of N-1Mam,(N-1M)3am,…,(N-1M)iam when i is odd, and when i is even that is a linear combination of am,(N-1M)2am,(N-1M)4am,…,(N-1M)iam. As a result, (10) can be expressed as a compact form(11)a0⋯am-1am=N-1Mmam⋯N-1MamamΓ,where Γ is a lower triangular matrix of (m+1)×(m+1). As the diagonal elements of the matrix Γ are nonzero, equality (11) implies that there exists an invertible linear transformation between the two groups of vectors, and we conclude that (12)span⁡a0,a1,…,am=Km+1N-1M;am, where Km+1(N-1M;am) is the Krylov subspace defined by N-1M and am. It means that the time domain reduction based on Chebyshev polynomials projects the original system onto a low order Krylov subspace. This is similar to that of moment matching methods in the frequency domain, just with a different starting vector.Themth Chebyshev approximation to the system output y(t) can be obtained directly from (8). There holds y(t)≈cTa0T0(t)+cTa1T1(t)+⋯+cTamTm(t), and Chebyshev coefficients of y(t) are referred to as Yi=cTai for i=0,1,…,m. It follows from (11) that Yi satisfy(13)Y0⋯Ym-1Ym=cTN-1Mmam⋯cTN-1MamcTamΓ,which implies that preserving Yi in the reduced model amounts to the preservation of the terms cT(N-1M)iam.We extend previous works in order to consider the derivatives of the time impulse. For a given input functionu(t), assume that it has the mth Chebyshev approximation u(t)≈∑i=0muiTi(t) with scalar coefficients ui. Plugging the approximation of x(t) and u(t) into the state equation results in (14)Mx˙t≈N∑i=0maiTit+b∑i=0muiTit. If M is nonsingular, multiplying the above equation on the left by cTM-1 leads to the Chebyshev expansion (15)y˙t≈∑i=0mcTM-1Nai+cTM-1buiTit. That is, Chebyshev coefficients of the derivative are Yi(1)=cTM-1Nai+cTM-1bui for i=0,1,…,m. With the aid of (11), it has the matrix form(16)Y01⋯Ym-11Ym1=cTN-1Mm-1am⋯cTamcTN-1M-1amΓ+cTM-1bu0⋯um-1um.We proceed to the higher order derivatives. Let thejth derivative of y(t) and u(t) be y(i)(t) and u(j)(t), respectively, and let u(j)(t) have the approximation u(j)(t)≈∑i=0mui(j)Ti(t) with scalar coefficients ui(j) for i,j=0,1,…,m. Generally, there holds(17)Mxjt=Nxj-1t+buj-1t,yjt=cTM-1Nxj-1t+cTM-1buj-1t.Likewise, by induction it follows from (16) and (17) that the Chebyshev coefficients Yi(j) of the high order derivatives y(j)(t) satisfy(18)Y0j⋯Ym-1jYmj=cTM-1Nm-jam⋯cTM-1N-j+1amcTM-1N-jamΓ+cTM-1Nj-1M-1bu0⋯um-1um⋯+cTM-1bu0j-1⋯um-1j-1umj-1.Note that such an explicit expression for Chebyshev coefficients looks rather cumbersome, but coefficient matrices of the system are involved in a cyclic manner, similar to that of moments in the frequency domain [25].With the above detection on the Chebyshev coefficients, in the next subsection we introduce the derivative information into the construction of reduced models, leading to a two-sided reduction method in the time domain.


## 3.2. Main Property

Now we are in a position to present the derivative-extended time domain reduction in the framework of two-sided projection methods. For reduced model (6), let the mth Chebyshev approximation of the state and the output be (19)x~t≈∑i=0ma~iTit,y~t≈∑i=0mY~iTit. Here, we have Y~i=c~Ta~i for i=0,1,…,m. The Chebyshev approximation of the high order derivatives is denoted by(20)y~jt≈Y~0jT0t+⋯+Y~m-1jTm-1t+Y~mjTmt for j=1,2,…,m. Obviously, Chebyshev coefficients Y~i and Y~i(j) also share the same characteristic structure as shown in (13) and (18), just by replacing M,N,b,c,am with M~,N~,b~,c~,a~m.Given the inputu(t) as the impulse function, we aim to generate a reduced model such that it can preserve a desired number of Chebyshev coefficients of the impulse response and its derivatives of original systems. More precisely, we desire(21)Yi=Y~i,Yij=Y~ijfor i=0,1,…,m,j=1,2,…,m. Fortunately, it follows from (13) and (18) that the achievement of (21) can be guaranteed by the following equalities:(22)cTN-1Miam=c~TN~-1M~ia~m,cTM-1Njam=c~TM~-1N~ja~m,cTM-1Nj-1M-1b=c~TM~-1N~j-1M~-1b~for i=0,1,…,m,j=1,2,…,m.Further, one can observe that the terms in (22) possess a similar structure to moments defined in the frequency domain. For more details on moments we refer the reader to [4]. Besides, all terms in (22) tie in together. Thanks to the observation, which provides us with the insight into the definition of projection matrices V and W, we define V and W, respectively, as an orthonormal basis matrix of the following Krylov subspaces:(23)span⁡V=Km+1N-1M;am,span⁡W=Km+1M-TNT;M-Tc.Theorem 1.
LetV and W be defined as (23). If the matrices N and M are invertible, then the equalities in (22) hold; that is, the reduced model preserves the first m+1 Chebyshev coefficients of the impulse response and its jth derivatives of the original system for j=1,2,…,m, as shown in (21).Proof.
By the definition ofV, there exists a vector a^m∈Rm+1 such that am can be expressed as am=Va^m. We first prove the equality(24)N-1Miam=VN~-1M~ia^m,fori=0,1,…,m.For ease of presentation, let Xi=(N-1M)iam, and then Xi solves the linear equation (25)I-MN⋱⋱-MNX0X1⋮Xm=am0⋮0. Since Xi lies in the subspace span⁡{V}, there exists a vector X~i∈Rm+1 such that Xi=VX~i. Plugging Xi into the above linear equation and multiplying each row on the left by WT except the first row by VT lead to (26)I-M~N~⋱⋱-M~N~X~0X~1⋮X~m=a^m0⋮0.Solving the above equation iteratively, we get X~i=(N~-1M~)ia^m, which implies equality (24).
Next, we provea~m=a^m. Multiplying the first equation in (9) on the left by WT yields(27)WTMa0=WTNa0-14WTNa1+∑i=2mWTNai-1i+1i2-1+WTb.Combining (24) with (11) leads to(28)a0⋯am-1am=VN~-1M~ma^m⋯N~-1M~a^ma^mΓ.Plugging (28) into (27) leads to a linear equation on a^m, which happens to be the one for a~m corresponding to the reduced model. Owing to the unique solution to the linear equation, we get a~m=a^m.
Consequently, there exists(29)N-1Miam=VN~-1M~ia~mfori=0,1,…,m.The preservation of cT(N-1M)iam is obtained by multiplying the equation on the left by cT. Besides, with the aid of Krylov-based moment matching techniques [4], there hold(30)M-TNTjM-Tc=WWTM-TNTjM-Tc,cTM-1NM-1jW=c~TM~-1N~M~-1j.Combining (29) with (30) demonstrates that (31)c~TM~-1N~ja~m=cTM-1NM-1j-1WWTNVa~m=cTM-1Njam,c~TM~-1N~j-1M~-1b~=cTM-1NM-1j-1WWTb=cTM-1Nj-1M-1b, which concludes the proof.Remark 2.
Notice that there arem2+m conditions imposed in (21), while there are 3m+1 conditions in (22), much less than the former, which ensures the preservation of all desired Chebyshev coefficients by the two-sided projection. Moreover, the left projection matrix can be exploited for other purposes, say preserving the stability, if one just considers the preservation of Chebyshev coefficients of the impulse response.Remark 3.
Clearly, the third equality in (22) is concerned with Markov parameters, which are defined as Taylor coefficients of the transfer function in the frequency domain when it is expanded around the infinity [4]. Thus a series of Markov parameters are also guaranteed at the same time by Theorem 1, and the proposed method can be regarded as a kind of time-frequency hybrid reduction in some sense.


## 3.3. Structure-Preserving Implementation

Generally, the topological structure of the system is closely related to its practical physical structure and should be taken into account in model reduction. To this end, we partition the projection matricesV and W given in (23) according to the structure of coupled systems as follows: (32)V=V1V2⋮Vl,W=W1W2⋮Wl, where the submatrices Vi,Wi∈Rni×(m+1) for i=1,2,…,l. The resulting projection matrices are defined as(33)V^=diag⁡V1,V2,…,Vl,W^=diag⁡W1,W2,…,Wl.Starting from (5), such block-wise diagonal projection matrices immediately lead to a reduced model, composed of(34)M~ix~˙it=N~ix~it+b~iuit,y~it=c~iTx~it,where M~i=WiTMiVi,N~i=WiTNiVi,b~i=WiTbi, and c~T=cTVi. Together with the coupled relationships(35)uit=ki1y~1t+⋯+kily~lt+giut,y~t=r1y~1t+⋯+rly~lt,we obtain the whole structure-preserving reduced model.The above procedure is equivalent to projecting each subsystem onto the subspaces spanned byWi and Vi from the perspective of subsystems. So the reduction can be implemented based on each subsystem instead of reducing the entire system directly. In fact, it follows from (5) that (36)M-Tc=Mc-TCcTRT,M-TNTM-Tc=Mc-TNcTMc-TCcTRT+Mc-TCcTKTBcTMc-TCcTRT, and then we get(37)span⁡M-Tc=span⁡Mc-TCcTRT⊆span⁡Mc-TCcT,span⁡M-TNTM-Tc⊆span⁡Mc-TNcTMc-TCcT∪span⁡Mc-TCcT.Recycling (37) in the definition of W leads to (38)span⁡W=Km+1M-TNT;M-Tc⊆Km+1Mc-TNcT;Mc-TCcT. By (3), if we partition the matrix CcT=diag⁡{c1,c2,…,cl}=C1TC2T⋯ClTT with the matrix Ci∈Rni×l, Krylov vectors of the subspace Km+1(Mc-TNcT;Mc-TCcT) profit from the block-diagonal structure of Mc,Nc and have the following explicit expression: (39)Mc-TNcTjMc-TCcT=M1-TN1TjM1-TC1M2-TN2TjM2-TC2⋯Ml-TNlTjMl-TCl for j=0,1,…,m. Further, as all columns of the matrix Ci are zeros, except one column being ci, the diagonal elements Wi defined in (33) can be chosen as follows:(40)span⁡Wi=Km+1Mi-TNiT;Mi-TCi=Km+1Mi-TNiT;Mi-Tci.We point out that Wi arises only relying on the ith subsystem, while full system (23) is still required and then (33) is used to obtain Vi for the subsystems. The nice structure similar to (40) is not available for V.In the following, we sketch the procedure of the proposed derivative-extended reduction method for coupled systems.Algorithm 4 (derivative-extended reduction method for coupled systems).
Input. The inputs are coefficient matrices of (1) and (2): Mi,Ni,bi,ci;kij,ri,gi.
Output. The outputs are coefficient matrices of (34): M~i,N~i,b~i,c~i.(1)
Calculateam by solving linear equation (9).(2)
Construct the matrixVi by (23) and (33).(3)
Construct the column orthonormal matrixWi by (40).(4)
Construct subsystem (34) for i=1,2,…,l.(5)
Combine all subsystems with coupled relationships (35).Remark 5.
In Algorithm4, one can calculate the vector am by solving the first linear equation of (9) with much less effort, instead of solving the whole linear system directly. Once the vector am is found, the matrices Vi and Wi with orthogonal columns can be obtained by using the well-known Arnoldi algorithm [8].Remark 6.
Algorithm4 cannot ensure the preservation of the stability of original systems in theory. Note that even though all subsystems are stable, the whole couple system may be still unstable. The stability of reduced models can be enforced by constructing a new left projection matrix at the expense of preserving Chebyshev coefficients of the derivatives. We refer the reader to [23], where the stability of coupled systems was discussed in detail based on the closed form.In the end, we point out that while we base our discussion mainly on the single-input single-output (SISO) systems, however, all the results presented in this paper can be easily extended to the multi-input multioutput (MIMO) case. The extension is nonessential in principle and here we omit the details.


## 4. Simulation Examples

In this section, two numerical examples from real world are simulated to show the performance of our approach. We compare our results with the time domain reduction method presented in [18] as well as the typical frequency domain reduction method [4]. Our simulation is carried out in Matlab Version 7.5.0.342 (R2007b) on a PC with Intel(R) Core (TM) CPU 2.60 GHz and 4 GB RAM.Example 1.
We use the benchmark model, Earth Atmospheric Example, which results from the modeling of the atmospheric storm track [26]. This model is dense and is of order 598. It is a general linear system and there is no coupled structure for this model. This example is mainly used to manifest the superior approximation accuracy when merging the derivative information into reduced models.
For comparison, the reduced model “Frequency-Sys” is produced to match the Markov parameters by the frequency domain reduction methods in [4], the reduced model “Time-Sys” by the time domain reduction method given in [18], and the last one, denoted as “Derivative-extended-Sys,” by our derivative-extended time domain method. All reduced models are of order 12. For the given input u(t)=e-2tsin⁡(5t), Figure 2 shows the output and the relative error for each reduced model in the time domain. It can be seen that the transient response of the system is well approximated by all reduced models, while our approach shows a much better accuracy, especially for the long term behavior of the system for this example.Figure 2
Outputs and relative errors for the inputut=e-2tsin⁡(5t).Example 2.
This coupled system is adapted from the standard example appearing in [22], where it describes a PI-controller used to control the temperature of a 1D beam. The coupled system is of order 2001 and comprises two subsystems.
For this example, the reduced models “Frequency-Sys” and “Time-Sys” are produced by the frequency domain reduction [4] and the time domain method [18], respectively, for comparison. “Derivative-extended-Sys” is generated by the proposed Algorithm 4 in this paper, and the topological structure of couple systems is retained naturally. Reduced order coupled systems of order 6 are first produced, and Figure 3 depicts the outputs and relative errors for the input function u(t)=e-t. Clearly, owing to preserving much more time domain information, the derivative-extended approach provides a much better approximation to the transient behavior of the system, while it shows more erratic error than the other two methods. Actually, the two-sided projection methods we adopted pave the way for retaining more information of the original systems but may result in ill-conditioned system matrices compared with one-sided methods [27].
To this end, one can proceed to define a new orthonormal matrixQi by enlarging the subspaces spanned by Wi and Vi as follows:(41)span⁡Qi=span⁡Wi∪span⁡Vifori=1,2,…,l,and then Qi can be employed to reduce each subsystem in the framework of one-sided methods along with roughly doubling the reduced order at most. With m=6, we obtain a reduced coupled system of order 15 by (41), which is denoted by “Enlarged-space-Sys.” To provide a fair comparison, three reduced models “Frequency-Sys,” “Time-Sys,” and “Derivative-extended-Sys” of order 16 are also carried out in our simulation, and the relative error of each reduced model is shown in Figure 4. We observe that the relatively higher reduced order largely improves the approximation accuracy of “Frequency-Sys” and our approach still has the best approximation to the system transient behavior, although it levels out to similar error to the others as time progresses.Figure 3
Outputs and relative errors for the inputut=e-t.Figure 4
Relative errors of reduced models with different reduced orders.


## 5. Conclusions

We have presented a derivative-extended time domain reduction method for coupled systems based on Chebyshev polynomials. The derived reduced model matches much more time domain information of the original system due to the specific structure of Chebyshev coefficients, and the coupled structure is also preserved naturally during model reduction. Theoretical analysis and numerical simulation verify the accuracy of the proposed method.



---

*Source: 8305925-2017-07-06.xml*