Combined Preorder and Postorder Traversal Algorithm for the Analysis of Singular Systems by Haar Wavelets

Abstract

An efficient computational method is presented for state space analysis of singular systems via Haar wavelets. Singular systems are those in which dynamics are governed by a combination of algebraic and differential equations. The corresponding differential-algebraic matrix equation is converted to a generalized Sylvester matrix equation by using Haar wavelet basis. First, an explicit expression for the inverse of the Haar matrix is presented. Then, using it, we propose a combined preorder and postorder traversal algorithm to solve the generalized Sylvester matrix equation. Finally, the efficiency of the proposed method is discussed by a numerical example.

1. Introduction

Wavelets are mathematical functions that cut up data into different frequency components and then study each component with a resolution matched to its scale. Wavelets are now being applied in many areas of science and engineering [1–4]. Much attention has been focused on the use of wavelet transforms to investigate dynamic systems. This is due to the powerful ability of wavelet transforms to decompose time series in time-frequency domain and wavelet basis functions. Chen and Hsiao [3, 4] derived a Haar operational matrix for integration and solved the lumped and distributed parameter systems by constructing operational matrices of various order. The main characteristic of this technique is that it converts a differential equation into an algebraic one with the result that the solution procedures are greatly reduced and simplified. This approach gives new insight into the use of the Haar wavelet method.

Singular systems (also referred to as descriptor or semistate systems) arise more naturally than do state-variable descriptions in the analysis of many sorts of systems. Examples occur in electrical networks, neural networks, control systems, chemical systems, economic systems, and so on (see [5, 6] and references therein). These systems are governed by a mixture of differential and algebraic equations. The complex nature of singular systems causes many difficulties in the analytical and numerical treatment of such systems.

Recently, Haar wavelet technique was applied to state analysis and observer design of singular systems [7]. This approach replaces the state function and the forcing function by the truncated Haar series, respectively. Then the state trajectories are obtained by solving a generalized Sylvester matrix equation. But there exists a trade-off between the resolution of the wavelets and the computation time. The accuracy of the solution can be achieved by increasing the resolution level, but this requires more computation time and very large memory.

In this paper, an efficient computational method is presented for state space analysis of singular systems via Haar wavelets. First, an explicit expression for the inverse of the Haar matrix is presented. This inverse matrix also has a recursive structure. By using this matrix, we propose a combined preorder and postorder traversal algorithm. Then, the full-order generalized Sylvester matrix equation should be solved in terms of the solutions of simple linear matrix equations. Finally, the efficiency of the proposed method is discussed by a numerical example.

2. Kronecker Product

Let 𝐀=[𝑎_𝑖𝑗] and 𝐁=[𝑏_𝑖𝑗] be 𝑛×𝑝 and 𝑟×𝑞 matrices, respectively. The Kronecker product of the matrices, denoted by 𝐀⊗𝐁(∈𝑅^{𝑛𝑟×𝑝𝑞}), is defined as 𝑎𝐀⊗𝐁=¹¹𝐁𝑎₁₂𝐁⋯𝑎_1𝑝𝐁𝑎₂₁𝐁𝑎₂₂𝐁𝑎_2𝑝𝐁𝑎⋮⋮_𝑛1𝐁𝑎_𝑛2𝐁⋯𝑎_𝑛𝑝𝐁. (2.1) The 𝑣𝑒𝑐 operator transforms a matrix 𝐀 of size 𝑛×𝑝 to a vector of size 𝑛𝑝×1 by stacking the columns of 𝐀. Some properties of the Kronecker product are given below [8]: (𝐀+𝐁)⊗𝐂=𝐀⊗𝐂+𝐁⊗𝐂,(𝐀⊗𝐁)𝐂=(𝐀𝐂⊗𝐁),(𝐀⊗𝐁)(𝐂⊗𝐃)=(𝐀𝐂⊗𝐁𝐃),(𝐀⊗𝐁)^𝑇=𝐀^𝑇⊗𝐁^𝑇. (2.2)

3. Haar Wavelets and Their Properties

Wavelets constitute a family of functions constructed from dilation and translation of a single function called the mother wavelet that generates orthogonal bases of 𝐿₂(𝑅). The simplest and most basic of the wavelet systems is the Haar wavelet which is a group of square waves with magnitudes of ±1 in certain intervals and zeros elsewhere [9]. The scaling function 𝜑₀(𝑡) and mother wavelet 𝜑₁(𝑡) are defined by, respectively, 𝜑₀𝜑(𝑡)=1,𝑡∈[0,1),0,𝑡∉[0,1),₁1(𝑡)=1,𝑡∈[0,21),−1,𝑡∈[2,1),0,𝑡∉[0,1). (3.1) Then, all the other basis functions 𝜑_𝑘(𝑡) are obtained by dilation and translation of the mother wavelet as follows: 𝜑_𝑘(𝑡)=𝜑₁(2^𝑛𝑡−𝑗)=1,𝑡∈[𝑡_𝑎,𝑡_𝑏),−1,𝑡∈[𝑡_𝑏,𝑡_𝑐),0,𝑡∉[𝑡_𝑎,𝑡_𝑐), (3.2) where 𝑘=2^𝑛+𝑗, integer 𝑛≥1 is a dilation parameter, integer 0≤𝑗<2^𝑛 is a shift parameter, and the intervals are given by 𝑡_𝑎=𝑚/2^𝑛, 𝑡_𝑏=(0.5+𝑗)/2^𝑛, and 𝑡_𝑐=(1+𝑗)/2^𝑛. Since the support of the Haar wavelet is [0,1), any square integrable function 𝑦(𝑡)∈𝐿₂[0,1) can be written as an infinite linear combination of Haar functions 𝑦(𝑡)=^∞_𝑘=0𝑐_𝑘𝜑_𝑘(𝑡),𝑡∈[0,1), (3.3) where the Haar coefficients are determined by 𝑐_𝑘=⟨𝑦(𝑡),𝜑_𝑘(𝑡)⟩=2^𝑛1₀𝑦(𝑡)𝜑_𝑘(𝑡)𝑑𝑡, (3.4) where ⟨⋅,⋅⟩ denotes the inner product. In practical applications, Haar series are truncated to 𝑚 terms, that is, 𝑦(𝑡)≅^𝑚−1_𝑘=0𝑐_𝑘𝜑_𝑘(𝑡)=𝐂^T_𝑚𝐡_𝑚(𝑡), (3.5) where Haar functions coefficient vector 𝐂_𝑚 and Haar functions vector 𝐡_𝑚 are defined as 𝐂_𝑚≜[𝑐₀𝑐₁⋯𝑐_𝑚−1]^T and 𝐡_𝑚(𝑡)≜[𝜑₀(𝑡)𝜑₁(𝑡)⋯𝜑_𝑚−1(𝑡)]^T.

Integrals of the Haar functions with respect to variable 𝑡 form ramp and triangular waveforms standing with uniform slope, respectively, at the positions of the corresponding rectangular functions. The group of these integrals can be expressed as follows: ¹₀𝜑₀(𝑡)𝑑𝑡=𝑡,𝑡∈[𝑡_𝑎,𝑡_𝑏),₁₀𝜑_𝑘(𝑡)𝑑𝑡=𝑡−𝑡_𝑎,𝑡∈[𝑡_𝑎,𝑡_𝑏),−𝑡+𝑡_𝑐,𝑡∈[𝑡_𝑏,𝑡_𝑐),0,𝑡∉[𝑡_𝑎,𝑡_𝑐). (3.6) Then, the Haar matrix 𝐇_𝑚 is defined as 𝐇_𝑚(𝑡)≜𝐡_𝑚(𝑡₀)𝐡_𝑚(𝑡₁)⋯𝐡_𝑚(𝑡_𝑚−1), (3.7) where 𝑖/𝑚≤𝑡_𝑖≤(𝑖+1)/𝑚.

Integration of the Haar function vector can be written as ^𝑡₀𝐡_𝑚(𝜏)𝑑𝜏≅𝐏_𝑚𝐡_𝑚(𝑡), (3.8) where 𝐏_𝑚 is the 𝑚-square operational matrix of integration which satisfies the following recursive formula [3]: 𝐏_𝑚=𝐏_𝑚/2−1𝐇2𝑚^𝑚/21𝐇2𝑚^−1𝑚/2𝟎_𝑚/2,𝐏₁1=[2], (3.9) where 𝟎_𝑚/2 is an 𝑚/2-square zero matrix. The Haar matrix 𝐇_𝑚 also has the following recursive formula [3]: 𝐇_𝑚=𝐇_𝑚/2⊗𝐈11_𝑚/2⊗1−1,𝐇₁=[1]. (3.10) Particularly, it was proven that the following relationship holds [3]: 𝐇_𝑚⁻¹=1𝑚𝐇^𝑇_𝑚𝐃_𝑚, (3.11) where 𝐃_𝑚=diag(1122⋯2^𝑝−1⋯2^𝑝−1_𝑚/2) and 𝑝=log₂𝑚. This diagonal matrix 𝐃_𝑚 also can be represented in the recursive form 𝐃_𝑚=𝐃_𝑚/2𝟎_𝑚/2𝟎_𝑚/2𝑚2𝐈_𝑚/2,𝐃₁=[1], (3.12) where 𝑚=2^𝑘,𝑘=1,2,…,𝐽, and 𝐽 is called a resolution scale or level.

We present the following lemma which will be used to decompose the generalized Sylvester matrix equation.

Lemma 3.1. Let 𝐇_𝑚 be a Haar matrix defined in (3.10). Then, its inverse matrix has the following recursive form: 𝐇_𝑚⁻¹=𝐇⁻¹_𝑚/2⊗𝐈0.50.5^𝑚/2⊗0.5−0.5. (3.13)

Proof. We assume that 𝐇_𝑚⁻¹ has the following recursive structure: 𝐇_𝑚⁻¹=𝐇⁻¹_𝑚/2⊗𝑎𝑏𝐈_𝑚/2⊗𝑐𝑑, (3.14) where 𝑎,𝑏,𝑐,𝑑 are constants to be determined. Now, we multiply 𝐇_𝑚 and 𝐇_𝑚⁻¹: 𝐇_𝑚𝐇_𝑚⁻¹=(𝐇_𝑚/2⊗11)𝐇⁻¹_𝑚/2⊗𝑎𝑏(𝐇_𝑚/2⊗1111)𝐈_𝑚/2⊗𝑐𝑑𝑐𝑑(𝐈_𝑚/2⊗1−11−1)𝐇⁻¹_𝑚/2⊗𝑎𝑏𝑎𝑏(𝐈_𝑚/2⊗1−11−1)𝐈_𝑚/2⊗𝑐𝑑𝑐𝑑. (3.15) Then, using the property of (𝐀⊗𝐁)(𝐂⊗𝐃)=𝐀𝐂⊗𝐁𝐃, we obtain 𝐇_𝑚𝐇_𝑚⁻¹=𝐈_𝑚/2⊗(𝑎+𝑏)𝐇_𝑚/2𝐇⊗(𝑐+𝑑)_−1𝑚/2⊗(𝑎−𝑏)𝐈_𝑚/2⊗(𝑐−𝑑). (3.16) Thus, 𝑎=𝑏=0.5,𝑐=0.5,𝑑=−0.5 satisfy  𝐇_𝑚𝐇_𝑚⁻¹=𝐇_𝑚⁻¹𝐇_𝑚=𝐈.

4. Singular Linear System

Consider a linear continuous-time singular system described by 𝐄𝐱̇(𝑡)=𝐀𝐱(𝑡)+𝐁𝐮(𝑡),𝐱(0)=𝐱₀, (4.1) where 𝐱(𝑡)∈𝑅^𝑝 denotes the vector of state variables, 𝐮(𝑡)∈𝑅^𝑞 denotes the vector of manipulated inputs, 𝐄, 𝐀 are 𝑝×𝑝 matrices, 𝐄 is generally singular, and 𝐁 is a 𝑝×𝑞 matrix. Without loss of generality, we assume that rank(𝐀)=𝑝 and (4.1) is regular, that is, det(𝜆𝐄−𝐀)≠0. Regularity means that the solution 𝐱(𝑡) is uniquely determined by the given initial value 𝐱₀ and input 𝐮(𝑡).

If the input function vector 𝐮(𝑡) is square integrable in the interval [0,1), then it can be represented in a Haar function basis 𝐡_𝑚(𝑡) as 𝐮(𝑡)=𝐆𝐡_𝑚(𝑡), (4.2) where 𝐆∈𝑅^𝑞×𝑚 is a Haar coefficient matrix and can be obtained by the method described in Section 3. Likewise, 𝐱̇(𝑡) is expanded in Haar function basis 𝐱̇(𝑡)=𝐕𝐡_𝑚(𝑡), (4.3) where 𝐕∈𝑅^𝑝×𝑚 is the unknown matrix to be determined. From the definition of the Haar function, the initial state can be represented as follows: 𝐱₀=𝐱₀0⋯0𝐡_𝑚(𝑡). (4.4) Integrating (4.3) from 0 to 𝑡, we have 𝐱(𝑡)=𝐕𝐏_𝑚𝐡_𝑚(𝑡)+𝐱₀. (4.5) Integrating (4.1) and using (3.8) and (4.4), after canceling 𝐡_𝑚(𝑡), we obtain 𝐄𝐕−𝐀𝐕𝐏_𝑚=𝐐, (4.6) where we define 𝐐≜[𝐀𝐱₀0⋯0]+𝐁𝐆. Thus, the differential matrix equation (4.1) has been transformed to a generalized Sylvester matrix equation that must be solved for 𝐕. Equation (4.6) can be solved by using Kronecker product as in [6] (𝐈_𝑚⊗𝐄+𝐏^T_𝑚⊗𝐀)vec(𝐕)=vec(𝐐), (4.7) where 𝐈_𝑚 is a unit matrix. Equation (4.7) can be solved by LU factorization. However, the coefficient matrix 𝐈_𝑚⊗𝐄+𝐏^T_𝑚⊗𝐀 has dimension 𝑝𝑚×𝑝𝑚, making this approach impractical except for small systems. There are other methods for solving the Sylvester matrix equation (4.6), for example, the Bartels-Stewart algorithm, Krylov subspace method, and matrix sign function method (see [10] and references therein). In [3, 11], recursive algorithms were derived to solve the equations of type 𝐕−𝐀𝐕𝐏_𝑚=𝐐 for linear systems. It should be noted that the algorithm in [3] is not applicable to a generalized Sylvester matrix equation (4.6), since E is a singular matrix.

4.1. Decomposition and Recursive Binary Tree

Under the assumption that 𝐀 is a nonsingular matrix, (4.6) can be written as the following Sylvester equation: 𝐀⁻¹𝐄𝐕−𝐕𝐏_𝑠=𝐀⁻¹𝐐. (4.8) To decompose (4.8), we split 𝐕 and 𝐀⁻¹𝐐 by columns: 𝐀_𝐸𝐕₁⁽¹⁾𝐕₁⁽²⁾−𝐕₁⁽¹⁾𝐕₁⁽²⁾𝐏_𝑚/2−1𝐇2𝑚^𝑚/21𝐇2𝑚^−1𝑚/2𝟎_𝑚/2=𝐐_𝑎𝐐_𝑏, (4.9) where 𝐀_𝐸≜𝐀⁻¹𝐄, 𝐀⁻¹𝐐=[𝐐_𝑎𝐐_𝑏], 𝐕=[𝐕₁⁽¹⁾𝐕₁⁽²⁾] with 𝐐_𝑎,𝐐_𝑏,𝐕₁⁽¹⁾,𝐕₁⁽²⁾∈𝑅^𝑝×𝑚/2. Here 𝐕_𝑘^(𝑟) denotes the matrix that is decomposed at level 𝑘 with 𝑟={2^𝑘,2^𝑘−1}. Then, we obtain the following reduced-order matrix equations: 𝐀_𝐸𝐕₁⁽¹⁾−𝐕₁⁽¹⁾𝐏_𝑚/2−1𝐕2𝑚¹⁽²⁾𝐇⁻¹_𝑚/2=𝐐_𝑎.𝐀_𝐸𝐕₁⁽²⁾+1𝐕2𝑚¹⁽¹⁾𝐇_𝑚/2=𝐐_𝑏. (4.10)(4.11) Since 𝐄 is a singular matrix, 𝐀_𝐸 is also singular. Thus, we postmultiply by 𝐇⁻¹_𝑚/2 both sides of (4.11) to express 𝐕₁⁽¹⁾ in terms of 𝐕₁⁽²⁾𝐕₁⁽¹⁾=−2𝑚𝐀_𝐸𝐕₁⁽²⁾𝐇⁻¹_𝑚/2+2𝑚𝐐_𝑏𝐇⁻¹_𝑚/2. (4.12) Substituting (4.12) into (4.10) yields −2𝑚𝐀²_𝐸𝐕₁⁽²⁾𝐇⁻¹_𝑚/2+2𝑚𝐀_𝐸𝐐_𝑏𝐇⁻¹_𝑚/2+2𝑚𝐀_𝐸𝐕₁⁽²⁾𝐇⁻¹_𝑚/2𝐏_𝑚/2−2𝑚𝐐_𝑏𝐇⁻¹_𝑚/2𝐏_𝑚/2−1𝐕2𝑚¹⁽²⁾𝐇⁻¹_𝑚/2=𝐐_𝑎. (4.13) Therefore, the original problem is decomposed into a reduced-order generalized Sylvester matrix equation (4.13) and a matrix algebraic equation (4.12). Again postmultiplying by 𝐇_𝑚/2 both sides of (4.13), we have(−2𝑚𝐀²_𝐸−12𝑚𝐈)𝐕₁⁽²⁾+2𝑚𝐀_𝐸𝐕₁⁽²⁾𝐇⁻¹_𝑚/2𝐏_𝑚/2𝐇_𝑚/2=𝐐_𝑎𝐇_𝑚/2−2𝑚𝐀_𝐸𝐐_𝑏+2𝑚𝐐_𝑏𝐇⁻¹_𝑚/2𝐏_𝑚/2𝐇_𝑚/2. (4.14) In (4.14), we define 𝐂_𝑚/2≜𝐇⁻¹_𝑚/2𝐏_𝑚/2𝐇_𝑚/2. (4.15) Then, the matrix 𝐂_𝑚/2 is an upper triangular matrix and has the following recursive form: 𝐂_𝑚/2=12𝐂_𝑚/42𝑚𝟏_𝑚/4𝟎_𝑚/412𝐂_𝑚/4,𝐂₁1=[2], (4.16) where 𝟏_𝑚/4 denotes 𝑚/4-square matrix with all elements being 1 (see Appendix ).

Substituting (4.16) into (4.14) and splitting 𝐕₁⁽²⁾ and the right-hand side of (4.14) by columns yields 𝐀_ℎ𝐕₂⁽³⁾𝐕₂⁽⁴⁾+2𝑚𝐀_𝐸𝐕₂⁽³⁾𝐕₂⁽⁴⁾12𝐂_𝑚/42𝑚𝟏_𝑚/4𝟎_𝑚/412𝐂_𝑚/4=𝐓₂⁽³⁾𝐓₂⁽⁴⁾, (4.17) where 𝐀_ℎ≜(−2𝑚𝐀²_𝐸−1𝐐2𝑚𝐈),_𝑎𝐇_𝑚/2−2𝑚𝐀_𝐸𝐐_𝑏+2𝑚𝐐_𝑏𝐂_𝑚/2≜𝐓₁⁽²⁾=𝐓₂⁽³⁾𝐓₂⁽⁴⁾,𝐕₁⁽²⁾=𝐕₂⁽³⁾𝐕₂⁽⁴⁾. (4.18) Thus, (4.17) is decomposed into two matrix equations with dependent and independent subsystems. 𝐀_ℎ𝐕₂⁽³⁾+𝑚𝐀_𝐸𝐕₂⁽³⁾𝐂_𝑚/4=𝐓₂⁽³⁾.𝐀_ℎ𝐕₂⁽⁴⁾+𝑚𝐀_𝐸𝐕₂⁽⁴⁾𝐂_𝑚/4=𝐓₂⁽⁴⁾−4𝐀_𝐸𝐕₂⁽³⁾𝟏_𝑚/4. (4.19)(4.20) In (4.19) and (4.20), we first solve for 𝐕₂⁽³⁾ and then after updating the right-hand side of (4.20) with respect to 𝐕₂⁽³⁾, solve for 𝐕₂⁽⁴⁾. Since (4.19) and (4.20) have the same form as (4.17) and 𝐂_𝑚/4 is still an upper triangular matrix, they can be decomposed into two subsystems in which the dimension has been reduced by half, respectively. Therefore, we recursively decompose each equation into two equations until no further decomposition is possible in which all 𝐕_𝐽^(𝑟),𝐓_𝐽^(𝑟)(𝑟=2^𝐽−1+1,…,2^𝐽) are column vectors. This procedure constructs the binary tree as shown in Figure 1.

Figure 1 

Binary tree for resolution scale 𝐽=4.

A binary tree is a rooted tree in which each node has at most two children, designated as a left child and a right child. A full binary tree is a binary tree in which each node has exactly two children or none. A perfect (or complete) binary tree is a full binary tree in which all leaves have the same depth [12]. In Figure 1, the binary tree in the dotted box is a perfect binary tree of depth 𝐽−1. An external node (or leaf node) is a node with no children. For instance, the nodes labeled 1, 9, 10, 11, 12, 13, 14, 15, and 16 in Figure 1 are external nodes.

Matrix equations corresponding to all external nodes of the perfect binary tree are classified into two types of equations described as follows: 𝐀_ℎ𝐕_𝐽^(𝑟)+4𝐀_𝐸𝐕_𝐽^(𝑟)𝐂₁=𝐓_𝐽^(𝑟),𝑟=2^𝐽−1+1,2^𝐽−1+3,…,2^𝐽𝐀−1(𝑟isodd),_ℎ𝐕_𝐽^(𝑟)+4𝐀_𝐸𝐕_𝐽^(𝑟)𝐂₁=𝐓_𝐽^(𝑟)−4𝐀_𝐸𝐕_𝐽^(𝑟−1)𝟏₁,𝑟=2^𝐽−1+2,2^𝐽−1+4,…,2^𝐽(𝑟iseven). (4.21) Note that in equation (4.21), 𝐂₁=1/2, 𝟏₁=1. Thus, they become simple linear matrix equations as follows: (𝐀_ℎ+2𝐀_𝐸)𝐕_𝐽^(𝑟)=𝐓_𝐽^(𝑟),if𝑟isodd,(𝐀_ℎ+2𝐀_𝐸)𝐕_𝐽^(𝑟)=𝐓_𝐽^(𝑟)−4𝐀_𝐸𝐕_𝐽^(𝑟−1),if𝑟iseven. (4.22)

4.2. Combined Preorder and Postorder Traversal Algorithm

Visiting all the nodes in a tree in some particular order is known as a tree traversal. A preorder traversal visits the root of a subtree, then the left and right subtrees recursively. A postorder traversal visits the left and right subtrees recursively, then the root node of the subtree [12]. For example, the preorder and postorder traversals of the binary tree shown in Figure 1 are as follows:

During the decomposition of (4.14), the right-hand side of the right child is split after updating it recursively as follows: 𝐓_𝑘^(𝑟)−4𝐀_𝐸𝐕_𝑘^(𝑟−1)𝟏_𝑚/2^𝑘=𝐓^(2𝑟−1)_𝑘+1𝐓^(2𝑟)_𝑘+1. (4.23) This splitting and updating sequence is a preorder traversal of the perfect binary tree from root node ②. The unknown matrix 𝐕₁⁽²⁾ is obtained by merging all column vectors 𝐕_𝐽^(𝑟)(𝑟=2^𝐽−1+1,…,2^𝐽). This sequence is a postorder traversal of the perfect binary tree from root node ②. To update (4.23), we need 𝐕_𝑘^(𝑟−1) which is obtained from the left child. Hence, to solve (4.22), it is necessary to update, split, and solve by using the following combined preorder and postorder traversal method.

The pseudocode of the proposed algorithm is as in Algorithm 1.

Step 1. 𝙸𝚗𝚒𝚝𝚒𝚊𝚕𝚒𝚣𝚎𝐀ℎ,𝐀𝐸,𝐓. Step 2. 𝙾𝚋𝚝𝚊𝚒𝚗𝐕1(2) 𝙸𝚗𝚙𝚞𝚝:𝚁𝚎𝚜𝚘𝚕𝚞𝚝𝚒𝚘𝚗𝚜𝚌𝚊𝚕𝚎𝐽 𝑊𝑎𝑣𝑒𝑆𝑜𝑙𝑣𝑒𝑟(𝐽) { for(𝑟=2𝐽−1+1;𝑟<2𝐽;𝑟=𝑟+2) { Solvefor𝐕𝐽(𝑟)thesystem(𝐀ℎ+2𝐀𝐸)𝐕𝐽(𝑟)=𝐓𝐽(𝑟) Update𝐓𝐽(𝑟+1)accrodingto𝐓𝐽(𝑟+1)=𝐓𝐽(𝑟+1)−4𝐀𝐸𝐕𝐽(𝑟) Solvefor𝐕𝐽(𝑟+1)thesystem(𝐀ℎ+2𝐀𝐸)𝐕𝐽(𝑟+1)=𝐓𝐽(𝑟+1) 𝑊𝑎𝑣𝑒𝑇𝑟𝑒𝑒(1,𝐽,𝑟+1) } } 𝙸𝚗𝚙𝚞𝚝:𝚛𝚗𝚘𝚒𝚜𝚊𝚗𝚞𝚖𝚋𝚎𝚛𝚘𝚏𝚛𝚎𝚌𝚞𝚛𝚜𝚒𝚟𝚎𝚌𝚊𝚕𝚕. 𝚁𝚎𝚜𝚘𝚕𝚞𝚝𝚒𝚘𝚗𝚜𝚌𝚊𝚕𝚎𝐽 𝑟𝚒𝚜𝚊𝚗𝚘𝚍𝚎𝚗𝚞𝚖𝚋𝚎𝚛. 𝑊𝑎𝑣𝑒𝑇𝑟𝑒𝑒(rno,𝐽,𝑟) { if(𝐽−rno≤0) return Merge:𝐕(𝑟/2)𝐽−rno=𝐕(𝑟−1)𝐽−rno+1𝐕(𝑟)𝐽−rno+1 𝑟if(2iseven) 𝑟𝑊𝑎𝑣𝑒𝑇𝑟𝑒𝑒(rno+1,𝐽,2); else UpdateandSplit:𝐓(𝑟−1)𝐽−rno+1𝐓(𝑟)𝐽−rno+1=𝐓(𝑟/2+1)𝐽−rno−4𝐀𝐸𝐕(𝑟/2)𝐽−rno𝟏𝑚/2𝐽−rno } Step 3. 𝚂𝚘𝚕𝚟𝚎𝐕1(1)𝚏𝚛𝚘𝚖(4.12).

Algorithm 1 

For example, at resolution scale 𝐽=4, the proposed combined preorder and postorder traversal method is illustrated in Figure 2.

Figure 2 

The combined preorder and postorder traversal for resolution scale 𝐽=4.

In Figure 2, nodes 2, 3, 5, and 9 of preorder traversal are done at Step 1 and the remaining nodes are processed at Step 2. The computational efficiency of the proposed method is discussed in the next section.

5. An Illustrated Example

In this section, an example is presented to illustrate the proposed algorithm. We consider a singular linear system of (4.1) with 00𝐄=1000010000000001,𝐀=−3301.000101.00621.4−28.2700−327.112.721,𝐁=52.65−23.69, (5.1) and 𝐗₀=[00.51.00]^T. And we assume that 𝑢(𝑡) is a unit step function. In the cases of 𝐽=4 and 8, the simulation results are depicted in Figures 3 and 4, respectively.

Figure 3 

Case for resolution scale 𝐽=4.

Figure 4 

Case for resolution scale 𝐽=8.

From these figures, it is clear that the solution accuracy is improved when the resolution scale is increased. However, it requires more computational time.

In (4.7), the LU factorization of 𝐈_𝑚⊗𝐄+𝐏^T_𝑚⊗𝐀 involves 𝑂(𝑚³𝑝³) flops. The cost of the proposed algorithm is the sum of the cost of WinSolver, 𝑂((𝑚/2)𝑝³+(𝑚/2)𝑝²), and the cost of WinTree, 𝑂(∑^𝐽−2_𝑘=1(2^𝐽−1𝑝²+(2^{𝐽+2𝑘−2}−2^𝐽−2)𝑝)) (see Appendix ). Since 𝑚=2^𝐽, the costs of WinSolver and 𝑂(𝑚³𝑝³) can be rewritten as 𝑂(2^𝐽−1𝑝³+2^𝐽−1𝑝²) and 𝑂(2^3𝐽𝑝³), respectively. Thus, the total cost of the proposed algorithm is 𝑂(2^𝐽−1𝑝³+2^𝐽−1𝑝²+^𝐽−2_𝑘=1(2^𝐽−1𝑝²+(2^{𝐽+2𝑘−2}−2^𝐽−2)𝑝))flops. (5.2) Table 1 and Figure 5 show that the computational cost of the proposed algorithm is significantly less than the Kronecker product method, and that the flop counts are increasing rapidly with resolution scale. As the resolution scale grows, the flop counts of WinTree is increasing more rapidly than that of WinSolver since the sizes of matrices 𝟏_𝑚,𝐓_𝑚, and 𝐕_𝑚 increase exponentially.

Table 1 

Flop counts for various sizes of the matrix A and resolution scales.

Figure 5 

Log plot of flop counts for the Kronecker product method and the proposed method.

6. Conclusions

An efficient computational method was presented for state space analysis of singular systems via Haar wavelets. The problem was formulated as a generalized Sylvester matrix equation. We presented an explicit expression for the inverse of the Haar matrix and a combined preorder and postorder traversal algorithm to solve the problem more effectively. The full-order generalized Sylvester matrix equation was solved in terms of the solutions of simple linear matrix equations by the proposed algorithm. The efficiency of the proposed method was demonstrated by a numerical example.

Appendices

A. Formula for 𝐂_𝑚

In this appendix, we derive a formula for 𝐂_𝑚. By using (3.13), (3.9), and (3.10), we can write 𝐂_𝑚=𝐇_𝑚⁻¹𝐏_𝑚𝐇_𝑚=𝐇⁻¹_𝑚/2⊗𝐈0.50.5^𝑚/2⊗0.5𝐏−0.5^𝑚/2−1𝐇2𝑚^𝑚/21𝐇2𝑚^−1𝑚/2𝟎_𝑚/2𝐇_𝑚/2⊗𝐈11_𝑚/2⊗=𝐇1−1[(⁻¹_𝑚/2⊗𝐏0.50.5^𝑚/2−1𝐇2𝑚^𝑚/2+𝐈_𝑚/2⊗0.51−0.5𝐇2𝑚^−1𝑚/2𝟎_𝑚/2×𝐇_𝑚/2⊗𝐈11_𝑚/2⊗=𝐇1−1_−1𝑚/2⊗𝐏0.50.5^𝑚/2−1𝐇2𝑚^−1𝑚/2⊗𝐇0.50.5^𝑚/2+1𝐈2𝑚^𝑚/2⊗𝐇0.5−0.5^−1𝑚/2𝟎_𝑚/2𝐇_𝑚/2⊗𝐈11_𝑚/2⊗=𝐇1−1_−1𝑚/2⊗𝐏0.50.5^𝑚/2−1𝐇2𝑚^−1𝑚/2⊗𝐇0.50.5^𝑚/2𝐇_𝑚/2⊗𝐈11_𝑚/2⊗+11−1𝐈2𝑚^𝑚/2⊗𝐇0.5−0.5^−1𝑚/2𝟎_𝑚/2𝐇_𝑚/2⊗𝐈11_𝑚/2⊗=𝐇1−1_−1𝑚/2⊗𝐏0.50.5^𝑚/2(𝐇_𝑚/2⊗111)−𝐇2𝑚^−1𝑚/2⊗𝐇0.50.5^𝑚/2(𝐈_𝑚/2⊗)+11−1𝐈2𝑚^𝑚/2⊗𝐇0.5−0.5^−1𝑚/2(𝐇_𝑚/2⊗11). (A.1) Since (𝐀⊗𝐁)𝐂=(𝐀⊗𝐁)(𝐂⊗1)=(𝐀𝐂⊗𝐁) and (𝐀⊗𝐁)(𝐂⊗𝐃)=(𝐀𝐂⊗𝐁𝐃), the above equation is rewritten as 𝐂_𝑚=(𝐇⁻¹_𝑚/2𝐏_𝑚/2)⊗0.50.5(𝐇_𝑚/2⊗111)−2𝑚(𝐇⁻¹_𝑚/2𝐇_𝑚/2)⊗0.50.5(𝐈_𝑚/2⊗)+11−12𝑚(𝐈_𝑚/2𝐇⁻¹_𝑚/2)⊗0.5−0.50.5−0.5(𝐇_𝑚/2⊗11)=(𝐇⁻¹_𝑚/2𝐏_𝑚/2𝐇_𝑚/2−1)⊗0.50.511112𝑚(𝐈_𝑚/2+1)⊗0.50.51−11−12𝑚(𝐈_𝑚/2)⊗0.5−0.511=𝐂_𝑚/2⊗12121212+1𝐈2𝑚^𝑚/2⊗+−0.50.5−0.50.50.50.5−0.5−0.5=𝐂_𝑚/2⊗12121212+1𝐈2𝑚^𝑚/2⊗=101−102𝐂_𝑚/21𝑚𝟏_𝑚/2𝟎_𝑚/212𝐂_𝑚/2. (A.2)

B. Flop Counts of the Combined Preorder and Postorder Traversal Algorithm

In this appendix, we show that the computational cost for the combined preorder and postorder traversal algorithm described in Section 4.2 can be obtained as follows:

(1) WinSolve
𝚂𝚘𝚕𝚟𝚎𝚏𝚘𝚛𝐕^(𝑟)_𝐽𝚝𝚑𝚎𝚜𝚢𝚜𝚝𝚎𝚖(𝐀_ℎ+2𝐀_𝐸)𝐕_𝐽^(𝑟)=𝐓_𝐽^(𝑟)∶𝑂(𝑝³).𝚄𝚙𝚍𝚊𝚝𝚎𝐓_𝐽^(𝑟+1)𝚊𝚌𝚌𝚘𝚛𝚍𝚒𝚗𝚐𝚝𝚘𝐓_𝐽^(𝑟+1)=𝐓_𝐽^(𝑟+1)−4𝐀_𝐸𝐕_𝐽^(𝑟)∶𝑂(𝑝(2𝑝−1)+𝑝)=𝑂(2𝑝²).𝚂𝚘𝚕𝚟𝚎𝚏𝚘𝚛𝐕_𝐽^(𝑟+1)𝚝𝚑𝚎𝚜𝚢𝚜𝚝𝚎𝚖(𝐀_ℎ+2𝐀_𝐸)𝐕_𝐽^(𝑟+1)=𝐓_𝐽^(𝑟+1)∶𝑂(𝑝³). The total iteration number of “for(𝑟=2^𝐽−1+1;𝑟<2^𝐽;𝑟=𝑟+2)” is 𝑚/4. Thus, WinSolve involves 𝑂((𝑚/4)(𝑝³+2𝑝²+𝑝³))=𝑂((𝑚/2)(𝑝³+p²))flops.

(2) WinTree
𝚄𝚙𝚍𝚊𝚝𝚎𝚊𝚗𝚍𝚜𝚙𝚕𝚒𝚝𝐓^(𝑟/2+1)_𝐽−rno−4𝐀_𝐸𝐕^(𝑟/2)_𝐽−rno𝟏_𝑚/2^𝐽−rno(see Table 2).

Table 2 

Therefore, the computational cost for WinTree can be calculated by 𝑂(^𝐽−2_𝑘=1(2^𝑘+1𝑝²+2^𝑘(2^2𝑘−1)𝑝)×2^{𝐽−𝑘−2})=𝑂(^𝐽−2_𝑘=1(2^𝐽−1𝑝²+(2^{𝐽+2𝑘−2}−2^𝐽−2)𝑝)). (A.1)