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

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.


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][2][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][7][8][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][15][16][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][21][22][23] and the references therein.
This paper is organized as follows.Section 2 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.

Model Reduction of Coupled Systems
We focus on the coupled system composed of  coupled linear time-invariable (LTI) subsystems where   ,   ∈ R   ×  ,   ,   ∈ R   , and   (),   (), and   () are the state, input, and output of the subsystems, respectively,  = 1, 2, . . ., .These subsystems are interconnected with each other by the following algebraic relationships: where   ,   ,   ∈ R ( = 1, 2, . . ., ) are constants.Besides, () and () are the input and output of the whole coupled system, respectively.For simplicity, we assume the zero initial condition   (0) = 0 for  = 1, 2, . . ., .An example of a coupled system connected by the inputs and outputs is shown in Figure 1.The whole coupled system (1) and ( 2) is an LTI system.In fact, by defining the state T along with the following matrices we get the following closed-form realization of the coupled system: where coefficient matrices are expressed as The whole system is of order  =  1 +  2 + ⋅ ⋅ ⋅ +   .
For the above linear system, with the properly chosen ,  ∈ R × , a reduced model resulting from projection methods can be defined as M ẋ () = Ñ x () + b () , ỹ () = cT x () , (6) where M =  T , Ñ =  T , b =  T , and c =  T .The current time domain reduction techniques typically use one-sided projection methods, that is, adopting  =  in the above definition.If the different  and  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]   () = cos ( arccos ) , − 1 ≤  ≤ 1,  = 0, 1, . . ., (7) and are orthonormal regarding the weight function 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 th Chebyshev polynomial approximation where   ∈ R  are Chebyshev coefficient vectors for  = 0, 1, . . .,  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  0 ,  1 , . . .,   satisfy the linear equation Note that the recurrence relationship of Chebyshev polynomial plays a vital role in the extraction of (9).Further, the projection matrix  is defined as span{} = span{ 0 ,  1 , . . .,   } along with  = , leading to reduced model (6) [12,18].In this paper we use the notation span{} to denote the column space spanned by the columns of the matrix , and so on.Unfortunately, the topological structure of couple systems is no longer recognizable after the above reduction procedure.More importantly, the choice of  =  appears not to be an advisable choice to achieve superior accurate reduced models in contrast to two-sided projection methods [5,22].

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.

Chebyshev Coefficients.
Let us pay attention back to linear equation (9).Since its coefficient matrix is a block-wise Hessenberg matrix, if   is available and  is nonsingular, the rest of coefficients   can be expressed as So instead of solving the linear equation directly, one can plug (10) into the first equation of ( 9) to get   first and then calculate other   sequently with much less effort.In fact, our examination indicates that, for the time domain reduction, only the vector   is essential for our purpose.It follows from (10) by induction that, for 1 ≤  ≤ ,  − is a linear combination of  −1   , ( −1 ) 3   , . . ., ( −1 )    when  is odd, and when  is even that is a linear combination of   , ( −1 ) 2   , ( −1 ) 4   , . . ., ( −1 )    .As a result, (10) can be expressed as a compact form where Γ is a lower triangular matrix of (+1)×(+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 where K +1 ( −1 ;   ) is the Krylov subspace defined by  −1  and   .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.The th Chebyshev approximation to the system output () can be obtained directly from (8).There holds () ≈  T  0  0 () +  T  1  1 () + ⋅ ⋅ ⋅ +  T     (), and Chebyshev coefficients of () are referred to as   =  T   for  = 0, 1, . . ., .It follows from (11) that   satisfy which implies that preserving   in the reduced model amounts to the preservation of the terms  T ( −1 )    .We extend previous works in order to consider the derivatives of the time impulse.For a given input function (), assume that it has the th Chebyshev approximation () ≈ ∑  =0     () with scalar coefficients   .Plugging the approximation of () and () into the state equation results in If  is nonsingular, multiplying the above equation on the left by  T  −1 leads to the Chebyshev expansion That is, Chebyshev coefficients of the derivative are  (1)    =  T  −1   +  T  −1   for  = 0, 1, . . ., .With the aid of (11), it has the matrix form [ (1)   0 ⋅ ⋅ ⋅  (1)  −1  (1)   ] We proceed to the higher order derivatives.Let the th derivative of () and () be  () () and  () (), respectively, and let  () () have the approximation () with scalar coefficients  ()  for ,  = 0, 1, . . ., .Generally, there holds Likewise, by induction it follows from ( 16) and ( 17) that the Chebyshev coefficients  ()  of the high order derivatives  () () satisfy 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.

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 th Chebyshev approximation of the state and the output be Here, we have Ỹ = cT ã for  = 0, 1, . . ., .The Chebyshev approximation of the high order derivatives is denoted by for  = 1, 2, . . ., .Obviously, Chebyshev coefficients Ỹ and Ỹ() also share the same characteristic structure as shown in ( 13) and (18), just by replacing , , , ,   with M, Ñ, b, c, ã .
Given the input () 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 for  = 0, 1, . . ., ,  = 1, 2, . . ., .Fortunately, it follows from ( 13) and ( 18) that the achievement of ( 21) can be guaranteed by the following equalities: for  = 0, 1, . . ., ,  = 1, 2, . . ., .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  and , we define  and , respectively, as an orthonormal basis matrix of the following Krylov subspaces: Theorem 1.Let  and  be defined as (23).If the matrices  and  are invertible, then the equalities in (22) hold; that is, the reduced model preserves the first  + 1 Chebyshev coefficients of the impulse response and its th derivatives of the original system for  = 1, 2, . . ., , as shown in (21).
Proof.By the definition of , there exists a vector â ∈ R +1 such that   can be expressed as   = â  .We first prove the equality For ease of presentation, let   = ( −1 )    , and then   solves the linear equation Since   lies in the subspace span{}, there exists a vector X ∈ R +1 such that   =  X .Plugging   into the above linear equation and multiplying each row on the left by  T except the first row by  T lead to Solving the above equation iteratively, we get X = ( Ñ−1 M)  â , which implies equality (24).Next, we prove ã = â .Multiplying the first equation in (9) on the left by  T yields Combining ( 24) with ( 11) leads to Plugging ( 28) into ( 27) leads to a linear equation on â , which happens to be the one for ã corresponding to the reduced model.Owing to the unique solution to the linear equation, we get ã = â .Consequently, there exists The preservation of  T ( −1 )    is obtained by multiplying the equation on the left by  T .Besides, with the aid of Krylovbased moment matching techniques [4], there hold Combining ( 29) with (30) demonstrates that which concludes the proof.
Remark 2. Notice that there are  2 +  conditions imposed in (21), while there are 3 + 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.

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 matrices  and  given in (23) according to the structure of coupled systems as follows: where the submatrices   ,   ∈ R   ×(+1) for  = 1, 2, . . ., .
The resulting projection matrices are defined as Starting from ( 5), such block-wise diagonal projection matrices immediately lead to a reduced model, composed of where M =  T      , Ñ =  T      , b =  T    , and cT =  T   .Together with the coupled relationships we obtain the whole structure-preserving reduced model.The above procedure is equivalent to projecting each subsystem onto the subspaces spanned by   and   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 and then we get Recycling (37) in the definition of  leads to By ( 3), if we partition the matrix T with the matrix   ∈ R   × , Krylov vectors of the subspace K +1 ( −T   T  ;  −T   T  ) profit from the block-diagonal structure of   ,   and have the following explicit expression: for  = 0, 1, . . ., .Further, as all columns of the matrix   are zeros, except one column being   , the diagonal elements   defined in (33) can be chosen as follows: span We point out that   arises only relying on the th subsystem, while full system ( 23) is still required and then (33) is used to obtain   for the subsystems.The nice structure similar to (40) is not available for .
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).
Remark 5.In Algorithm 4, one can calculate the vector   by solving the first linear equation of ( 9) with much less effort, instead of solving the whole linear system directly.Once the vector   is found, the matrices   and   with orthogonal columns can be obtained by using the well-known Arnoldi algorithm [8].
Remark 6. Algorithm 4 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.

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 () =  −2 sin(5), 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.
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 () =  − .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 matrix   by enlarging the subspaces spanned by   and   as follows:   and then   can be employed to reduce each subsystem in the framework of one-sided methods along with roughly doubling the reduced order at most.With  = 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.

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.

Figure 1 :
Figure 1: A diagram of a system coupled through inputs/outputs.

Figure 4 :
Figure 4: Relative errors of reduced models with different reduced orders.