Refinement Methods for State Estimation via Sylvester-Observer Equation

Abstract

We present new iterative methods based on refinement process for solving large sparse Sylvester-observer equations applied in state estimation of a continuous-time system. These methods use projection methods to produce low-dimensional Sylvester-observer matrix equations that are solved by the direct methods. Moreover, the refinement process described in this paper has the capability of improving the results obtained by any other methods. Some numerical results will be reported to illustrate the efficiency of the proposed methods.

1. Introduction

Consider the continuous-time linear system 𝑥̇(𝑡)(𝑡)(𝑡),𝑦(𝑡)(𝑡),=𝐴𝑥+𝐵𝑢=𝐶𝑥 (1.1) where 𝐴∈𝑅^𝑛×𝑛, 𝐵∈𝑅^𝑛×𝑚, and 𝐶∈𝑅^𝑟×𝑛.

We well know that all the state-feedback problems, such as the feedback stabilization, the LQR, and the state-feedback 𝐻_∞-control problems, require the state vector 𝑥(𝑡) explicitly [1]. In most practical situations, the states are not fully accessible and the designer only knows the output and the input vectors. The unavailable states somehow need to be estimated accurately from the knowledge of the matrices 𝐴, 𝐵, and 𝐶, the output vector 𝑦(𝑡), and the input vector 𝑢(𝑡). There are two common procedures for state estimation: one via eigenvalue assignment (EVA) and the other via solution of the Sylvester-observer equation. The main step in state estimation via solution of the Sylvester-observer equation is solving the Sylvester-observer equation of the form 𝑋𝐴−𝐹𝑋=𝐺𝐶, (1.2) where 𝑋, 𝐴, 𝐹∈𝑅^𝑛×𝑛, 𝐺∈𝑅^𝑛×𝑟, and 𝐶∈𝑅^𝑟×𝑛. Sylvester-observer equations (1.2) play vital roles in a number of applications such as control and communications theory [1], model reduction [2–4], numerical solution of matrix differential Riccati equations [5], and image processing [6].

The analytical solution of the matrix equation (1.2) has been considered by many authors; see [7, 8]. Direct methods for solving the matrix equation (1.2) are attractive if the matrices are of small size. These methods are based on the Schur decomposition, by which the original equation is transformed into a form that is easy to be solved by a forward substitution. Moreover, the matrices of the large practical problem are very sparse. Since, the standard methods for solving the Sylvester equations destroy the sparsity of the problems, they are only applicable for the matrices of small size; see [9–11]. Iterative projection methods for solving large Sylvester-observer matrix equations have been developed during the past years; see [12–16]. These methods use the Arnoldi process to compute an orthonormal basis of certain Krylov subspace. In this paper, we extend the idea to propose a new projection method for solving (1.1) based on Weighted block Krylov subspace process. These methods are based on the reduction of the large sparse Sylvester-observer equation to the low-dimensional problem by orthogonal similarity that is solved by the direct methods. Moreover, refinement process presented in Sections 4 and 5 has the capability of improving the results obtained by any other methods.

The remainder of the paper is organized as follows. In Section 2, we describe some fundamental results about control theory. Then, we discuss how the states of a continuous-time system can be estimated in Section 3. In Sections 4 and 5, we introduce two new iterative methods for solving large sparse Sylvester-observer equation. These methods are based on the reduction of the large sparse Sylvester-observer equation to the low-dimension problem. Section 6 is devoted to some numerical tests. Some concluding remarks are given in Section 7.

2. Some Fundamental Results

The two basic concepts in control theory are controllability and observability of a control system. In this section, we will state some well-known facts about controllability and observability for convenient use later in the paper. For an excellent account of controllability and observability and related results, see [1].

Definition 2.1. The system (1.1) is said to be controllable if, starting from any initial state 𝑥(0), the system can be driven to any final state 𝑥_𝑓 in some finite time 𝑡_𝑓, choosing the input vector 𝑢(𝑡), 0≤𝑡≤𝑡_𝑓, appropriately.

Observability is the dual concept of controllability.

Definition 2.2. The system (1.1) is said to be observable if there exists 𝑡₁>0 such that the initial state 𝑥(0) can be uniquely determined from the knowledge of 𝑢(𝑡) and 𝑦(𝑡) for all0≤𝑡≤𝑡₁.

Remark 2.3. Since matrix 𝐶 does not have any role in the definition of controllability, the controllability of the system (1.1) (see [1]), is often referred to as the controllability of the pair (𝐴,𝐵). Similarly, since matrix 𝐵 does not have any role in the definition of observability, the observability of the system (1.1) is often referred to as the observability of the pair (𝐴,𝐶).

Some well-known criteria of controllability and observability are now stated in the next two theorems. The proofs of Theorems 2.4 and 2.5 can be found in [1]. In the following, 𝐴 is 𝑛×𝑛, 𝐵 is 𝑛×𝑚(𝑚≤𝑛), and 𝐶 is 𝑟×𝑛(𝑟≤𝑛).

Theorem 2.4. The pair (𝐴,𝐵) is controllable if and only if any one of the following equivalent conditions holds.
(1) The controllability matrix 𝐶_𝑀=𝐵,𝐴𝐵,𝐴²𝐵,…,𝐴^𝑛−1𝐵, (2.1) has rank 𝑛. (2) Rank (𝐴−𝜆𝐼,𝐵)=𝑛 for every eigenvalue 𝜆 of 𝐴. (3) Let (𝜆,𝑥) be an eigenpair of 𝐴, that is, 𝐴𝑦=𝜆𝑦; then 𝐶𝑦≠0.

Proof. See [1].

Theorem 2.5. The pair (𝐴,𝐶) is observable if and only if any one of the following equivalent conditions holds
(1) The observability matrix 𝑂_𝑀=𝐶𝐶𝐴𝐶𝐴²⋮𝐶𝐴^𝑛−1 (2.2) has rank 𝑛. (2) The matrix 𝐶𝜆𝐼−𝐴 has rank 𝑛 for every eigenvalue 𝜆 of 𝐴. (3) Let (𝜆,𝑦) be an eigenpair of 𝐴, that is, 𝐴𝑦=𝜆𝑦; then 𝐶𝑦≠0.

Proof. See [1].

Definition 2.6. A matrix 𝐴 is called a stable matrix if all of the eigenvalues of 𝐴 have negative real parts.

3. State Estimation

In this section, we discuss how the states of a continuous-time linear system (1.1) can be estimated. The discussions here apply equally to the discrete-time systems possibly with some minor changes. So we concentrate on the continuous-time case only. We describe two common procedures for state estimation: one via eigenvalue assignment (EVA) and the other,via solution of the Sylvester-observer equation.

3.1. State Estimation via Eigenvalue Assignment

Now consider the linear time-invariant system (1.1). Let 𝑥(𝑡) be an estimate of the state vector 𝑥(𝑡). Obviously, we would like to construct the vector 𝑥(𝑡) in such a way that the error 𝑒(𝑡)=𝑥(𝑡)−𝑥(𝑡) approaches zero as fast as possible, for all initial states 𝑥(0) and for every input 𝑢(𝑡). The following theorem shows that the problem of state estimation can be solved by finding a matrix 𝐾 such that the matrix 𝐴−𝐾𝐶 has a suitable desired spectrum. See [1].

Theorem 3.1. If (𝐴,𝐶) is observable, then the states x(𝑡) of the system (1.1) can be estimated by 𝑥̇(𝑡)=(𝐴−𝐾𝐶)𝑥(𝑡)+𝐾𝑦(𝑡)+𝐵𝑢(𝑡), (3.1) where 𝐾 is constructed such that 𝐴−𝐾𝐶 is a stable matrix. The error 𝑒(𝑡)=𝑥(𝑡)−𝑥(𝑡) is governed by 𝑒̇(𝑡)=(𝐴−𝐾𝐶)𝑒(𝑡), (3.2) and 𝑒(𝑡)→0 as 𝑡→∞.

Proof. See [1].

3.2. State Estimation via Sylvester Equation

There is another approach for state estimation. Knowing 𝐴, 𝐵, 𝐶, 𝑢(𝑡), and 𝑦(𝑡), let us construct the system 𝑧̇(𝑡)=𝐹𝑧(𝑡)+𝐺𝑦(𝑡)+𝑃𝑢(𝑡), (3.3) where 𝐹 is 𝑛×𝑛, 𝐺 is 𝑛×𝑟, and 𝑃 is 𝑛×𝑚, in such a way that for some constant 𝑛×𝑛 nonsingular matrix 𝑋 the error vector 𝑒(𝑡)=𝑧(𝑡)−𝑋𝑥(𝑡)→0 for all 𝑥(0), 𝑧(0), and for every input 𝑢(𝑡). The vector 𝑧(𝑡) will then be an estimate of 𝑋𝑥(𝑡). The system (3.3) is then said to be the state observer for the system (1.1). D. Luenberger originated the idea and is hence referred to in control theory as the Luenberger observer; see [1].

We now show that the system (3.3) will be a state observer if the matrices 𝑋, 𝐹, 𝐺, and 𝑃 satisfy certain requirements.

Theorem 3.2. The system (3.3) is a state observer of the system (1.1), that is, 𝑧(𝑡) is an estimate of 𝑋𝑥(𝑡) in the sense that the error 𝑒(𝑡)=𝑧(𝑡)−𝑋𝑥(𝑡)→0 as 𝑡→∞ for any initial conditions 𝑥(0), 𝑧(0), and 𝑢(𝑡) if
(1)  𝑋𝐴−𝐹𝑋=𝐺𝐶,
(2)  𝑃=𝑋𝐵,
(3)  𝐹 is stable.

Proof. See [1].

Definition 3.3. The matrix equation 𝑋𝐴−𝐹𝑋=𝐺𝐶, (3.4) where 𝑋, 𝐴, 𝐹∈𝑅^𝑛×𝑛, 𝐺∈𝑅^𝑛×𝑟, and 𝐶∈𝑅^𝑟×𝑛, is called the Sylvester-observer equation.

Theorem 3.2 suggests the following method for the observer design; see [1].

Algorithm 3.4 (full-order observer design via Sylvester-observer equation). One has the following.
Inputs
The system matrices 𝐴, 𝐵, and 𝐶 of order 𝑛×𝑛, 𝑛×𝑚, and 𝑟×𝑛, respectively.Output
An estimate 𝑥(𝑡) of the state vector 𝑥(𝑡).Assumptions. (𝐴,𝐶) is observable.Step 1. Find a nonsingular solution of the Sylvester-observer equation (1.2) by choosing 𝐹 as a stable matrix and choosing 𝐺 in such a way that the resulting solution 𝑋 is nonsingular.Step 2. Compute 𝑃=𝑋𝐵.Step 3. Construct the observer 𝑧(𝑡) by solving the system of differential equations 𝑧(𝑡)=𝐹𝑧(𝑡)+𝐺𝑦(𝑡)+𝑃𝑢(𝑡),𝑧(0)=𝑧₀. (3.5)Step 4. Find an estimate 𝑥(𝑡) of 𝑥(𝑡): 𝑥(𝑡)=𝑋⁻¹𝑧(𝑡).

3.3. Characterization of Nonsingular Solutions of the Sylvester Equation

In this section, we describe some necessary conditions for a unique solution of the Sylvester equation to have such properties. The following theorem was proved by Bhattacharyya and de Souza. The proof here has been taken from [1].

Theorem 3.5. Let 𝐴, 𝐹, 𝐺, and 𝐶, respectively, be of order 𝑛×𝑛, 𝑛×𝑛, 𝑛×𝑟, and 𝑟×𝑛. Let 𝑋 be a unique solution of the Sylvester-observer equation (1.2). Then, necessary conditions for 𝑋 to be nonsingular are that (𝐴,𝐶) is observable and (𝐹,𝐺) is controllable.

Proof. See [1].

Corollary 3.6. If 𝐺 is 𝑛×1 and 𝐶 is 1×𝑛, then necessary and sufficient conditions for the unique solution 𝑋 of (1.2) to be nonsingular are that (𝐹,𝐺) is controllable and (𝐴,𝐶) is observable.

Proof. See [1].

Remark 3.7. According to Theorem 3.5 and Corollary 3.6, the controllability of (𝐹,𝐺) and observability of (𝐴,𝐶) guarantee the existence of a nonsingular solution of the Sylvester-observer equation (1.2) in Step 1 of Algorithm 3.4. Moreover, there are other choices for 𝐹 and 𝐺 in Step 1 of Algorithm 3.4 provided (𝐹,𝐺) is controllable. Also, we can use Theorems 2.4 and 2.5 for analyzing the controllability of (𝐹,𝐺) and the observability of (𝐴,𝐶).

Example 3.8. In this example we show how the Sylvester-observer equation (1.2) can be applied in state estimation of a continuous time system (1.1). In this sense, at first we use the MATLAB function ode23 for directly solving (1.1) with 𝑢(𝑡) as the unit step function and 𝑥(0)=(6,0)^𝑇. Then, we apply Algorithm 3.4 for computing the estimate 𝑥(𝑡). Also, the differential equation (3.12) was solved with 𝑧(0)=0 and MATLAB function ode23. The comparison of the state 𝑥(𝑡) and estimate 𝑥(𝑡) for this example is shown in Figure 1. The solid line corresponds to the exact state and the dotted line corresponds to the estimate state. Let 10𝐴=1111,𝐵=,𝐶=10. (3.6) According to criteria 1 of Theorem 2.5, (𝐴,𝐶) is observable. Thus, we can use Algorithm 3.4 for state estimation of (1.1).

(a)

(b)

Figure 1 

The (a) first and (b) second variables of the state 𝑥(𝑡) and estimate 𝑥(𝑡) for Example 3.8.

Step 5. Choose 13𝐺=,𝐹=diag(−1,−3). (3.7) According to criteria 1 of Theorem  2.4, (𝐹,𝐺) is controllable. Thus, by Corollary 3.6, the nonsingular solution 𝑋 of 𝑋𝐴−𝐹𝑋=𝐺𝐶 is 𝑋=0.6667−0.33330.8000−0.2000, (3.8) (computed by MATLAB function lyap).

Step 6. One has𝑃=𝑋𝐵=0.66670.8000. (3.9)

Step 7. An estimate 𝑥(𝑡) of 𝑥(𝑡) is 𝑥(𝑡)=𝑋⁻¹𝑧(𝑡)=−1.52.5𝑧−65¹(𝑡)𝑧₂(𝑡)=−1.5𝑧₁+2.5𝑧₂−6𝑧₁+5𝑧₂, (3.10) where 𝑧(𝑡)=𝑧₁(𝑡)𝑧₂(𝑡), (3.11) is given by 𝑧(𝑡)=𝑧(𝑡)+13𝑦(𝑡)+𝑢(𝑡)(0)−100−30.66670.8000,𝑧=𝑧₀. (3.12)

Remark 3.9. According to Algorithm 3.4 the most important step is Step 1. In the case that 𝑛 is small, there are many reliable algorithm s for solving the Sylvester-observer equation and the states of a continuous-time system can be estimated. However, for large and sparse systems solving the Sylvester-observer equation by the available methods can be costly. In Sections 4 and 5, we introduce two iterative refinement methods for solving large sparse Sylvester-observer equations.

4. Block Refinement Method

As we already mentioned, so far many numerical methods have been developed by different authors; see [1, 13, 17]. For example, the Hessenberg-Schur method is now widely used as an effective computational method for the Sylvester-observer equation. But numerical stability of this method has not been investigated. As the iterative methods are very efficient for the solution of computational problems, we thought it will be good idea to create an iterative method for solving the Sylvester-observer equation 𝑋𝐴−𝐹𝑋=𝐺𝐶 where 𝐴, 𝐹, 𝑋∈𝑅^𝑛×𝑛, 𝐺∈𝑅^𝑛×𝑟, and 𝐶∈𝑅^𝑟×𝑛. In this section we propose to show that the obtained approximate solution of the Sylvester-observer equation by any method can be improved, in other words the accuracy can be increased. If this idea is applicable, then we have found an iterative method for solving of the Sylvester-observer equation.

Theorem 4.1. Let 𝑋₀ be the approximate solution obtained by an arbitrary method for the matrix equation (1.2), and let 𝑅(𝑋₀)=𝐺𝐶−(𝑋₀𝐴−𝐹𝑋₀) be the corresponding residual. If m steps of the block Arnoldi process for matrices 𝐴^𝑇 and 𝐹 have been run (𝑚≪𝑛) and 𝑌_𝑚∈𝑅^𝑚×𝑚 is the solution of the low-dimensional Sylvester equation 𝑌_𝑚𝐴^𝑇_𝑚−𝐹_𝑚𝑌_𝑚=𝑊^𝑇_𝑚𝑅𝑋₀𝑉_𝑚, (4.1) then 𝑅𝑋₁=𝑅𝑋₀+𝑊_𝑚𝑌_𝑚𝑉^𝑇_𝑚=0. (4.2)

Proof. Let 𝑉_𝑚=[𝑣₁,…,𝑣_𝑚] and 𝑊_𝑚=[𝑤₁,…,𝑤_𝑚] be two orthogonal bases constructed by the block Arnoldi process for the matrices 𝐴^𝑇 and 𝐹, respectively. Thus, we have 𝑉^𝑇_𝑚𝑉_𝑚=𝐼_𝑚,𝑊^𝑇_𝑚𝑊_𝑚=𝐼_𝑚. (4.3) Also, the square block Hessenberg matrices 𝐴_𝑚 and 𝐹_𝑚 (𝑚=𝑟ℓ, where 𝑟 and ℓ are the dimensions of blocks) whose nonzero entries are the scalars 𝑎_𝑖𝑗 and 𝑓_𝑖𝑗, constructed by the Block Arnoldi process, can be expressed as 𝐴_𝑚=𝑉^𝑇_𝑚𝐴^𝑇𝑉_𝑚,𝐹_𝑚=𝑊^𝑇_𝑚𝐹𝑊_𝑚. (4.4)
If we set 𝑋₁=𝑋₀+𝑊_𝑚𝑌_𝑚𝑉^𝑇_𝑚, (4.5) then the corresponding residual 𝑅(𝑋₁)=𝐺𝐶−(𝑋₁𝐴−𝐹𝑋₁) satisfies 𝑅𝑋₁=𝐺𝐶−𝑋₀+𝑊_𝑚𝑌_𝑚𝑉^𝑇_𝑚𝐴−𝐹𝑋₀+𝑊_𝑚𝑌_𝑚𝑉^𝑇_𝑚,=𝑅𝑋₀−𝑊_𝑚𝑌_𝑚𝑉^𝑇_𝑚𝐴+𝐹𝑊_𝑚𝑌_𝑚𝑉^𝑇_𝑚,=𝑅𝑋₀−𝑊_𝑚𝑌_𝑚𝐴^𝑇_𝑚𝑉^𝑇_𝑚+𝑊_𝑚𝐹_𝑚𝑌_𝑚𝑉^𝑇_𝑚,=𝑅𝑋₀−𝑊_𝑚𝑌_𝑚𝐴^𝑇_𝑚−𝐹_𝑚𝑌_𝑚𝑉^𝑇_𝑚. (4.6) Since 𝑌_𝑚 is the solution of (4.1) and by using (4.3), we have 𝑅𝑋₁=𝑅𝑋₀−𝑊_𝑚𝑊^𝑇_𝑚𝑅₀𝑉_𝑚𝑉^𝑇_𝑚=0. (4.7)

According to Theorem 4.1, we can develop an iterative method for solving the Sylvester-observer equation when the matrices 𝐴, 𝐹, 𝐺, and 𝐶 are large and sparse. For achieving this idea, if we choose 𝑚<𝑛, then instead of solving 𝑋𝐴−𝐹𝑋=𝐺𝐶 we can solve (4.1). In other words, in this method, first we transform the initial Sylvester-obsever equation to another Sylvester equation with less dimensions and then in each iteration step solve this matrix equation and extend the obtained solution to the solution of the initial equation by (4.5). The algorithm is as follows.

Algorithm 4.2 (block refinement method). (1) Start: choose an initial approximate solution 𝑋₀, and a tolerance 𝜀.
(2) Select two numbers 𝑟 and ℓ for dimensions of block and set 𝑚=𝑟ℓ(𝑚<𝑛).
(3) Compute 𝑅(𝑋₀)=𝐺𝐶−(𝑋₀𝐴−𝐹𝑋₀).
(4) Construct the orthonormal bases 𝑉_𝑚 and 𝑊_𝑚∈𝑅^𝑛×𝑚 by the block Arnoldi process, such that 𝐴_𝑚=𝑉^𝑇_𝑚𝐴^𝑇𝑉_𝑚,𝐹_𝑚=𝑊^𝑇_𝑚𝐹𝑊_𝑚. (4.8)
(5) Solve the low-dimensional Sylvester-observer equation 𝑌_𝑚𝐴^𝑇_𝑚−𝐹_𝑚𝑌_𝑚=𝑊^𝑇_𝑚𝑅(𝑋₀)𝑉_𝑚 with the Hessenberg-Schur method.
(6) Set 𝑋₁=𝑋₀+𝑊_𝑚𝑌_𝑚𝑉^𝑇_𝑚.
(7) Compute 𝑅(𝑋₁)=𝐺𝐶−(𝑋₁𝐴−𝐹𝑋₁).
(8) If ‖𝑅(𝑋₁)‖_𝐹/‖𝑅(𝑋₀)‖_𝐹≤𝜖 stop, otherwise set 𝑋₀=𝑋₁, 𝑅(𝑋₀)=𝑅(𝑋₁) and go to step (3).

Remark 4.3. By choosing 𝑚<𝑛, Algorithm 4.2 reduces the original large sparse Sylvester-observer equation (1.2) to a low-dimensional Sylvester-observer equation (4.1). In step (5), we solve this low-dimensional matrix equation by any direct method such as the Hessenberg-Schur method. Then, in step (6) by using relation (4.5), we extend the obtained solution to the solution of the original matrix equation. Also, according to Theorem 4.1, Algorithm 4.2 is the convergence for any initial matrix 𝑋₀.

5. Weighted Block Refinement Method

In this section we discuss a new iterative method based upon a modified block refinement method. The new process uses instead of the Euclidean scalar product another one, denoted by (⋅,⋅)_𝐷 where 𝐷 is a chosen diagonal matrix. The idea of changing the inner product is to accelerate the convergence of the components of the residual which are far away from zero. To achieve this, an appropriate weight is associated to each term of the inner product. A natural choice of these weights is the entries of the first residual. The following method is based on reduction of 𝐴 and 𝐹 to the Hessenberg matrix with the use of weighted block Arnoldi process. Before giving a complete description of the new algorithm, let us define the 𝐷-scalar product.

If 𝑢 and 𝑣 are two vectors of 𝑅^𝑛, their 𝐷-scalar product is (𝑢,𝑣)_𝐷=𝑣^𝑇𝐷𝑢=^𝑛_𝑖=1𝑑_𝑖𝑢_𝑖𝑣_𝑖, (5.1) where 𝐷=diag(𝑑₁,𝑑₂,…,𝑑_𝑛) is a diagonal matrix.

This inner product is well defined if and only if the matrix 𝐷 is positive definite, that is, 𝑑_𝑖>0, for all 𝑖∈{1,…,𝑛}.

In this case, we can define the 𝐷-norm ‖⋅‖_𝐷 associated with this inner product by ‖𝑢‖_𝐷=(𝑢,𝑢)_𝐷=𝑢^𝑇D𝑢=^𝑛_𝑖=1𝑑_𝑖𝑢²_𝑖∀𝑢∈𝑅^𝑛. (5.2)

Theorem 5.1. Let 𝑋₀ be the approximate solution obtained by an arbitrary method for the matrix equation (1.2), and let 𝑅(𝑋₀)=𝐺𝐶−(𝑋₀𝐴−𝐹𝑋₀) be the corresponding residual. If 𝑚(𝑚≪𝑛) steps of the weighted block Arnoldi process by the diagonal matrices 𝐷 and 𝐷, respectively, for matrices 𝐴 and 𝐹 have been run and 𝑌_𝑚∈𝑅^𝑚×𝑚 is the solution of the low-dimensional Sylvester equation 𝑌_𝑚𝐴_𝑚−𝐹^𝑇_𝑚𝑌_𝑚=𝑉^𝑇_𝑚𝑅𝑋₀𝑈_𝑚, (5.3) then 𝑉^𝑇_𝑚𝑅𝑋₁𝑈_𝑚=𝑉^𝑇_𝑚𝑅𝑋₀+𝐷𝑉_𝑚𝑌_𝑚𝑈^𝑇_𝑚𝐷𝑈_𝑚=0. (5.4)

Proof. By using the weighted Arnoldi process, we generate the bases 𝑈_𝑚=[𝑢₁,…,𝑢_𝑚] and 𝑉_𝑚=[𝑣₁,…,𝑣_𝑚] that are, respectively, 𝐷-orthonormal and 𝐷-orthonormal; thus it holds that 𝑈^𝑇_𝑚𝐷𝑈_𝑚=𝐼_𝑚,𝑉^𝑇_𝑚𝐷𝑉_𝑚=𝐼_𝑚, (5.5) where 𝑈_𝑚,𝑉_𝑚∈𝑅^𝑛×𝑚 and 𝐷,𝐷∈𝑅^𝑛×𝑛 are two diagonal matrices.
Moreover, the square Hessenberg matrices 𝐴_𝑚 and 𝐹_𝑚 whose nonzero entries are the scalars 𝑎_𝑖𝑗 and 𝑓_𝑖𝑗, constructed by the weighted Arnoldi process, can be expressed in the form 𝐴_𝑚=𝑈^𝑇_𝑚𝐷𝐴𝑈_𝑚,𝐹_𝑚=𝑉^𝑇_𝑚𝐷𝐹𝑉_𝑚. (5.6) Now, we set 𝑋₁=𝑋₀+𝐸₀, (5.7) where 𝐸₀=𝐷𝑉_𝑚𝑌_𝑚𝑈^𝑇_𝑚𝐷 and 𝑌_𝑚∈𝑅^𝑚×𝑚 is the solution of the Sylvester-obsever equation (5.3). Thus, the new residual matrix becomes 𝑅𝑋₁=𝐺𝐶−𝑋₁𝐴−𝐹𝑋₁,=𝐺𝐶−𝑋₀+𝐸₀𝐴−𝐹𝑋₀+𝐸₀,=𝑅𝑋₀−𝐸₀𝐴−𝐹𝐸₀,=𝑅𝑋₀−𝐷𝑉_𝑚𝑌_𝑚𝑈^𝑇_𝑚𝐷𝐴−𝐹𝐷𝑉_𝑚𝑌_𝑚𝑈^𝑇_𝑚𝐷. (5.8) Multiplying the above relation on the left by 𝑉^𝑇_𝑚 and on the right by 𝑈_𝑚, we have 𝑉^𝑇_𝑚𝑅𝑋₁𝑈_𝑚=𝑉^𝑇_𝑚𝑅𝑋₀𝑈_𝑚−𝑉^𝑇_𝑚𝐷𝑉_𝑚𝑌_𝑚𝑈^𝑇_𝑚𝐷𝐴𝑈_𝑚−𝑉^𝑇_𝑚𝐹𝐷𝑉_𝑚𝑌_𝑚𝑈^𝑇_𝑚𝐷𝑈_𝑚. (5.9) Now, by using (5.3), (5.5), and (5.6) we get 𝑉^𝑇_𝑚𝑅𝑋₁𝑈_𝑚=𝑉^𝑇_𝑚𝑅𝑋₀𝑈_𝑚−𝑌_𝑚𝐴_𝑚−𝐹^𝑇_𝑚𝑌_𝑚=0. (5.10)

In order to get 𝑌_𝑚∈𝑅^𝑚×𝑚, we need to solve the low-dimensional Sylvester equation (5.3). According to the results, we can develop an iterative method for solving of the Sylvester-observer equation. The algorithm is as follows.

Algorithm 5.2 (weighted block refinement (WBR) method). (1) Start: choose an initial solution 𝑋₀, new dimension 𝑚 lesser than 𝑛 and a tolerance 𝜖.
(2) Compute 𝑅(𝑋₀)=𝐺𝐶−(𝑋₀𝐴−𝐹𝑋₀).
(3) Construct the 𝐷-orthonormal basis 𝑈_𝑚∈𝑅^𝑛×𝑚 and 𝐷-orthonormal basis 𝑉_𝑚∈𝑅^𝑛×𝑚 by the weighted Arnoldi process, such that 𝐴_𝑚=𝑈^𝑇_𝑚𝐷𝐴𝑈_𝑚,𝐹_𝑚=𝑉^𝑇_𝑚𝐷𝐹𝑉_𝑚. (5.11)
(4) Solve the low-dimensional Sylvester-observer equation 𝑌_𝑚𝐴_𝑚−𝐹^𝑇_𝑚𝑌_𝑚=𝑉^𝑇_𝑚𝑅(𝑋₀)𝑈_𝑚 with the Hessenberg-Schur method.
(5) Set 𝑋₁=𝑋₀+𝐷𝑉_𝑚𝑌_𝑚𝑈^𝑇_𝑚𝐷.
(6) Compute residual matrix 𝑅(𝑋₁)=𝐺𝐶−(𝑋₁𝐴−𝐹𝑋₁).
(7) If ‖𝑅(𝑋₁)‖_𝐹/‖𝑅(𝑋₀)‖_𝐹≤𝜖 stop, otherwise set 𝑋₀=𝑋₁, 𝑅(𝑋₀)=𝑅(𝑋₁) and go to step (2).

Remark 5.3. By choosing 𝑚<𝑛, Algorithm 5.2 reduces the original large sparse Sylvester-observer equation (1.2) to a low-dimensional Sylvester-observer equation (5.3). In step (4), we solve this low-dimensional matrix equation by any direct method such as the Hessenberg-Schur method. Also, according to Theorem 5.1, Algorithm 5.2 is the convergence for any initial matrix 𝑋₀.

6. Numerical Experiments

In this section, we present some numerical examples to illustrate the effectiveness of our new iterative methods for solving large sparse Sylvester-obsever equation. In Examples 6.1 and 6.2, we apply Algorithms 2 and 3 for solving matrix equation (1.2). In Example 6.3, we compare the Hessenberg-Schur method described in [1] with our new algorithms for solving large sparse Sylvester-obsever equation. In order to show the efficiency of our algorithms, we choose the matrices 𝐴 and 𝐶 arbitrary in these three examples. But in Example 6.4, we use four matrices from MATLAB matrix collection with the large estimation of condition numbers.

The initial approximate solution is 𝑋₀=0_𝑛×𝑛. The error is monitored by means of the test𝑅𝑋_1𝐹𝑅𝑋_0𝐹≤𝜖, (6.1) with the value of 𝜖 depending on the examples. The time is given in seconds for all examples. All numerical tests are performed in MATLAB software on a PC with 2.20 GHz with main memory 2 GB.

Example 6.1. For the first test we use two arbitrary matrices 𝐴 and 𝐶. We choose the matrices 𝐹 and 𝐺 completely satisfying the controllability requirement of the pair (𝐺,𝐶). Now We apply block refinement method for solving Sylvester-observer equation 𝑋𝐴−𝐹𝑋=𝐺𝐶 with 𝑛=200. Also, we take 𝜖=10⁻⁶. In Table 1, we report the results for different values of 𝑚. In Table 1, the results show that by increasing the values of 𝑚 and 𝑙, the number of iterations decreases. The last column of Table 1 also shows the decreasing of time consumption. Note that the fourth and fifth columns of this table are the errors of the orthogonalization method. The desired accuracy has been chosen as 10⁻⁶, but the model works well with any choice of 10^−𝑡: 𝐴=101.2.42.82.3.80……01.8101.2.42.82.3.80⋱01.61.8⋱⋱⋱⋱⋱⋱0⋮1.641.6⋱⋱⋱⋱⋱⋱.801.31.64⋱⋱⋱⋱⋱⋱2.3.81.611.3⋱⋱⋱⋱⋱⋱.82.301.61⋱⋱⋱⋱⋱⋱.42.8⋮0⋱⋱⋱⋱⋱⋱1.2.42⋮⋮01.611.31.641.61.8101.200…01.611.31.641.61.810_𝑛×𝑛𝐶=(1,0,0,…,0). (6.2)

Table 1 

Implementation of block refinement method to solve the Sylvester equation with different values of 𝑚.

Example 6.2. Consider Example 6.1 again. We apply the weighted Block refinement method for solving 𝑋𝐴−𝐹𝑋=𝐺𝐶 and take 𝜖=10⁻⁶. In Table 2, we report the results for different values of 𝑚.

Table 2 

Implementation of weighted block refinement method to solve the Sylvester-observer equation with different values of 𝑚.

Example 6.3. According to the results in Tables 1 and 2, we see that the weighted block refinement method in comparison with block refinement method works better. Now, consider that 𝐴 and 𝐶 are the same matrices used in Example 6.1. We apply our two iterative methods with 2 iterations and the Hessenberg-Schur method to solve the Sylvester-observer equation when the dimensions of the matrices are large. Results are shown in Table 3.

Table 3 

Implementation of new Iterative methods with 2 iterations and the Hessenberg-Schur method for solving the Sylvester equation.

Example 6.4. In this example we show that the convergence of our proposed algorithms independent of the matrices structure. In this sense, we use four matrices from MATLAB collection for the matrix 𝐴. The first matrix is a sparse, random finite element matrix with the condition number 1.84𝐸+03. The second matrix is a symmetric, positive semidefinite (SPD) Toeplitz matrix that is composed of the sum of 800 rank 2 SPD Toeplitz matrices with the condition number 6.18𝐸+04. The third matrix is a row diagonally dominant matrix with the condition number 1.35𝐸+010. The last matrix is a sparse singular, row diagonally dominant matrix resulting from discrediting the Neumann problem with the usual five-point operator on a regular mesh. The estimated condition number is 5.56𝐸+017. For all of these examples, the matrix 𝐶 is 𝐶=sprand(1,𝑛,𝑑), where is a random, sparse matrix with approximately 𝑑⋅𝑛⋅𝑛 uniformly distributed nonzero entries with 𝑑=0.5. We choose the matrices 𝐹 and 𝐺 completely satisfying the controllability and observability requirement of the pairs (𝐹,𝐺) and (𝐺,𝐶). We apply the Hessenberg-Schur method and weighted block refinement method (Algorithm 5.2) with 3 iterations for solving the Sylvester-observer equation 𝑋𝐴−𝐹𝑋=𝐺𝐶. The results are shown in Table 4.
It is also obvious from Table 4 that the performance of the weighted block refinement method is much better than that of the Hessenberg-Schur method, specifically for the ill-conditioned matrices.

Table 4 

Effectiveness of the Hessenberg-Schur method and WBR algorithm with 3 iterations for randomly generated matrices.

7. Comments and Conclusion

In this paper, we propose two new iterative algorithm s for solving the large sparse Sylvester-obsever matrix equations. The existing projection methods use the Arnoldi process, but the methods described in this paper are based on the weighted block Arnoldi process. Moreover, the refinement process presented in Sections 4 and 5 has the capability of improving the results obtained by an arbitrary method. The numerical examples show the efficiency of the proposed schemes.