Abstract
Gram matrix is an important tool in system analysis and design as it provides a description of the input-output behavior for system; its partial derivative matrix is often required in some numerical algorithms. It is essential to study computation of these matrices. Analytical methods only work in some special circumstances; for example, the system matrix is diagonal matrix or Jordan matrix. In most cases, numerical integration method is needed, but there are two problems when compute using traditional numerical integration method. One is low accuracy: as high accuracy requires extremely small integration step, it will result in large amount of computation; and another is stability and stiffness issues caused by the dependence on the property of system matrix. In order to overcome these problems, this paper proposes an efficient numerical method based on the key idea of precise integration method (PIM) for the Gram matrix and its partial derivative of linear time-invariant systems. Since matrix inverse operation is not required in this method, it can be used with high precision no matter the system is normal or singular. The specific calculation algorithm and block diagram are also given. Finally, numerical examples are given to demonstrate the correctness and validity of this method.
1. Introduction
Gram matrix contains elements which are scalar products of repeated integrals of impulse response of the system. It provides a description of the input-output behavior of the system and is therefore a potential tool in system analysis and design. Not only it determines controllability and observability of system [1], but also in several applications [2–10], it is essential. Such as in [2], Gram matrix was used in identification of linear system, in [3, 4], methods based on Gram matrix for observability restoration were presented, in [5, 6], Gram matrix was applied in state estimation and observability analysis, and in [4, 7–10] problems of state estimation using Gram matrix associated with Jacobian matrix for power system were studied. Great attention has been gained in the theoretical study and practical engineering, so it is very necessary to study its calculation.
If the system matrix is in a special form, for example, is a diagonal matrix or a Jordan matrix, or when the order of the system under consideration is not large, analytical solution of Gram matrix will be obtained [1]. But in practice, most of the system matrices are not in conformity with the forms above; the order is large or even system matrix is singular in some cases; analytical method is essentially useless; some other methods can be used to solve it. A variety of methods permit a both fast and accurate computation of the Gram matrix of system [11–16]. In [11], a method for evaluation of the Gram matrix from the coefficients of the system transfer function is proposed, but it is tedious and cumbersome; the methods proposed in [12–14] computed Gram matrix using its properties in time-domain, but as the authors stated, the technique though very elegant in principle it gives poor results due to two reasons. Firstly, it is difficult to obtain accurate measurements of impulse response date since the data obtained is usually contaminated with noise. And secondly, numerical integration techniques again introduce errors due to truncation and round off. Reference [15] presents an efficient algorithm for orthogonal decomposition of derivatives and antiderivatives of functions with rational Laplace transforms. A method for computing extended Gram matrix was proposed based on the new theorem related to Rough expansion. But it is relatively applied in a narrow scope and not suited for irrational resonant system. In [16], a method for evaluation of the Gram matrix via a Kautz model is proposed. However, the efficiency or accuracy of the computation heavily depends on the choice of Kautz parameters, and it is also only suitable for stable systems. And traditional numerical integration methods [17], such as trapezoidal integration, Simpson integration, and Gaussian integral, also can be used to compute Gram matrix. Unfortunately, they are always with low accuracy; to improve accuracy, integration step needs to be divided as small as possible; thus, computational effort will greatly increase and the solution also heavily depends on the behavior of system matrix, which will result in some troubles, such as stability and stiffness issues.
To overcome the above shortcomings and improve the accuracy of time step integration, Zhong and Williams proposed the precise integration method originally in [18]. It has been widely used in solving matrix integrals because of its high precision, stability, and fast calculation. On application to linear time-invariant dynamic homogeneous systems [19], the method can give precise numerical results approaching the exact solution at the integration points. In [20–23], it is applied to non-homogeneous dynamic systems. And Duhamel integral term is approximated using polynomial functions, sin/cos functions, exponential functions, and their combinations firstly. The corresponding recursive format is derived [20]. The accuracy of the algorithm is restricted by the inherent errors in matrix inversion and simulation of the applied loading. References [21, 22] deal with the Duhamel term using other methods, but it is still restricted by the inherent errors in matrix inversion. Reference [23] takes the idea of precise integration of exponential matrix into Duhamel integration to avoid the matrix inversion. Reference [24] reviews the precise integration method; it is not only used to solve time-invariant system, but also for time-variant, nonlinear system, and two-point boundary value problems. References [25–27] discuss the numerical stability and accuracy of precise time integration method in detail.
Unfortunately, in the above research, precise time integration method is only used to calculate the dynamic response of system or further to calculate two-point boundary value problems. It is necessary to notice that there is only one term of matrix exponent function in the above integrands [18–24]; the way how to deal with the part of matrix exponent function is the difficulty of matrix integration methods. But for Gram matrix and its partial derivative matrix, there are two terms of matrix exponent functions multiplied and separated by a matrix. Considering that general matrix operations cannot be directly combined as they are not commutative, so existing precise integration method cannot be applied to solve Gram matrix and its partial derivative matrix directly. Therefore, based on the key ideas of precise integration, this paper gives the specific computing method and derives the incremental form of the additional theorem for the computation of Gram matrix and its partial derivative matrix.
The rest of this paper is organized as follows: In Section 2, the formulation of Gram matrix and its partial derivative matrix are stated, and then the precise integration method to compute these matrices is described in detail. In Section 3, some representative numerical examples are presented and simulation results are discussed. Finally, conclusions are presented in Section 4.
2. Materials and Methods
2.1. Gram Matrix
A linear time-invariant system can be given in matrix/vector form as
in which a dot above means the differentiation with respect to time , is an -dimensional state vector, is a -dimensional input vector, and is a -dimensional output vector. , , and are , and dimensional system matrix, input matrix, and output matrix, respectively. For any time , the corresponding Gram matrix contains controllability matrix and observability matrix ; they are defined as
If and are nonsingular separately, the relevant system is completely controllable and observable. Considering that and have the similar form, therefore take as example to illustrate the computation of Gram matrix using precise integration method. To facilitate the derivation, supposing that , , , then Gram matrix corresponding to (2) is equivalent to the following forms:
In case the system corresponding to (1) is stable, then with the increasing of time , will converge to a constant value, that is to say ; the built-in function gram () in MATLAB can be used to get this constant.
2.2. Partial Derivative of Gram Matrix
If the elements of system matrix are functions of a variable, the resulting Gram matrix is also a function of this variable; then there is a partial derivative matrix to this variable, which is often required in some numerical optimization algorithm. It is also necessary to study the computation of this partial derivative matrix. Assume that the system matrix is a function of ( and are independent); that is to say
For a specific system, expressions of Gram matrix are generally complicated enough, so it is rather cumbersome to solve the partial derivatives matrix directly. Consider Therefore, make the following processing:
According to the following relationship:
(6) is equivalent to
Two different integration terms are contained in (9), and these two terms can be summarized to the computation of integration as follows:
in which the dimensions of and are equivalent. As long as (10) is solved, the partial derivative of Gram matrix corresponding to (9) will be solved. Therefore, take (10) as example to illustrate the computation of partial derivative matrix in the following.
2.3. Precise Integration Method for Computation of Gram Matrix
Now the problem is to compute (4) using precise integration method. A basic time step, denoted as , is necessary to be given for the numerical integration. The following relationship is satisfied:
According to (4), it can be obtained that between two Gram matrices, at any time and the next interval , the following relationship is satisfied:
Therefore, the problem is transformed to compute . Next, the precise integration method for will be given. According to the key point of precise integration method, the corresponding additional theorem is derived:
Equation (14) is the main key point of precise integration method. According to this equation, take as precise interval to divide the basic interval : is the weighting coefficient of matrix ; generally the computational results are accurate enough when [18, 29]; that is to say . is a small time interval, and is an extremely small time interval. Hence for the interval, the truncated Taylor expansion is applied with high precision:
Considering that is very small, the first five term series expansion should be enough. Theoretically, the solution at basic interval will be obtained through merge operations times as described by (14). But it can be seen from (17) that and both depart from the unit matrix to a very small extent. And is extremely small; if it is added to the unit matrix , it will become an appended part and its precision will be seriously dropped in the round-off operation in computer arithmetic. In order to avoid the resulting loss of accuracy, the increment and the unit matrix should be separated, and is kept in the memory rather than the matrix . So, (14) is rewritten into the following form of incremental part :
So at each step, attention will be focused on the incremental, rather than the whole amount; the precision loss caused by computer rounding operation will be avoided. And this is the second key point of precise integration method.
2.4. Precise Integration Method for Computation of Partial Derivative Matrix
The strategy of precise integration method will be used to deal with . The corresponding additional theorem is derived firstly:
In the above expression, , , and can refer to the preview results. For , finite Taylor expansion can be used to approximate The incremental form of the additional theorem is
Now, the numerical solution of is obtained using precise integration method according to (21); therefore, the exact numerical solution of partial derivatives of Gram matrix can be obtained. Since the dimensions of and are equivalent, their values do not affect the specific verification. Therefore, in the numerical examples of the next section, calculations are based on .
2.5. Programming Flow Chart of the Method
Take the computation of Gram matrix at basic interval as example, the programming flow chart is given in Figure 1.
The amount of computation mainly focuses on matrix multiplication [29] and contains two parts: one is the computation of initial basic interval, and another istimes merging during the addition theorem. The amount of computation is about , in which is the weighting coefficient of matrix (see (15)), q represents the order of Taylor expansion series (see (16), (17)), m is number of columns of matrix , and is order of matrix .
3. Numerical Examples
In this section, we verify the method presented in the preceding section with three numerical examples.
Example 1. The system matrix is , ; the integration time gradually increases with increments of 0.2 s; compute Gram matrices and .
According to
, , as these matrices are all of one-dimensional; the following results can be obtained using the analytical method:
Using the precise integration method above, taking , numerical solutions of and can be obtained. Compared with the analytical solution and calculating the errors, the results are shown in Tables 1 and 2.
From the data of tables, it can be seen that the difference between the result using method described in this paper and using analytical method is very small, so the numerical solution can be seen as the exact solution on the computer. Meanwhile, as the system is stable (poles are in the left half-plane), the steady-state value of Gram matrix can be calculated by using built-in functions gram () in MATLAB. In this example, this steady state value is 0.5. From Table 1, it can be seen that is increasing with the increase of time and gradually approaching to the steady state value. When the time t iterates to = 10 s, almost equals to the steady-state value. It also verifies the validity and effectiveness of the proposed method in this paper from another way.
Example 2. The system matrix is , ; the integration timeis taken as 0.2 s and 0.4 s, respectively, compute Gram matrixand.
According to
, , ; these matrix are of two-dimensional, and system matrix is singular that is different from the previous example. The following results can be obtained using analytical method:
With the precise integration method, taking , numerical solutions of and can be obtained. Compared with the analytical solution, the results are shown in Tables 3 and 4.
System matrix is singular and the corresponding system is also unstable (the corresponding poles are and , resp., is in the right half-plane); built-in functions gram () in MATLAB cannot be used to calculate the steady value. However, from the data in Tables 3 and 4, it can be seen that in this case methods described herein can also obtain very favorable effect, compared with the analytical method, up to ten or more significant digits. Therefore, in many cases, when other methods are unable to solve, the advantages of the method proposed in this paper are very obvious.
Example 3. In this example, the method is applied to a positioning gantry system, as shown in Figure 2. The model specification and parameters described in the following refer to paper [28]. In this model, a mass is connected to a fixed surface by a linear spring with constant . A flexible inelastic belt is connected to the mass and wraps around a pulley with radius , which is mounted on a DC motor with armature resistance and motor constant . It is assumed that the rotor inertia of the motor and the inertia of the pulley are negligible. The motor will be actuated by a voltage signal. The displacement of the mass from its original position is . The system can be modeled by the following equations [28]:
where the terms , , and are given by [28]
According to (4), using the values shown in Table 2.1 in page 29 of [28], Gram matrix can be calculated with the method proposed previews, and its results were also given by (3.21)–(3.25) in page 47-48 of [28]. Both of them are listed in Table 5.
From the data of Table 5, it can be seen that for positioning gantry system, the difference between the results is very small. This indicates that this method is effective.
4. Conclusions
This paper proposes a method using the key idea of precise integration for linear time-invariant systems to deal with the Gram matrix and its partial derivative matrix; accuracy can be compared with the exact solution in computer, the high efficiency due to making full use of the homogeneity of linear system to time, namely time-invariant properties of dynamic equation to the translation of time coordinate, and suitable for algorithm. To achieve high accuracy of results, the time is divided as small as possible. And since the method focuses on calculating the part of incremental, avoid numerical morbidity caused by large number eating decimal value. Besides, this method does not require matrix inversion, so no matter the system is singular or unstable, it can also compute correctly and effectively.
The approach is suitable for linear time-invariant system, but nonlinear and time-varying systems are much more in real life. The computation of Gram matrix and its partial derivative cannot be neglected in these cases. Methods described above bring a basis for these difficult tasks. Some kind of approximation needs to be introduced in precise integration method, such as perturbation method or making use of analytical functions to approximate. As some technique problems still remain, further research and discussions are needed on this issue.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
This project is supported by National Natural Science Foundation of China (nos. 51305321 and 51035006) and the Fundamental Research Funds for the Central Universities (nos. K5051304021 and K5051304026). These supports are gratefully acknowledged. Many thanks are due to the editor, associate editor, and reviewers for their informative, valuable, and constructive comments.