The Robust Pole Assignment Problem for Second-Order Systems

Abstract

Pole assignment problems are special algebraic inverse eigenvalue problems. In this paper, we research numerical methods of the robust pole assignment problem for second-order systems. The problem is formulated as an optimization problem. Depending upon whether the prescribed eigenvalues are real or complex, we separate the discussion into two cases and propose two algorithms for solving this problem. Numerical examples show that the problem of the robust eigenvalue assignment for the quadratic pencil can be solved effectively.

1. Introduction

Pole assignment problems are special algebraic inverse eigenvalue problems [1, 2]. The nature model for the vibrating systems, arising in a wide range of applications, especially in the design and analysis of vibration structures, such as bridges, buildings, and airplanes, can be described by the second-order differential equation. The properties of systems of second-order differential equation are governed by its associated quadratic eigenvalue problem (QEP). To avoid unwanted oscillations of the vibratory system, pole assignment can change the poles of the system by choosing a control force and improve its stability.

Consider the following second-order matrix differential equation:𝑀𝑧(̈𝑡)+𝐶𝑧(̇𝑡)+𝐾𝑧(𝑡)=0, (1.1) where the dots denote differentiation with respect to time, and 𝑀,𝐶, and 𝐾 are 𝑛×𝑛 real symmetric matrices; 𝑀 is positive definite (denoted by 𝑀>0). Separation of variables𝑍(𝑡)=𝑥𝑒^𝜆𝑡, (1.2) where 𝜆∈𝐂, 𝑥∈𝐂^𝐧 is a constant vector, in (1.1), leads to the following quadratic eigenvalue problem:𝑄(𝜆)𝑥=(𝜆²𝑀+𝜆𝐶+𝐾)𝑥=0, (1.3) where 𝑄(𝜆)=𝜆²𝑀+𝜆𝐶+𝐾 is quadratic pencil. In general, the 2𝑛 eigenvalues of 𝑄(𝜆) are named poles of system (1.1). In engineering, the dynamics of equation (1.1) can be modified by applying a control force 𝐵𝑢(𝑡), where 𝐵∈𝑅^𝑛×𝑚,𝑢(𝑡)∈𝑅^𝑚. The relation (1.1) now becomes𝑀𝑧(̈𝑡)+𝐶𝑧(̇𝑡)+𝐾𝑧(𝑡)=𝐵𝑢(𝑡). (1.4) The special choice𝑢(𝑡)=𝐹^𝑇𝑧(̇𝑡)+𝐺^𝑇𝑧(𝑡), (1.5) where 𝐹 and 𝐺 are 𝑛×𝑚 real matrices, is called state feedback control, and (1.1) becomes𝑀𝑧(̈𝑡)+(𝐶−𝐵𝐹^𝑇)𝑧(̇𝑡)+(𝐾−𝐵𝐺^𝑇)𝑧(𝑡)=0. (1.6) The associated quadratic eigenvalue problem becomes𝑄_𝑐(𝜆)𝑥=(𝜆²𝑀+𝜆(𝐶−𝐵𝐹^𝑇)+(𝐾−𝐵𝐺^𝑇))𝑥=0, (1.7) where 𝑄_𝑐(𝜆)=𝜆²𝑀+𝜆(𝐶−𝐵𝐹^𝑇)+𝐾−𝐵𝐺^𝑇 is a closed-loop pencil.

The quadratic eigenvalue assignment problem is stated explicitly as follows.

Problem QEA
Given real matrices 𝑀,𝐶,𝐾∈𝑅^𝑛×𝑛 and 𝐵∈𝑅^𝑛×𝑚, and a set of 2𝑛 complex numbers ℒ={𝜆₁,𝜆₂,…,𝜆_2𝑛}, closed under complex conjugation, find real matrices 𝐹,𝐺∈𝑅^𝑛×𝑚, such that the eigenvalues of 𝑄_𝑐(𝜆) are equal to 𝜆_𝑗(𝑗=1,2,…,2𝑛).

Conditions for the existence of solutions to problem QEA are known as in the following theorem.

Theorem 1.1 (see [3]). Solutions 𝐹,𝐺 to problem QEA exist for every set ℒ of self-conjugate complex numbers if and only if the system (1.1) is completely controllable, that is rank[𝜆²𝑀+𝜆𝐶+𝐾,𝐵]=𝑛,∀𝜆∈𝐂. (1.8)

In a realistic situation, it is desirable to choose the feedback to ensure that the eigenstructure of closed-loop system is as robust, or insensitive to perturbation in the system matrices 𝑀,𝐶−𝐵𝐹^𝑇, and 𝐾−𝐵𝐺^𝑇, as possible to the following inverse eigenvalue problem, known as the robust quadratic eigenvalue assignment problem.

Problem RQEA
Given real matrices 𝑀,𝐶,𝐾,𝐵, and a set of ℒ as in problem QEA, find real matrices 𝐹,𝐺∈𝑅^𝑛×𝑚 and matrix 𝑋∈𝐶^𝑛×2𝑛, such that 𝑋=[𝑋^𝑇,(Λ𝑋)^𝑇]^𝑇 is nonsingular, satisfying 𝑀𝑋Λ²+(𝐶−𝐵𝐹^𝑇)𝑋Λ+(𝐾−𝐵𝐺^𝑇)𝑋=0,Λ=diag{𝜆₁,𝜆₂,…,𝜆_2𝑛}, (1.9) and the eigenstructure of the closed-loop system (1.6) is as robust as possible.

We remark that the requirement that the matrix Λ is diagonal, together with the invertibility of 𝑋. These require the prescribed eigenvalues to be distinct. In the next section, we derive conditions for the solution of problem RQEA.

In the majority of methods that have been proposed for solving problem RQEA, the second-order control system (1.1) is rewritten as a first-order system. There are two difficulties in using this approach. The first is the linear system which has a double dimension of the original quadratic system and, hence, the computational work used to solve the problem is greater than necessary. The second difficulty arises because all the exploitable properties such as definiteness, and sparsity, of the coefficient matrices 𝑀,𝐶, and 𝐾, usually offered by a practice problem, will be completely destroyed. So it is natural to wonder if solutions of the robust pole assignment problem can be obtained without resorting to a first-order reformulation.

Over the past years, many techniques for pole assignment without linearization have been proposed. Chu and Datta [4] gave a method to solve Problem RQEA. Nichols and Kautsky [5] recently have proposed a numerical method to solve this problem without linearization; the measures of robustness are subject to structured perturbations.

Datta and Sarkissian [6] proposed a direct partial modal approach to solve the partial eigenvalue assignment problem for second-order systems. It is “direct,” because the solutions are obtained directly in the second-order system without any types of reformulations; It is “partial modal,” because only a part of spectral data is needed for the solution. This method do not, however, ensure the robustness of the closed-loop system. Qian [7], Qian and Xu [8] recently have proposed a direct method to solve the robust partial pole assignment problem for second-order systems which seems more efficient and reliable. Bai et al. ([9]) gave a new optimization approach for solving this problem.

According to new measures ([10, 11]), we suggest in this paper two numerical methods for solving the problem RQEA. Numerical results show that this problem can be solved effectively.

We begin by presenting the solutions to problem RQEA without linearization. In Section 3, we describe a new measure of robustness and formulate the problem as an optimization problem. Two numerical methods are developed in Section 4, and numerical results are given in Section 5.

2. Solution to Problem RQEA

Without loss of generality, We assume in this paper that the system (1.1) is completely controllable, and 𝐵 is of full column rank. The next theorem provides necessary and sufficient conditions for the existence of solutions to problem RQEA.

Theorem 2.1 (see [5]). Let 𝑋∈𝐂^{𝐧×𝟐𝐧} be such that 𝑋=[𝑋^𝑇,(Λ𝑋)^𝑇]^𝑇 is nonsingular, where Λ=diag{𝜆₁,𝜆₂,…,𝜆_2𝑛}, then there exist real matrices 𝐹,𝐺, satisfying condition (1.9) of problem RQEA if and only if 𝑈^𝑇₁(𝑀𝑋Λ²+𝐶𝑋Λ+𝐾𝑋)=0, (2.1) where 𝐵=𝑈₀,𝑈₁𝑍0, (2.2) with 𝑈=[𝑈₀,𝑈₁] orthogonal and 𝑍 nonsingular. The matrices 𝐹,𝐺 are given by [𝐺^𝑇,𝐹^𝑇]=𝑍⁻¹𝑈^𝑇₀[𝑀𝑋Λ²+𝐶𝑋Λ+𝐾𝑋]𝑋⁻¹. (2.3)

An immediate consequence of Theorem 2.1 is the following.

Corollary 2.2. Let 𝑥_𝑗 be the right eigenvector of 𝑄_𝑐(𝜆) corresponding to the prescribed eigenvalue 𝜆_𝑗∈ℒ, then 𝑥_𝑗∈𝒩{𝑈^𝑇₁(𝜆²_𝑗𝑀+𝜆_𝑗𝐶+𝐾)}=𝒮_𝑗(𝑗=1,2,…,2𝑛), (2.4) where 𝒩{⋅} denotes right nullspace.

For every 𝜆_𝑗, find an orthogonal basis, comprised by columns of matrix 𝒮_𝑗 for the space 𝒮_𝑗,𝑗=1,2,…,2𝑛. Observe that the matrix 𝒮_𝑗 can be obtained by the QR decomposition of (𝑈^𝑇₁(𝜆²_𝑗𝑀+𝜆_𝑗𝐶+𝐾))^𝑇.

Theorem 2.3. If the system (1.1) is completely controllable, then dim(𝒮_𝑗)=𝑚,(𝑗=1,2,…,2𝑛), (2.5) where dim denotes the dimension of space 𝒮_𝑗.

Proof. From Theorem 1.1, we can get rank[𝐵,𝜆²_𝑗𝑀+𝜆_𝑗𝐶+𝐾]=𝑛, (2.6) then 𝑈𝑛=rank^𝑇0𝑈^𝑇₁[𝐵,𝜆²_𝑗𝑀+𝜆_𝑗𝐶+𝐾]=rank𝑍𝑈^𝑇₀(𝜆²_𝑗𝑀+𝜆_𝑗𝐶+𝐾)0𝑈^𝑇₁(𝜆²_𝑗𝑀+𝜆_𝑗𝐶+𝐾)=𝑚+rank[𝑈^𝑇₁(𝜆²_𝑗𝑀+𝜆_𝑗𝐶+𝐾)]. (2.7) So we can have rank[𝑈^𝑇₁(𝜆²_𝑗𝑀+𝜆_𝑗𝐶+𝐾)]=𝑛−𝑚, from which (2.5) easily follows.

According to the previous theory, we can write the following:𝑥_𝑗=𝑆_𝑗𝑤_𝑗∈𝒮_𝑗,(𝑗=1,2,…,2𝑛), (2.8) where 𝑤_𝑗∈𝐂^𝑚, and if 𝜆_𝑗=𝜆_𝑘, then 𝑥_𝑗=𝑥_𝑘. Because in the case where 𝜆_𝑗 is complex the associated eigenvector is a complex vector, in order to ensure that the computed feedback matrices are real, the eigenvector corresponding to the conjugate eigenvalue 𝜆_𝑗 must be taken to be the conjugate vector 𝑥_𝑗.

3. The Measures of Robustness

In this section, we present the measures of robustness for the second-order system. The eigenvectors of the closed-loop system can be selected to minimize the measure of robustness.

Consider a matrix 𝐴∈𝑅^2𝑛×2𝑛, if 𝐴 is nondefective, namely, then there exists a nonsingular matrix 𝑌=[𝑦₁,𝑦₂,…,𝑦_2𝑛]∈𝐂^{𝟐𝐧×𝟐𝐧}, such that𝐴𝑌=𝑌Λ, (3.1) where Λ=diag(𝜆₁,𝜆₂,…,𝜆_2𝑛), let𝑍^𝐻=[𝑧₁,𝑧₂,…,𝑧_2𝑛]^𝐻=𝑌⁻¹, (3.2) then the sensitivity of the eigenvalues of 𝜆_𝑗 of 𝐴 to perturbations in the components of 𝐴 depends upon the magnitude of the condition number 𝑐_𝑗, where𝑐_𝑗=‖𝑦_𝑗‖₂‖𝑧_𝑗‖₂≥1. (3.3) Hence, every reasonable measure of the magnitude of the vector 𝑐=(𝑐₁,𝑐₂,…,𝑐_2𝑛)^𝑇 is a reflection of the robustness of the system.

If 𝜆₁ is a multiple eigenvalue of 𝐴 and𝜆₁=𝜆₂⋯=𝜆_𝑟,𝜆₁≠𝜆_𝑗,for𝑗=𝑟+1,𝑟+2,…,2𝑛, (3.4) then the sensitivity of eigenvalue 𝜆₁ depends upon the magnitude of the number𝑐₁̃=max{𝑐₁,𝑐₂,…,𝑐_𝑟}. (3.5) Therefore, it is natural to take𝒱=(2𝑛)⁻¹‖𝑌‖_𝐹‖𝑌⁻¹‖_𝐹∶=(2𝑛)⁻¹𝜅_𝐹(𝑌) (3.6) as global measure of the sensitivity of eigenvalue.

In essence, the aim of the robust eigenvalue assignment problem is to select eigenvectors 𝑦_𝑗, such that ‖𝑦_𝑗‖₂=1 and the vectors 𝑦_𝑗 are as orthogonal as possible to each other. Therefore, in this paper we first consider the following measure 𝒱_ℎ(𝑑):𝒱_ℎ∑(𝑑)=[^{1≤𝑖<𝑗≤2𝑛}𝑑²_𝑖𝑗|𝑦^𝐻_𝑖𝑦_𝑗|²∑_{1≤𝑖<𝑗≤2𝑛}𝑑²_𝑖𝑗]^1/2, (3.7) where𝑑=[𝑑₁₂,𝑑₁₃,…,𝑑_1,2𝑛,𝑑₂₃,𝑑₂₄,…,𝑑_2,2𝑛,…,𝑑_{2𝑛−1,2𝑛}]^𝑇,𝑑_𝑖𝑗>0,∀𝑖,𝑗. (3.8)

Observe that if 𝑀 is nonsingular, then the QEP (1.7) can be formulated as the following standard eigenvalue problem:𝐴𝑥≡𝜆𝑥𝑀0𝐼₋₁(𝐾−𝐵𝐺^𝑇)𝑀⁻¹(𝐶−𝐵𝐹^𝑇)𝑥𝑥𝜆𝑥=𝜆𝜆𝑥. (3.9) Let 𝑋 be the matrix comprised by the right eigenvectors of (3.9), then𝑥𝑋=¹𝑥₂…𝑥_2𝑛𝜆₁𝑥₁𝜆₂𝑥₂…𝜆_2𝑛𝑥_2𝑛∶=[𝑥₁,𝑥₂,…,𝑥_2𝑛]. (3.10) By (2.8), we have𝑥_𝑗=𝐼𝜆_𝑗𝐼𝑥_𝑗=𝐼𝜆_𝑗𝐼𝑆_𝑗𝑤_𝑗. (3.11) From (3.7) it follows that[𝒱_ℎ(𝑑)]²=∑_{1≤𝑖<𝑗≤2𝑛}𝑑²_𝑖𝑗|𝑥^𝐻_𝑖𝑥_𝑗|²∑_{1≤𝑖<𝑗≤2𝑛}𝑑²_𝑖𝑗=_{1≤𝑖<𝑗≤2𝑛}𝛿_𝑖𝑗|𝑥^𝐻_𝑖𝑥_𝑗|²=_{1≤𝑖<𝑗≤2𝑛}𝛿_𝑖𝑗|𝑥^𝐻_𝑖[𝐼,𝜆_𝑖𝐼𝐼]𝑥𝜆𝑗𝐼^𝑗|²=_{1≤𝑖<𝑗≤2𝑛}𝛿_𝑖𝑗|𝑤^𝐻_𝑖𝑆^𝐻_𝑖[𝐼,𝜆_𝑖𝐼𝐼]𝑆𝜆𝑗𝐼^𝑗𝑤_𝑗|²≡𝑓(𝑤), (3.12) where𝑤=[𝑤^𝐻₁,𝑤^𝐻₂,…,𝑤^𝐻_2𝑛]^𝐻,𝑤_𝑗=[𝑤_1𝑗,𝑤_2𝑗,…,𝑤_𝑚𝑗]^𝐻∈𝐂^𝐦𝛿(𝑗=1,2,…,2𝑛),_𝑖𝑗=𝑑²_𝑖𝑗∑_{1≤𝑖<𝑗≤2𝑛}𝑑²_𝑖𝑗>0,∀𝑖,𝑗. (3.13)(3.14) Thus, we must solve an unconstrained optimization problemmin𝑓(𝑤), (3.15) where 𝑓(𝑤) and 𝑤 are defined by (3.12) and (3.13), 𝑤∈C^N and 𝑁=2𝑚𝑛.

4. Numerical Methods

Unfortunately, it is difficult to solve the optimization problem (3.15). Depending upon whether the prescribed eigenvalues are real or complex, we separate the discussion into two cases.

Case 1. Assume that the prescribed eigenvalues are real.
In this case, where the eigenvectors are real vectors, (3.12) becomes 𝑓(𝑤)=_{1≤𝑖<𝑗≤2𝑛}𝛿_𝑖𝑗(𝑤^𝑇_𝑖𝑆^𝑇_𝑖[𝐼,𝜆_𝑖𝐼𝐼]𝑆𝜆𝑗𝐼^𝑗𝑤_𝑗)². (4.1) Let 𝒯={𝑤=[𝑤^𝑇₁,…,𝑤^𝑇_2𝑛]^𝑇∈𝑅^𝑁∶𝑤_𝑗∈𝑅^𝑚,‖𝑤_𝑗‖₂𝐴=1∀𝑗},_𝑖𝑗=𝛾_𝑖𝑗𝑆^𝑇_𝑖[𝐼,𝜆_𝑖𝐼𝐼]𝑆𝜆𝑗𝐼^𝑗,𝑟_𝑖𝑗(𝑤)=𝑤^𝑇_𝑖𝐴_𝑖𝑗𝑤_𝑗,𝐴_𝑖𝑗(𝑤)=𝑟_𝑖𝑗(𝑤)𝐴_𝑖𝑗𝛼,∀𝑖≠𝑗,𝐴(𝑤)=²𝐼−𝐴₁₂(𝑤)…−𝐴_1,2𝑛(𝑤)−𝐴₁₂(𝑤)^𝑇𝛼²𝐼…−𝐴_2,2𝑛(𝑤)⋮⋮⋮⋮−𝐴_1,2𝑛−1(𝑤)^𝑇⋯𝛼²𝐼−𝐴_{2𝑛−1,2𝑛}(𝑤)−𝐴_1,2𝑛(𝑤)^𝑇⋯−𝐴_{2𝑛−1,2𝑛}(𝑤)^𝑇𝛼²𝐼,𝛼>0, (4.2) where 𝛾_𝑖𝑗=𝛿_𝑖𝑗, then minimizing the measure 𝒱_ℎ(𝑑) may be reduced to solving the following nonlinear programming problem with constrains: max𝑤^𝑇𝐴(𝑤)𝑤,subjectto𝑤∈𝒯. (4.3)

Letℬ={𝑤=[𝑤^𝑇₁,…,𝑤^𝑇_2𝑛]^𝑇∈𝑅^𝑁∶𝑤_𝑗∈𝑅^𝑚,‖𝑤_𝑗‖₂≤1∀𝑗}, (4.4) we can prove easily that the programming problem (4.3) is equivalent to the following programming problem with inequality constrains:max𝑤^𝑇𝐴(𝑤)𝑤,subjectto𝑤∈ℬ. (4.5)

Let𝐷={𝐷=diag(𝜏₁𝐼^(𝑚),…,𝜏_2𝑛𝐼^(𝑚))∶𝜏_𝑖≥0,∀𝑖}. (4.6) We consider a multiparameter eigenvalue problem𝐴(𝑤)𝑤=𝐷𝑤,𝐷∈𝐷,𝑤∈𝒯. (4.7) Using the Kuhn-Tucker optimality ([12, 13]), we can prove that every solution to problem (4.5) is necessarily a solution to problem (4.7).

In this section, we develop an algorithm, called Horst algorithm ([14]), to solve (4.7), for finding a global optimal solution 𝑤 of problem (4.3).

Algorithm 1. Initialization Step
Let 𝜀>0 be the termination scalar. Given 𝑀,𝐶,𝐾∈𝑅^𝑛×𝑛,𝐵∈𝑅^𝑛×𝑚,Λ=diag(𝜆₁,𝜆₂,…,𝜆_2𝑛)∈𝐂^2𝐧×2𝐧, choose an initial vector 𝑤⁽⁰⁾=[𝑤₁^𝑇(0),𝑤₂^𝑇(0),…,𝑤^𝑇(0)_2𝑛]^𝑇∈𝒯^𝑁.
Main Step
   (1) Find the decomposition (2.2) of 𝐵 and an orthonormal basis, comprised by the columns of the matrix 𝑆_𝑗, for the subspaces 𝒮_𝑗,𝑗=1,2,…,2𝑛, defined by (2.4). Let 𝑘=0.
(2) For 𝑖=1,2,…,2𝑛, we compute 𝑧_𝑖^(𝑘+1)=𝛼²𝑤_𝑖^(𝑘)−^𝑖−1_𝑗=1𝑤_𝑖^𝑇(𝑘)𝐴_𝑖𝑗𝑤_𝑗^(𝑘+1)𝐴_𝑖𝑗𝑤_𝑗^(𝑘+1)−^2𝑛_𝑗=𝑖+1𝑤_𝑖^𝑇(𝑘)𝐴_𝑖𝑗𝑤_𝑗^(𝑘)𝐴_𝑖𝑗𝑤_𝑗^(𝑘),𝜏_𝑖^(𝑘)=‖𝑧_𝑖^(𝑘+1)‖₂,𝑤_𝑖^(𝑘+1)=𝑧_𝑖^(𝑘+1)𝜏_𝑖^(𝑘)if𝜏_𝑖^(𝑘)𝑤>0,_𝑖(𝑘)if𝜏_𝑖^(𝑘)=0. (4.8)
(3) Set Δ𝑤^(𝑘)=𝑤^(𝑘+1)−𝑤^(𝑘). If ‖Δ𝑤^(𝑘)‖₂≤𝜀, then 𝑤^∗=𝑤^(𝑘+1) is an approximation optimal solution, go to (4); if ‖Δ𝑤^(𝑘)‖₂>𝜀, then replace 𝑘 by 𝑘+1 and repeat step (2). 
  (4) Let 𝑋=[𝑥₁,𝑥₂,…,𝑥_2𝑛], where 𝑥_𝑗 is defined by (2.8), and construct feedback matrices 𝐹,𝐺 by solving [𝐺^𝑇,𝐹^𝑇]=𝑍⁻¹𝑈^𝑇₀[𝑀𝑋Λ²+𝐶𝑋Λ+𝐾𝑋]𝑋⁻¹. (4.9)

Case 2. Assume that prescribed eigenvalues are complex.
In this case, it is almost impossible to solve optimization problem (3.15) without costing too much. Based on this ideal, we consider 𝜅_𝐹(𝑋)=‖𝑋‖_𝐹‖𝑋⁻¹‖_𝐹, where 𝑋 is defined by (3.10), and use some similar techniques developed by Kautsky et al. in 1985 ([15]).
Obviously, 𝜅_𝐹(𝑋)=‖𝑋‖_𝐹‖𝑋⁻¹‖_𝐹 achieves the minimum if and only if 𝑋 is unitary.

Let𝑊_𝑗𝑈=𝑁(^𝑇1(𝜆²_𝑗𝑀+𝜆_𝑗𝜆𝐶+𝐾)_𝑗𝑈^𝑇₁(𝜆²_𝑗𝑀+𝜆_𝑗𝐶+𝐾)), (4.10) then 𝑥_𝑗∈𝑊_𝑗. Let the columns of matrix 𝑊_𝑗 comprise an orthonormal basis, for the subspace 𝑊_𝑗, and𝑋_𝑗=span{𝑥₁,…,𝑥_𝑗−1,𝑥_𝑗+1,…,𝑥_2𝑛}, (4.11) and let 𝑦_𝑗 be a normalized vector orthogonal to 𝑋_𝑗. The objective here is to choose vectors 𝑥_𝑗(𝑗=1,2,…,2𝑛), such that each vector is as orthogonal as possible to the space 𝑋_𝑗. Observe that 𝑥_𝑗 orthogonal to 𝑋_𝑗 is equivalent to choose 𝑥_𝑗∈𝑊_𝑗, such that the angle between 𝑥_𝑗 and 𝑦_𝑗 is minimized.

Obviously, 𝑦_𝑗 must satisfy[𝑥₁,…,𝑥_𝑗−1,𝑥_𝑗+1,…,𝑥_2𝑛]^𝐻𝑦_𝑗=0. (4.12) The coefficient matrix of equation (4.12) is a (2𝑛−1)×2𝑛 matrix, so we can easily solve 𝑦_𝑗 from equation (4.12), and let 𝑥_𝑗 be the normalized projection of 𝑦_𝑗 onto the space 𝑊_𝑗, then𝑥_𝑗=𝑊_𝑗𝑊^𝐻_𝑗𝑦_𝑗‖𝑊^𝐻_𝑗𝑦_𝑗‖₂. (4.13) This method is then to sweep through the columns of 𝑋 replacing the 𝑗th column in turn with the normalized projection 𝑥_𝑗. Because the prescribed eigenvalue 𝜆_𝑗 is complex, two columns need to be altered simultaneously. Repeat the previous procedure until satisfying the stopping criterion or reaching the maximal amount of iteration. We can choose‖𝑋^𝐻𝑋−𝐼‖_𝐹<𝜀, (4.14) as the stopping criterion.

It is worthwhile to point out that the procedure is not guaranteed to converge, that is, the criterion (4.14) may not be satisfied. To ensure the end of iteration, we need to set a maximal amount of iteration 𝑘_max. Now we develop an algorithm, called orthogonal vector method, for solving the problem RQEA.

Algorithm 2.
Initialization Step
Let 𝜀>0 be the termination scalar, 𝑘_max is the maximal amount of iteration. Given 𝑀,𝐶,𝐾∈𝑅^𝑛×𝑛,𝐵∈𝑅^𝑛×𝑚,and Λ=diag(𝜆₁,𝜆₂,…,𝜆_2𝑛)∈𝐶^2𝑛×2𝑛.
Main Step
(1) Find the decomposition (2.2) of 𝐵 and an orthonormal basis, comprised by the columns of the matrix 𝑊_𝑗, for the subspaces 𝑊_𝑗,𝑗=1,2,…,2𝑛, defined by (4.10).
(2) Select an initial matrix 𝑋=[𝑥₁,𝑥₂,…,𝑥_2𝑛], defined by (3.10), such that 𝑥_𝑗∈𝑊 and 𝑋 is nonsingular. Let 𝑘=0.
(3) If 𝑘<𝑘_max, do for 𝑗=1,2,…,2𝑛,  compute the solution 𝑦_𝑗 of equation (4.12) and normalize 𝑦_𝑗;  compute 𝑥_𝑗=𝑊_𝑗𝑊^𝐻_𝑗𝑦_𝑗/‖𝑊^𝐻_𝑗𝑦_𝑗‖₂;  compute res₁=‖𝑋^𝐻𝑋−𝐼‖_𝐹;  update the 𝑗th column of 𝑋 with 𝑥_𝑗; compute res₂=‖𝑋^𝐻𝑋−𝐼‖_𝐹.  If |res₂−res₁|<𝜀, then go to (4), else continue;
(4) let 𝑘=𝑘+1; repeat (3) until convergence.
(5) let the first 𝑛 rows of 𝑋 be 𝑋 construct feedback matrices 𝐹,𝐺 by solving [𝐺^𝑇,𝐹^𝑇]=𝑍⁻¹𝑈^𝑇₀[𝑀𝑋Λ²+𝐶𝑋Λ+𝐾𝑋]𝑋⁻¹. (4.15)

Although this method does not guarantee convergence, it is simple to implement, and numerical results show that it often leads to better conditioned closed-loop systems.

5. Numerical Results

To illustrate the performance of the present algorithms, in this section we give some numerical examples, which were carried out using MATLAB 6.1.

Example 5.1. In the first example we examine a case from [5], where the matrix 𝑀 is very ill conditioned. The system matrices 𝑀,𝐶,𝐾, and 𝐵 are defined by 𝑀=500000011011.00001,𝐶=0,𝐾=−4040040−8040040−40,𝐵=123234. (5.1) The system is undamped and the open-loop eigenvalues are {±48999.0𝑖,±4.4726𝑖,±0.051635𝑖}. (5.2) The desired closed-loop eigenvalues 𝜆_𝑗,𝑗=1,2,…,6 are given by ℒ={−1,−2,−3,−4,−5,−6}. (5.3)

With Algorithm 1, we choose an initial vector 𝑤⁽⁰⁾=[𝑤₁^𝑇(0),𝑤₂^𝑇(0),…,𝑤^𝑇(0)_2𝑛]^𝑇 with 𝑤₁⁽⁰⁾=[1,0]^𝑇,𝑤₂⁽⁰⁾=[0,1]^𝑇,𝑤₃⁽⁰⁾=1/√2[1,−1]^𝑇,𝑤₄⁽⁰⁾=1/√2[−1,1]^𝑇,𝑤₅⁽⁰⁾=[1,0]^𝑇, and 𝑤₆⁽⁰⁾=[0,1]^𝑇, and take 𝛾_𝑖𝑗=1/√6(𝑖,𝑗=1,2,…,6),𝛼²=0.7. The corresponding matrix 𝑋 has condition number 𝜅_𝐹(𝑋)=1.2626×10⁵, and after three iterations, it is reduced to 𝜅_𝐹(𝑋)=3.8416×10⁴. The computed feedback matrices are given by𝐹^𝑇=,𝐺−312.9824−7.3540−6.6748−0.0233−0.0003−0.0003_𝑇=.−194.945438.2092−79.276119.9877−60.001659.9985 (5.4)

Example 5.2. In the second example, we examine a case from [16]. The system matrices 𝑀,𝐶,𝐾, and 𝐵 are defined by 𝑀=𝐼₄,𝐶=0.5000000000000000.5,𝐾=5−500−510−500−510−500−56,𝐵=10010000. (5.5) The open-loop eigenvalues are {3.525,−3.559,−0.059±3.732𝑖,−0.191±1.489𝑖,−0.233±2.692𝑖}. (5.6) The desired closed-loop eigenvalues 𝜆_𝑗,𝑗=1,2,…,8 are given by ℒ={−1±𝑖,−2±𝑖,−3±𝑖,−4±𝑖}. (5.7)

With Algorithm 2, we choose 𝜀=10⁻⁶,𝑘_max=8, and the initial matrix 𝑋 is generated by a random selection of vectors from each subspace. The corresponding matrix 𝑋 has condition number 𝜅_𝐹(𝑋)=5.8818×10³, and after eight iterations, it is reduced to 𝜅_𝐹(𝑋)=2.3927×10³. The computed feedback matrices are given by𝐹^𝑇=,𝐺−5.5706−15.6774−92.744957.97600.0050−13.42942.44490.0238_𝑇=.−5.3794−173.3879235.9619−93.4799−5.0455−46.2669113.8836−70.6403 (5.8)

Acknowledgment

This research was supported by the National Natural Science Foundation of China under Grant no. 10926130.