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.

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 [13]. 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 [69]. 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 [1417]. 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 [2023] 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.

2. Model Reduction of Coupled Systems

We focus on the coupled system composed of coupled linear time-invariable (LTI) subsystemswhere ,  , and  ,  , and are the state, input, and output of the subsystems, respectively, . These subsystems are interconnected with each other by the following algebraic relationships:where are constants. Besides, and are the input and output of the whole coupled system, respectively. For simplicity, we assume the zero initial condition for . 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 along with the following matriceswe get the following closed-form realization of the coupled system: where coefficient matrices are expressed asThe whole system is of order .

For the above linear system, with the properly chosen , a reduced model resulting from projection methods can be defined aswhere , and 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] 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 approximationwhere are Chebyshev coefficient vectors for and Chebyshev polynomials are rescaled as the definition interval is different than . By integrating the state equation of the system and plugging the approximation into the derived equality, it turns out that satisfy the linear equationNote that the recurrence relationship of Chebyshev polynomial plays a vital role in the extraction of (9). Further, the projection matrix is defined as along with , leading to reduced model (6) [12, 18]. In this paper we use the notation 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].

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 is available and is nonsingular, the rest of coefficients can be expressed asSo 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 ,   is a linear combination of when is odd, and when is even that is a linear combination of . As a result, (10) can be expressed as a compact formwhere is a lower triangular matrix of . 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 is the Krylov subspace defined by 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 , and Chebyshev coefficients of are referred to as for . It follows from (11) that satisfywhich implies that preserving in the reduced model amounts to the preservation of the terms .

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 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 leads to the Chebyshev expansion That is, Chebyshev coefficients of the derivative are for . With the aid of (11), it has the matrix form

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 . Generally, there holdsLikewise, by induction it follows from (16) and (17) that the Chebyshev coefficients of the high order derivatives satisfyNote 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 th Chebyshev approximation of the state and the output be Here, we have for . The Chebyshev approximation of the high order derivatives is denoted by for . Obviously, Chebyshev coefficients and also share the same characteristic structure as shown in (13) and (18), just by replacing with .

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 desirefor Fortunately, it follows from (13) and (18) that the achievement of (21) can be guaranteed by the following equalities:for

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 Chebyshev coefficients of the impulse response and its th derivatives of the original system for , as shown in (21).

Proof. By the definition of , there exists a vector such that can be expressed as . We first prove the equalityFor ease of presentation, let , and then solves the linear equation Since lies in the subspace , there exists a vector such that . Plugging into the above linear equation and multiplying each row on the left by except the first row by lead to Solving the above equation iteratively, we get , which implies equality (24).
Next, we prove . Multiplying the first equation in (9) on the left by yieldsCombining (24) with (11) leads toPlugging (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 existsThe preservation of is obtained by multiplying the equation on the left by . Besides, with the aid of Krylov-based moment matching techniques [4], there holdCombining (29) with (30) demonstrates that which concludes the proof.

Remark 2. Notice that there are conditions imposed in (21), while there are 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 matrices and given in (23) according to the structure of coupled systems as follows: where the submatrices for . The resulting projection matrices are defined asStarting from (5), such block-wise diagonal projection matrices immediately lead to a reduced model, composed ofwhere , and . Together with the coupled relationshipswe 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 getRecycling (37) in the definition of leads to By (3), if we partition the matrix with the matrix , Krylov vectors of the subspace profit from the block-diagonal structure of and have the following explicit expression: for Further, as all columns of the matrix are zeros, except one column being , the diagonal elements defined in (33) can be chosen as follows: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).
Input. The inputs are coefficient matrices of (1) and (2):
Output. The outputs are coefficient matrices of (34): (1)Calculate by solving linear equation (9).(2)Construct the matrix by (23) and (33).(3)Construct the column orthonormal matrix by (40).(4)Construct subsystem (34) for .(5)Combine all subsystems with coupled relationships (35).

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.

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 , 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 , 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.

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.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

This work was supported by the Natural Science Foundation of China (NSFC) under Grants 11401472 and 11371287, the Natural Science Foundation of Shaanxi Province (2017JM1019), and the Fundamental Research Funds for the Central Universities under Grant 3102015ZY072.