Linear quadratic Nash game-based tracker for multiparameter singularly perturbed sampled-data systems: digital redesign approach

國 立 成 功 大 學 電 機 工 程 學 系 碩 士 論 文 適用於多變數殊異擾動系統的線性二次 納許賽局追蹤器：數位再設計法則 Linear Quadratic Nash Game Based Tracker for Multiparameter Singularly Perturbed Sampled-data Systems ：Digital Redesign Approach

研 究 生：楊子儀 指導教授：蔡聖鴻 中華民國 九十五 年 六 月

適用於多變數殊異擾動系統的線性二次 納許賽局追蹤器：數位再設計法則 楊子儀* 蔡聖鴻** 國立成功大學電機工程學系

摘要 在本論文中，新建立了一個適用於多變殊異擾動系統的線性二次納許賽 局追蹤器。加入 LQR 設計方法應用於追蹤器設計的兩個二次估測方程的 常態交錯多變數立卡迪函數(GCMARE)需要被解出。首先，新提出了針對 GCMARE 的漸進擴展式，且提出基於牛頓法的演算法來解 GCMARE 且保證其 收斂性。然後，藉由數位再設計的方法設計出一個低增益且有高效能的 控制器。最後，為了增進追蹤器的效能，使用了渾沌進化演算法(CEPA) 來調整追蹤器的參數。最後以一個例子來驗證提出的方法是有效的。

*

作者

** 指導教授

I

Linear Quadratic Nash Game Based Tracker for Multiparameter Singularly Perturbed Sampled-data Systems: Digital Redesign Approach Zi-Yi Yang* Jason Sheng-Hong Tsai** Department of Electrical Engineering National Cheng Kung University, Tainan, Taiwan, R.O.C.

Abstract In this thesis, a tracker for the linear quadratic Nash game of multiparameter singularly perturbed sampled-data systems is newly established. A generalized cross-coupled multiparameter algebraic Riccati equation (GCMARE) for two quadratic cost functions is needed to be solved by applying the LQR design methodology for the tracker design. Firstly, the asymptotic expansions for the GCMARE are newly established, and the proposed algorithm based on the Newton’s method for solving the GCMARE guarantees the quadratic convergence. Then the low-gain sample-data controller with a high design performance is realized though the digital redesign method. Finally, for further improving the tracking performance, the chaos evolutionary programming algorithm (CEPA) is utilized to tune the parameters of the tracker. An example is presented to demonstrate the effectiveness on the proposed methodology.

*

The author

** The advisor

II

Acknowledgments

I wish to express my sincere gratitude to my advisor, Professor Jason S. H. Tsai, for his intelligent advice and valuable guidance. I also wish to thank Professors Yu-Pin Chang and Chia-Wei Chen for serving on my thesis committee. I wish to express the most the most sincere appreciation to my family and all my friends and members in the multivariable control Lab., for their supports and encouragement during the entire course of this thesis. Without their continual help, I would not be able to complete this thesis.

III

List of Contents Page

Chinese Abstract …….................................................................................... I Abstract ……...…………..……………………………………….. II Acknowledgments ………………...………………………………………. III List of Contents …………………..………………………………………. IV List of Tables ……………………...………………………………………. V List of Figures …………………...……………………………………….. VI

Content Chapter 1 Introduction ………...….……….…………………………………… 1-1 Chapter 2 Linear Quadratic Nash Game Based Tracker for Multiparameter Singularly Perturbed Sampled-data Systems： Digital Redesign Approach……………………………………………………………...2-1 2.1 The design of controller with digital redesign method …...………….………… 2-2 2.1.1

Digital redesign approach……………………………………………..……. 2-2

2.1.2

The matrix sign function…………………………………………………….. 2-9

2.1.3

Asymptotic structure .................................................................................... 2-11

2.1.4

Newton’s method …………………………………….………….………….. 2-17

2.2 The chaos-evolutionary-programming approach ………………………… 2-23 2.2.1

The chaos system selected for optimization …………………...……..…... 2-23

2.2.2

The chaos-evolutionary-programming approach ………...……………….. 2-23

2.3 An illustrative example……………………………………………...………… 2-29

Chapter 3 Conclusion………………….…...……………………..………………. 3-1 References ………………….…….……………………………...……………….….. R-1 Appendix ………………………………………………………………………….A-1

IV

List of Tables Page Table 2.1 Natural numbers in binary scale ..................................................................... 2-24 Table 2.2 Quasi-random sequences …………………………...………………………. 2-25

V

List of Figures Page Fig. 2.1 Continuous-time control system ……………….…………………...………. 2-3 Fig. 2.2 Digitally sampled-data control system …………………..…………….……. 2-4 Fig. 2.3 Time responses of ri (t ) , yci (t ) and ydi (t ) by proposed controllers without CEPA: (a) continuous-time system (b) sampled-data system ………….….. 2-32 Fig. 2.4 The time series of uci (t ) and udi (t ) without CEPA ..……………………… 2-33 Fig. 2.5 Time responses of ri (t ) , yci (t ) and ydi (t ) by proposed controllers with CEPA: (a) continuous-time system (b) sampled-data system ………….…… 2-34 Fig. 2.6 The time series of uci (t ) and udi (t ) with CEPA …………………………... 2-34

VI

Chapter 1 Introduction

Many feedback control problems of singularly perturbed systems have been investigated extensively (see e.g.,[1], [2], [9]). In order to obtain the controller, we must solve the algebraic Riccati equation (ARE) with small positive constant parameter ε > 0 , note that it is very difficult to solve directly the singularly perturbed ARE due to high dimension and numerical stiffness ([1], [10]). The linear quadratic Nash games and their applications have been studied intensively in many papers (see e.g., [5], [6]). Starr and Ho [5] derived the closed-loop perfect-state linear Nash equilibrium strategies for a class of analytic differential games. And then the linear quadratic Nash games for the multiparameter singularly perturbed systems (MSPS) and the singularly perturbed systems (SPS) have been studied by using composite controller design [2] , [7], [8]. The recursive algorithm for the solution of ARE of singularly perturbed systems have been developed in many literatures (see e.g. [13] ). From a practical point of view, it has been shown that the recursive algorithm is very effective to solve the ARE when the

1-1

system matrices are functions of a small perturbation parameter ε . However, the recursive algorithm converge only to the approximation solution. Moreover, such an algorithm is the linear convergence. In recent years, the numerical algorithm based on Newton’s method for solving the CARE for the SPS [11], and the cross-coupled Multiparameter algebraic Riccati equations (CMARE) [12] has been proposed. However there are some programs needed to be improved. Firstly, the assumptions of the fast state matrices are needed be nonsingular. Secondly, the convergence speed is very slow. In [14], the solution to the GCMARE have been newly proved without the nonsingularity assumption. It also improves the convergence rate with high speed. In this thesis, we develop an effective digital tracker for linear quadratic Nash games of multiparameter singularly perturbed systems. The main idea to achieve these is to design an high-gain analog state-feedback controller by the algorithm based on Newton’s method for solving the GCMARE. Then the tracker is designed based on some advanced digital redesign techniques equipped with a predictive feature, which yields an equivalent digital controller but with a low gain, without losing the high quality performance as much as possible. Finally we use the chaos evolutionary programming approach (CEPA) [22] to optimal the performance and show that is good and effectively.

Notation: The notations used in this thesis are fairly standard. The superscript T denotes matrix transpose. I n denotes the n × n identity matrix.

i denotes its Euclidean norm

for a matrix. det M denotes the determinant of M . vec M denotes an ordered stack of the columns of M [4]. Let M be an m × n matrix and m j its j -th column, then vec M is the mn×1 vector

1-2

⎡ m1 ⎤ ⎢m ⎥ vec M = ⎢ 2 ⎥ . ⎢ ⎥ ⎢ ⎥ ⎣ mn ⎦

⊗ denotes Kronecker product. U lm denotes a permutation matrix in Kronecker matrix sense [4] such that U lm vecM = vecM T , M ∈ ℜl ×m . Re(λ ) denotes a real part of λ ∈ C .

1-3

Chapter 2 Linear Quadratic Nash Game Based Tracker for Multiparameter Singularly Perturbed Sampled-data Systems: Digital Redesign Approach

In this chapter, a tracker for the linear quadratic Nash game of multiparameter singularly perturbed sampled-data systems is newly established. In Sec.2.1, design the controller with digital redesign. First, we use digital redesign method to design the controller in Sec.2.1.1. Then, we use matrix sign function to turn the parameter in Sec.2.1.2. Finally, In Sec.2.2, use the chaos evolutionary programming approach to optimal the tracker for the linear quadratic Nash game of multiparameter singularly perturbed sampled-data systems. In Sec.2.3, An illustrative example is presented to demonstrate the effectively on the proposed methodology.

2-1

2.1 The design of controller with digital redesign method Here we design the controller with digital. In Sec.2.1.1 we introduction the digital redesign method. In Sec. 2.1.2, The matrix sign function is used to tune the forward gain Ec . In order to solve generalized cross-coupled multiparameter algebraic Riccati equation, use the asymptotic structure and Newton’s method in Sec 2.1.3 and Sec. 2.1.4.

2.1.1 Digital redesign approach The design of the controller used here is based on the optimal tracker [15] approach, so we can use digital redesign approach [16] to develop a sample-data controller. Consider a linear controllable continuous-time system described by xc (t ) = Axc (t ) + Buc (t ),

xc (0) = x0 ,

yc = Cxc (t ) ,

(2.1a) (2.1b)

where xc (t ) ∈ ℜ n is the state vector, uc (t ) ∈ ℜ m is the control input, yc (t ) ∈ ℜ p is the measurement output and ( A, B, C ) are system matrices of appropriate dimensions. According to the preface, we desire the closed-loop control system to track a reference trajectory r(t ) , so we can use the optimal state-feedback control law to minimize the following performance index: ∞

J = ∫ {[Cxc (t ) − r (t )]T Qc [Cxc (t ) − r (t )] + uc T (t ) Rc uc (t )}dt ,

(2.2)

0

where Qc = QcT ≥ 0 , Rc = RcT > 0 and r (t ) is the given reference input. This optimal controller is given by uc (t ) = − K c xc (t ) + Ec r (t ), where K c = Rc−1B T P ,

2-2

(2.3)

Ec = − Rc−1B T [( A − BK c ) −1 ]T C T Qc ,

and P is the solution of the following Riccati Equation AT P + PA − PBRc−1 B T P + C T QcC = 0 .

(2.4)

K c ∈ ℜ m×n , Ec ∈ ℜ m×m have been designed to satisfy some specified goals, and

r (t ) ∈ℜm is a reference input vector. . Then the controlled system is xc (t ) = Ac xc (t ) + BEc r (t ),

(2.5)

where Ac = A − BK c . The block diagram of the optimal controlled system is shown in Fig.2.1.

r (t )

+

u c (t )

Ec −

x c (t )

xc (t ) = Axc (t ) + Bu(t )

y c (t ) C

Kc Fig. 2.1. Continuous-time control system

Let the state equation of a corresponding discrete-time equivalent model be xd (t ) = Axd (t ) + Bud (t ),

xd (0) = x0 ,

yd (t ) = Cxd (t ),

(2.6a) (2.6b)

where ud (t ) ∈ ℜ m is a piecewise constant input vector, satisfying ud (t ) = ud ( kT ) ,

for

kT ≤ t < ( k + 1)T ,

and T > 0 is the sampling period. Let the discrete-time state-feedback controller be ud ( kT ) = − K d xd ( kT ) + Ed r ∗ ( kT ) ,

(2.7)

where K d ∈ℜm×n and Ed ∈ℜm×m are respective digital state-feedback gain and digital

2-3

forward gain, which are to be determined from the analog gain, K c and Ec , in (2.3). We set r * (kT ) that is a sequence of reference input vector, which equals to r (kT + T ) in order to track the original reference input r (kT ) . r * ( kT ) Ed

u d (t )

+

Z.O.H.

−

y d (t )

x d (t )

x d (t ) = Axd (t ) + Bud (t )

C

T T

x d (kT ) Kd

Fig. 2.2. Digitally sampled-data control system

xd (t ) = Axd (t ) + B[ − K d xd ( kT ) + Ed r ∗ ( kT )] , xd (0) = x0 ,

yd (t ) = Cxd (t ),

(2.8)

for kT ≤ t < ( k + 1)T .

A zero-order-hold device is used for (2.7). The digital redesign problem is to find digital controller gains ( K d , Ed ) in (2.5) from the analog gains ( K c , Ec ) in (2.3), so that the closed-loop state xd (t ) in (2.4) can closely match the closed-loop state in (2.5) at all the sampling instants for a given r (t ) = r ( kT ) , k = 0,1, 2

.

The solution xc (t ) of (2.1a) at t = tv = kT + vT for 0 ≤ v ≤ 1 is found to be tv

xc (t ) = exp( A(tv − kT )) xc (kT ) + ∫ exp( A(tv − τ )) Buc (τ )dτ . kT

(2.9)

Let uc (tv ) be a piecewise-constant input. Then (2.9) reduces to xc (tv ) ≈ exp( AvT ) xc ( kT ) + ∫

kT + vT

kT

exp( A( kT + vT − τ )) Bdτ uc (tv )

= G xc ( kT ) + H uc (tv ), (v)

(v)

where G ( v ) = exp( A(tv − kT )) = exp( AvT ) = (exp( AT )) v = (G ) v ,

2-4

(2.10)

tv

vT

kT

0

H ( v ) = ∫ exp( A(tv − τ )) Bdτ = ∫

exp( Aτ ) Bdτ

⎧ (vT ) 2 (vT )3 + A2 + ...]B, ⎪[ I n vT + A =⎨ 2! 3! ⎪(G ( v ) − I ) A−1 B, n ⎩

for singular A

.

for nonsingular A

Also, the solution xd (t ) of (2.6a) at t = tv = kT + vT for 0 ≤ v ≤ 1 is obtained as tv

xd (tv ) = exp( A(tv − kT )) xd ( kT ) + ∫ exp( A(tv − τ )) Bdτ ud ( kT )

(2.11)

kT

= G xd (kT ) + H ud ( kT ). (v)

(v)

Thus, from (2.10) and (2.11) if follows that to obtain the predicted state xc (tv ) = xd (tv ) under

the

assumption

of xc ( kT ) = xd ( kT ) ,

it

is

necessary

to

have ud ( kT ) = uc (tv ) . This leads to the prediction-based digital controller ud ( kT ) = uc (tv ) = − K c xc (tv ) + Ec r (tv ) = − K c xd (tv ) + Ec r (tv ) ,

(2.12)

where the future state xd (tv ) needs to be predicted based on the available causal signals xd ( kT ) and ud (kT ) . Substituting the predicted state xd (tv ) in (2.11) to (2.12) and then solving it for ud (kT ) yields ud ( kT ) = ( I m + K c H ( v ) ) −1[ − K cG ( v ) xd ( kT ) + Ec r (tv )] .

(2.13)

As a result, the desired predicted digital controller (2.7) is found from (2.13) to be

ud (kT ) = − K d( v ) xd (kT ) + E d( v ) r ∗ (kT ) .

(2.14)

where for tracking purpose, r ∗ (kT ) = r (kT + T ) , and K d( v ) = ( I m + K c H ( v ) ) −1 K cG ( v ) , Ed( v ) = ( I m + K c H ( v ) ) −1 Ec .

In particular, if v = 1 then the pre-requisite xc (kT ) = xd (kT ) is ensured, so that for any k = 0,1, 2

, the new predicted controller is given by

ud (kT ) = − K d xd (kT ) + E d r ∗ (kT )

2-5

(2.15)

where K d = ( I m + K c H ) −1 K c G ,

(2.16a)

Ed = ( I m + K c H ) −1 Ec ,

(2.16b)

r * (kT ) = r (kT + T ) ,

(2.16c)

G = exp( AT ) ,

(2.17a)

in which

3 ⎧ T2 2 T [ + + + ...]B, I T A A ⎪ H =⎨ n 2! 3! ⎪ (G − I ) A−1B, ⎩ n

for singular A

.

(2.17b)

for nonsingular A

Now, consider a linear quadratic Nash game of Multiparameter singular perturbed system (MSPS) as follows 2

2

k =0

k =1

x0 (t ) = ∑ A0i xi (t ) + ∑ B0i ui (t ) , x0 (t ) = x00 ,

ε i xi (t ) = Ai 0 x0 (t ) + Aii xi (t ) + Bii ui (t ) ,

(2.18)

xi (t ) = xi0 ,

yi (t ) = Ci 0 x0 (t ) + Cii xi (t ) = Ci x(t ) ,

⎡ x0 (t ) ⎤ x(t ) = ⎢⎢ x1 (t ) ⎥⎥ , ⎢⎣ x2 (t ) ⎥⎦

i =1, 2 .

The optimal state-feedback control is to minimize the following performance index [15]

Ji =

1 2

∫ {[ y ∞

0

( t ) − ri ( t ) ] Q i [ y i ( t ) − ri ( t ) ] T

i

+ u ( t ) R ii u i ( t ) } d t ,

where x i ∈ ℜ

ni

, i = 0 , 1, 2

(2.19)

i = 1, 2

T i

are the state vectors and

ui ∈ ℜ mi , i = 1, 2 are the

control inputs. All the matrices are constant matrices of appropriate dimensions.

2-6

ε1 and ε 2 are two small positive singular perturbation parameters of the same order of magnitude such that 0 < k1 ≤ α ≡

ε1 ≤ k2 < ∞ . ε2

(2.20)

It is assumed that the limit of α exists as ε1 and ε 2 tend to zero (see, e.g.,[5],[6] ),that is

α = lim α .

(2.21)

ε1 → +0 ε2 → +0

It is worth pointing out that the matrices Aii , i = 1, 2 may be singular. In fact such systems arise in some real physical applications like a flexible space structure [11]. In this case, it should be noted that the composite design [5]–[8] cannot be applied. Let us introduce the partitioned matrices

Ae := Λ e− 1 A , B ie := Λ e− 1 B i , S ie := B ie R ii− 1 B ieT = Λ e− 1 S i Λ e− 1 , i = 1, 2 0 0 ⎤ ⎡ In0 ⎢ 0 ⎥⎥ , Λ e := ⎢ 0 ε1 I n1 ⎢⎣ 0 0 ε 2 I n 2 ⎥⎦ ⎡ B01 ⎤ ⎡ B02 ⎤ ⎢ ⎥ B1 := ⎢ B11 ⎥ , B2 := ⎢⎢ 0 ⎥⎥ , ⎢⎣ 0 ⎥⎦ ⎢⎣ B22 ⎥⎦

A01 A11 0

C11

C2 = [C20

0 C 22 ]

S 011 S111 0

0⎤ 0 ⎥⎥ , 0 ⎥⎦

⎡ Q 001 := C Q 1 C 1 = ⎢⎢ Q 0T0 1 ⎢⎣ 0

Q 011 Q1 1 1 0

0⎤ 0 ⎥⎥ , 0 ⎥⎦

T 1

−1 11

B

T 1

A02 ⎤ 0 ⎥⎥ , A22 ⎥⎦

C1 = [C10

⎡ S 001 = ⎢⎢ S 0T0 1 ⎢⎣ 0

S 1 := B 1 R

Q1 n

⎡ A00 A := ⎢⎢ A10 ⎢⎣ A20

0]

,

⎡ S 002 = ⎢⎢ 0 ⎢⎣ S 0T2 2

0 0 0

S 022 ⎤ 0 ⎥⎥ , S 2 2 2 ⎥⎦

⎡ Q 002 := C Q 2 C 2 = ⎢⎢ 0 ⎢⎣ Q 0T2 2

0 0 0

Q 022 ⎤ 0 ⎥⎥ . Q 2 2 2 ⎥⎦

S 2 := B 2 R

Q2n

−1 22

B

T 2

T 2

We now consider the linear quadratic Nash games for infinite horizon nonstandard MSPS (2.18) under the following basic assumptions (see, e.g., [3,5]).

2-7

Assumption 1: There exists a µ * > 0 such that the triplet ( Ae , Bie , Ci ), i = 1, 2 are stabilizable and detectable for all µ ∈ (0, µ * ], where µ := ε1ε 2 .

Assumption 2: The triplet ( Aii , Bii , Cii ), i = 1, 2 are stabilizable and detectable.

These conditions are quite natural since at least one control agent has to be able to control and observe unstable modes. Our purpose is to find a linear feedback controller (u1* , u2* ) such that

J i ( u i* , u *j ) ≤ J i ( u i , u *j ),

i , j = 1, 2, i ≠ j .

(2.22)

The Nash inequality shows that ui* regulates the state to zero with minimum output energy. The following lemma is already known [8].

Lemma 1: Under Assumption 1, there exists an admissible controller such that the inequality (2.22) holds iff the following full-order CMARE AeT X e + X e Ae + Q1n − X e S1e X e − X e S 2 eYe − Ye S 2 e X e = 0 ,

(2.23a)

AeT Ye + Ye Ae + Q2 n − Ye S 2 eYe − Ye S1e X e − X e S1eYe = 0 .

(2.23b)

Have stabilizing solutions X e ≥ 0 and Ye ≥ 0 , where ⎡ X 00 ⎢ X e = ⎢ ε1 X 10 ⎢ ⎣⎢ε 2 X 20

ε1 X 10T ε1 X 11 ε1ε 2 X 21

T ⎤ ε1 X 20

⎡ Y00 ⎢ ⎥ T ε1ε 2 X 21 ⎥ , Ye = ⎢ ε 1Y10 ⎢ ⎥ ε 2 X 22 ⎥⎦ ⎣⎢ ε 2Y20

ε 1Y10T ε 1Y11 ε 1ε 2 Y21

ε 1Y20T

⎤ ⎥ ε 1ε 2 Y21T ⎥ . ⎥ ε 2Y22 ⎦⎥

This optimal control is given by u (t ) = − K c x ( t ) + E c r (t ) ,

(2.24)

where the analog feedback gain K c ∈ R m × n , the forward gain E c ∈ R m × m ,

2-8

K c 1 = R 1−11 B 1Te X e , K c 2 = R 2−21 B 2T e Y e ,

(2.25a)

⎡K ⎤ K c = ⎢ c1 ⎥ , ⎣ K c2 ⎦ Ec1 = − R11−1 B1Te [( Ae − B1e K c1 ) −1 ]T C1T Q1 , Ec 2 = − R22−1 B2Te [( Ae − B2 e K c 2 ) −1 ]T C2T Q2 , ⎡E Ec = ⎢ c1 ⎣ 0

(2.25b)

0 ⎤ . Ec 2 ⎥⎦

2.1.2 The matrix sign function

Unfortunately, it’s easy to find that det( Ae − B1e K c1 ) = 0 and det( Ae − B2 e K c 2 ) = 0 , which implies the inverses of ( Ae − B1e K c1 ) and ( Ae − B2e K c 2 ) do not exist. Furthermore, since rank [ Ae − λi I n

B1e ] < n for λi = 0 and rank [ Ae − λi I n

B2e ] < n for λi = 0 , which

interprets the system is uncontrollable with respect to λi = 0 . Whenever the matrix ( Ae − B1e K c1 ) and ( Ae − B2e K c 2 ) in (25b) are singular matrices, the exact evaluations of Ec1 and Ec 2 in (25b) become difficult. To find the acceptable Ec1 and Ec 2 with the above-mentioned singular matrices ( Ae − B1e K c1 ) and ( Ae − B2e K c 2 ) , we propose a matrix sign function method [23] (see Appendix) as follows. First, we carry out the bilinear transformation with a circle of a radius ζ 1 from the origin of the coordinates to formulate matrices Ab1 and Ab 2 as Ab1 = ( Ae − B1e K c1 − ζ 1 I )( Ae − B1e K c1 + ζ 1 I ) −1 ,

(2.26a)

Ab 2 = ( Ae − B2 e K c 2 − ζ 1 I )( Ae − B2 e K c 2 + ζ 1 I ) −1 .

(2.26b)

The bilinear transformed matrices Ab1 and Ab 2 can be utilized to evaluate the block eigenvectors of ( Ae − B1e K c1 ) and ( Ae − B2e K c 2 ) associated with the null eigenvalues and the non-null eigenvalues of ( Ae − B1e K c1 ) and ( Ae − B2e K c 2 ) ,respectively. The matrix sign functions and the associated matrix sign minus and sign plus functions of matrices

2-9

Ab1 and Ab 2 can be expressed as

Sign( Ab1 ) = Ab1 ( 2 Ab21 )−1 ,

(2.27a)

Sign( Ab 2 ) = Ab 2 ( 2 Ab22 )−1 ,

(2.27b)

Sign − ( Ab1 ) =

1 ( I n − Sign( Ab1 )) , 2

(2.28a)

Sign − ( Ab 2 ) =

1 ( I n − Sign( Ab 2 )) , 2

(2.28b)

Sign + ( Ab1 ) =

1 ( I n + Sign( Ab1 )) , 2

(2.28c)

and

1 (2.28d) ( I n + Sign( Ab 2 )) , 2 respectively. A fast and stable algorithm for computing the matrix sign function (see Sign + ( Ab 2 ) =

Appendix A) is given as follows : For the order of the desired convergence rate l = 2 , one has G1 (l + 1) =

1 ⎡⎣G1 (l ) + G1−1 (l ) ⎤⎦ , G1 (0) = Ab1 , lim G1 = Sign( Ab1 ), for l = 0,1, 2, l →∞ 2

Then, Ae1 = ζ 1 × sign − ( Ab1 ) + ( Ae − B1e K c1 ) ,

(2.28a)

Ae 2 = ζ 1 × sign − ( Ab 2 ) + ( Ae − B2 e K c 2 ) .

(2.28b)

Notice that Sign − ( Ab1 ) and Sign − ( Ab 2 ) preserve the original eigenspaces of Ae − B1e K c1 and Ae − B2 e K c 2 , respectively. So the eigenvalues of Ae1 and Ae 2 preserve the original ones of Ae − B1e K c1 and Ae − B2 e K c 2 , respectively, except that the null eigenvalues of Ae − B1e K c1 and Ae − B2 e K c 2 are replaced by ζ 1 . Thus, we replace Ae − B1e K c1 and Ae − B2 e K c 2 by Ae1 and Ae 2 , respectively, in (2.25) to have

2-10

Ec1 = − R11−1 B1Te [ Ae−11 ]T C1T Q1 , Ec 2 = − R22−1 B2Te [ Ae−21 ]T C2T Q2 , ⎡E Ec = ⎢ c1 ⎣ 0

(2.30)

0 ⎤ . Ec 2 ⎥⎦

2.1.3 Asymptotic structure

In order to obtain the solutions of the CMARE (2.23), we introduce the following useful lemma [11], [12].

Lemma 2: The CMARE (2.23) is equivalent to the following GCMARE (2.31), respectively AT X + X T A + Q1n − X T S1 X − X T S 2Y − Y T S 2 X = 0 ,

(2.31a)

AT Y + Y T A + Q2 n − Y T S 2Y − Y T S1 X − X T S1Y = 0 ,

(2.31b)

where

⎡ X 00 ⎢ X = ⎢ X 10 ⎢ ⎢⎣ X 20

ε1 X 10T

⎡Y00 ⎢ Y = ⎢Y10 ⎢ ⎢⎣Y20

ε1Y10T

X 11

α X 21

Y11

α Y21

T ⎤ ε 2 X 20

⎥

−1

T ⎥, α X 21

⎥ ⎥⎦

X 22

X = Λ e− 1 X e , X ii = X iiT ,

ε 2Y20T ⎤

⎥ α Y ⎥ , Y = Λ e− 1Y e = Y T Λ e , Y ii = Y i iT , ⎥ Y22 ⎥ ⎦ −1

T 21

i = 0 , 1, 2 .

The GCMARE equation (2.29) can be partitioned as

⎡ f x1 A X + X A + Q1 n − X S1 X − X S 2Y − Y S 2 X = ⎢⎢ f x 2 ⎢⎣ f x 3 T

T

A Y + Y A + Q2n T

T

T

T

T

⎡ f y1 ⎢ − Y S 2 Y − Y S 1 X − X S 1Y = ⎢ f y 2 ⎢ f y3 ⎣ T

T

T

2-11

fx2 fx4 f x5 f y2 f y4 f y5

f x3 ⎤ f x 5 ⎥⎥ = 0 , f x 6 ⎥⎦

(2.32a)

f y3 ⎤ ⎥ f y5 ⎥ = 0 . f y 6 ⎥⎦

(2.32b)

Therefore, from (2.32) can obtain the (2.33) as ε i → +0 , i = 1, 2 , where X

lm

, Y lm ,

lm =00,10, 20,11,21,22 are the zeroth-order solutions T T T A00 X 00 + X 00 A00 + A10T X 10 + X 10T A10 + A20 X 20 + X 20 A20 − X 00 S001 X 00

f x1 ⇒

T T T − X 00 S002Y00 − X 10T S011 X 00 − X 00 S011 X 10 − X 10T S111 X 10 − X 20 S022Y00

(2.33a)

T T − X 00 S022Y20 − X 20 S 222Y20 − Y00 S002 X 00 − Y20T S022 X 00 − Y00 S022 X 20 T − Y20T S 222 X 20 + Q001 = 0,

fx2 ⇒

f x3 ⇒

T X 00 A01 + X 10T A11 + A10T X 11 + α A20 X 21 − ( X 00 S011 X 11 + X 10T S111 X 11 ) T T − α ( X 00 S022Y21 + X 20 S 222Y21 ) − α (Y00 S022 X 21 + Y20T S 222 X 21 ) + Q011 = 0,

T T X 00 A02 + X 20 A22 + A20 X 22 +

1

α

(X α

1

T A10T X 21 −

00

− ( X 00 S022Y22 + X S 222Y22 ) − (Y00 S022 X 22 + Y S T 20

T 20

T T S011 X 21 + X 10T S111 X 21 )

T 222

X 22 ) = 0,

T f x 4 ⇒ A11T X 11 + X 11 A11 − X 11S111 X 11 − α X 21 S 222Y21 − α Y21T S 222 X 21 + Q111 = 0,

T f x 5 ⇒ α X 21 A22 +

1

α

A11T −

T f x 6 ⇒ A22 X 22 + X 22 A22 −

1

α

1

α

(2.33b)

(2.33c)

(2.33d)

T T X 11S111 X 21 − α X 21 S222Y22 − α Y21T S222 X 22 = 0,

T X 21 S111 X 21 − X 22 S 222Y22 − Y22 S 222 X 22 = 0,

(2.33e)

(2.33f)

T T A00 Y00 + Y00 A00 + A10T Y10 + Y10T A10 + A20 Y20 + Y20T A20 − Y00 S 002Y00

f y1 ⇒

T T − Y00 S 001 X 00 − Y10T S 022 Y00 − Y00 S 022Y10 − Y20T S 222Y20 − Y10T S 011 X 00

(2.33g)

T − Y00 S 011 X 10 − Y10T S111 X 10 − X 00 S 001Y00 − X 10T S 011 Y00 − X 00 S 011Y10 T − X 10T S111 Y10 + Q002 = 0,

f y2 ⇒

T Y00 A01 + Y10T A11 + A10T Y11 + α A20 Y21 − α (Y00 S022Y21 + Y20T S 222Y21 )

T Y00 A02 + Y20T A22 + A20 Y22 +

f y3 ⇒

(2.33h)

T − (Y00 S011 X 21 + Y20T S111 X 11 ) − ( X 00 S011Y11 + X 10T S111 Y11 ) = 0,

−

1

α

(Y

00

1

α

A10T Y21T − (Y00 S022Y22 + Y20T S222Y22 )

S011 X + Y S111 X T 21

T 10

T 21

)−

1

α

(X

00

S Y +X S Y T 011 21

T T 10 111 21

f y 4 ⇒ A11T Y11 + Y11 A11 − α Y21T S 222Y21 + Y11S111 X 11 − X 11S111Y11 = 0 , f y 5 ⇒ α Y21T A22 +

1

α

A11T Y21T − α Y21T S222Y22 −

2-12

1

α

T Y11S111 X 21 −

)+Q

022

(2.33i) = 0, (2.33j)

1

α

X 11T S111Y21T = 0 , (2.33k)

T Y22 + Y22 A22 − Y22 S 222Y22 − f y 6 ⇒ A22

1

α

T Y21 S111 X 21 −

1

α

X 21 S111Y21T + Q222 = 0 .

(2.33l)

If Assumption 2 holds, there exist the matrices X 11 ≥ 0 and Y22 ≥ 0 such that the matrices A11 − S111 X 11 and A22 − S 222Y22 are stable, where A11T X11 + X11 A11 − X11S111 X11 + Q11 = 0 and T A22 X 22 + X 22 A22 − X 22 S 222 X 22 + Q22 = 0 . Now we choose X 11 and Y22 to be X 11 and

Y22 , respectively. Then, there exist λx ∈ C and λ y ∈ C such that

(A

− S111 X 11 ν x = λxν x ,

(A

− S 222 X 22 ν y = λ yν y ,

11

22

)

)

where υ x ∈ C n1 and υ y ∈ C

n2

Re ( λx ) < 0 ,

(2.34a)

Re ( λ y ) < 0 ,

(2.34b)

are any vectors.

Using (2.34) we can change (2.33f) and (2.33j) as follows (2.33f) ⇔ υ Ty ( A22 − S 222Y22 )T X 22υ y + υ Ty X 22 ( A22 − S 222Y22 )υ y − ⇔ 2λ yυ Ty X 22υ y −

1

α

1

α

T υ Ty X 21S111 X 21 υy = 0

T υ Ty X 21S111 X 21 υy = 0 ,

(2.35a)

(2.33j) ⇔ υ xT ( A11 − S111 X 11 )T Y11υ x + υ xT Y11 ( A11 − S111 X 11 )υ x − αυ xT Y21S 222Y21Tυ x = 0 ⇔ 2λxυ xT Y11υ x − αυ xT Y21S 222Y21Tυ x = 0 .

(2.35b)

Taking Re(λx ) < 0 and Re(λ y ) < 0 into account, we have X 22 = Y11 = 0 Then, from (2.33e) and (2.33k), we obtain T α X 21 ( A22 − S 222Y22 ) +

1

α Y21T ( A22 − S 222Y22 ) +

1

α α

T =0, ( A11 − S111 X 11 )T X 21

(2.36a)

( A11 − S111 X 11 ) T Y21T = 0 .

(2.36b)

Hence, the unique solution of (2.33f) and (2.33j) are given by X 21 = Y21 = 0 because of the stability of A11 − S111 X 11 and A22 − S 222 X 22 . Thus the parameter α does not appear

2-13

in (2.33), namely, it does not affect the (2.33) in the limit when ε1 and ε 2 tend to zero. Therefore, we obtain the zeroth-order equations (2.37) T T T A00 X 00 + X 00 A00 + A10T X 10 + X 10T A10 + A20 X 20 + X 20 A20 − X 00 S001 X 00 T T T X 00 − X 00 S011 X 10 − X 10T S111 X 10 − X 20 S022Y00 − X 00 S002Y00 − X 10T S011 T T S 222Y20 − Y00 S002 X 00 − Y20T S022 X 00 − Y00 S022 X 20 − X 00 S022Y20 − X 20

(2.37a)

− Y20T S222 X 20 + Q001 = 0, X 00 A01 + X 10T A11 + A10T X 11 − ( X 00 S 011 X 11 + X 10T S111 X 11 ) + Q011 = 0 ,

(2.37b)

T T X 00 A02 + X 20 A22 − ( X 00 S 022Y22 + X 20 S 222Y22 ) = 0 ,

(2.37c)

A11T X 11 + X 11 A11 − X 11S111 X 11 + Q111 = 0 ,

(2.37d)

T T A00 Y00 + Y00 A00 + A10T Y10 + Y10T A10 + A20 Y20 + Y20T A20 − Y00 S002Y00 T T Y00 − Y00 S022Y20 − Y20T S222Y20 − Y10T S011 X 00 − Y00 S001 X 00 − Y20T S022 T Y00 − Y00 S011Y10 − Y00 S011 X 10 − Y10T S111 X 10 − X 00 S001Y00 − X 10T S011

(2.37e)

− X 10T S111Y10 + Q002 = 0, Y00 A01 + Y10T A11 − (Y00 S011 X 11 + Y10T S111 X 11 ) = 0 ,

(2.37f)

T Y00 A02 + Y20T A22 + A20 Y22 − (Y00 S022Y22 + Y20T S 222Y22 ) + Q022 = 0 ,

(2.37g)

T A22 Y22 + Y22 A22 − Y22 S 222Y22 + Q222 = 0 .

(2.37h)

The Nash equilibrium strategies for the MSPS will be studied under the following basic assumption, so that one can apply the proposed method to the nonstandard MSPS. Assumption 3: The Hamiltonian matrices T i i i , i = 1, 2 are nonsingular, where

⎡ A Tiii := ⎢ ii ⎣ −Qiii

− Siii ⎤ . − AiiT ⎦⎥

(2.38)

Under Assumptions 2 and 3, we obtain the following zeroth-order equations: AsT X 00 + X 00 As + Qs1 − X 00 S s1 X 00 − X 00 S s 2Y00 − X 00 S s 2 X 00 = 0 ,

(2.39a)

AsT Y00 + Y00 As + Qs 2 − Y00 S s 2Y00 − Y00 S s1 X 00 − X 00 S s1Y00 = 0 ,

(2.39b)

2-14

A11T X 11 + X 11 A11 − X 11S111 X 11 + Q111 = 0 ,

(2.39c)

T A22 Y22 + Y22 A22 − Y22 S 222Y22 + Q222 = 0 ,

(2.39d)

X 10 = − Dx−11T DxT01 X 00 − Dx−11T N xT01 ,

(2.39e)

Y10 = − Dx−11T DxT01Y00 ,

(2.39f)

T X 20 = − Dy−22 DTy 02 X 00 ,

(2.39g)

T T Y20 = − Dy−22 DTy 02Y00 − Dy−22 N Ty 02 ,

(2.39h)

where 1 As = A00 − Dx 01 Dx−111 A10 − D02 Dy−22 A20 + ( S011 − D01 Dx−111 S111 ) Dx−11T N xT01 T 1 + ( S022 − D02 Dy−22 S222 ) Dy−22 N Ty 02 ,

T S s1 = S001 − Dx 01 Dx−111 S011 − S011 Dx−11T DxT01 + Dx 01 Dx−111 S111 Dx−11T DxT01 ,

1 T T 1 T S s 2 = S002 − Dy 02 Dy−22 S022 − S022 Dy−22 DTy 02 + Dy 02 Dy−22 S222 Dy−22 DTy 02 ,

Qs1 = Q001 − A10T Dx−11T N xT01 − N x11 Dx−111 A10 − N x 01 Dx−111 S111 Dx−11T N xT01 , T T 1 1 T Qs 2 = Q002 − A20 Dy−22 N Ty 02 − N y 02 Dy−22 A20 − N y 02 Dy−22 S222 Dy−22 N Ty 02 ,

Dx 01 = A01 − S011 X 11 , Dx11 = A11 − S111 X 11 ,

Dy 02 = A02 − S022Y22 , Dy 22 = A22 − S222Y22 .

Lemma 3: The matrices

As , S s1 , S s 2 , Qs1 , and Qs 2 do not depend on the matrices

X 11 , X 21 , Y21 , and Y22 . Therefore, their matrices can be calculated by using the following

2-15

Hamiltonian matrices:

⎡ As ⎢* ⎣

* ⎤ ⎡ A00 = − AsT ⎥⎦ ⎢⎣ *

* ⎤ −1 −1 − T011T111 T101 − T022T222 T202 , T ⎥ − A00 ⎦

− S s1 ⎤ ⎡ * −1 = T001 − T011T111 T101 , ⎢ −Q ⎥ * ⎦ ⎣ s1 ⎡ * ⎢ −Q ⎣ s2

(2.40b)

−Ss 2 ⎤ = T − T T −1 T , * ⎥⎦ 002 022 222 202

− S001 ⎤ ⎡ A , T001 = ⎢ 00 T ⎥ ⎣ −Q001 − A00 ⎦

(2.40a)

(2.40c)

− S011 ⎤ ⎡ A T011 = ⎢ 01 , T ⎥ ⎣ −Q011 − A10 ⎦

T − S111 ⎤ ⎡ A ⎤ − S011 ⎡ A , T111 = ⎢ 11 T101 = ⎢ 10T , T ⎥ T ⎥ ⎣ −Q111 − A11 ⎦ ⎣ −Q011 − A01 ⎦

⎡ A T002 = ⎢ 00 ⎣ −Q002

− S002 ⎤ ⎡ A , T022 = ⎢ 02 T ⎥ − A00 ⎦ ⎣ −Q022

− S022 ⎤ , T ⎥ − A20 ⎦

⎡ A T202 = ⎢ 20T ⎣ −Q022

T ⎤ − S022 ⎡ A , T222 = ⎢ 22 T ⎥ − A02 ⎦ ⎣ −Q222

− S222 ⎤ , T ⎥ − A22 ⎦

where * stands for an appropriate matrix. Moreover, we can change the form of the solutions X 10 , X 20 , Y10 , and Y20 .

0⎤ ⎡I −1 ⎡⎣ X 10 Y10 ⎤⎦ = ⎡⎣ X 11 − I n1 ⎤⎦ T111 T101 ⎢ n 0 ⎥, ⎣ X 00 Y00 ⎦

(2.41a)

⎡ 0 I n0 ⎤ −1 T202 ⎢ − I n 2 ⎤⎦ T222 ⎥. ⎣ X 00 Y00 ⎦

(2.41b)

⎡⎣ X 20 Y20 ⎤⎦ = ⎡⎣Y22

The following theorem will establish the relation between the solutions X and Y and the zeroth-order solutions X lm and Ylm , lm = 00,10, 20,11, 21, 22 for the reduced-order equations (2.38)

2-16

Lemma 4 : By [14], suppose that ⎡ Aˆ T ⊗ I + I ⊗ Aˆ − ⎡( S X )T ⊗ I + I ⊗ ( S X )T ⎤ ⎤ n0 n0 s n0 n0 s 2 00 ⎢ s ⎣⎢ s 2 00 ⎦⎥ ⎥ det ⎢ ⎥ ≠ 0, T T T ⎡ ⎤ ˆ ˆ ⎢ − ⎢( S s1Y00 ) ⊗ I n 0 + I n 0 ⊗ ( S s1Y00 ) ⎥ + As ⊗ I n 0 + I n 0 ⊗ As ⎥ ⎦ ⎣ ⎣ ⎦

(2.42)

where Aˆ s := As − S s1 X 00 − S s 2Y00 and the matrix Aˆ s is stable . Under Assumptions 1-3, the GCMARE (2.31) admits the stabilizing solutions X and Y such that these matrices possess a power series expansion at µ = 0 . That is ⎡ X 00 ⎢ X = ⎢ X 10 ⎢ X 20 ⎣ ⎡Y00 ⎢ Y = ⎢Y10 ⎢Y20 ⎣

0 X 11 0

0⎤ ⎥ 0⎥ + O ( µ ) = X + O ( µ ) , 0 ⎥⎦

0⎤ ⎥ 0 0 ⎥ + O( µ ) = Y + O( µ ) . 0 Y11 ⎥⎦

(2.43a)

0

(2.43b)

2.1.3 Newton’s method

In order to improve the convergence rate of the Lyapunov iterations [7], the following new algorithm based on the Newton’s method is used

φ ( n )T P ( n +1) + P ( n +1)T φ ( n ) − Θ( n ) P ( n +1) ℑ − ℑP ( n +1) Θ( n ) + Ξ ( n ) = 0, ⎧φ ( n ) T X ( n +1) + X ( n +1) T φ1( n ) − Θ (2n ) T Y ( n +1) − Y ( n +1) T Θ (2n ) + Ξ 1( n ) = 0 ⇔ ⎨ 1( n ) T ( n +1) + Y ( n +1) T φ 2( n ) − Θ 1( n ) T X ( n +1) − X ( n +1) T Θ 1( n ) + Ξ (2n ) = 0 ⎩φ 2 Y

where

φ

(n)

:= A − SP

(n)

⎡φ1( n ) 0 ⎤ − ℑSP ℑ = ⎢ , (n) ⎥ ⎣ 0 φ2 ⎦ (n)

2-17

n=0,1, (2.44)

⎡ 0 Θ( n ) := S ℑΡ ( n ) = ⎢ ( n ) ⎣Θ2 Ξ

(n)

:= Q + P

( n )T

⎡ Φ (00n )i ⎢ (n) Φ i( n ) := ⎢ Φ10 i ( n ⎢Φ 20)i ⎣ ⎡Ξ (00n )i ⎢ Ξ i( n ) := ⎢ Ξ (01n i) ⎢Ξ (02n )i ⎣

(n) Φ11 i

Φ (21n i) (n) Ξ 01 i

Ξ

(n) 11i (n) 12 i

Ξ

⎡A 0⎤ A := ⎢ ⎥, ⎣ 0 A⎦

0⎤ , S 2 ⎥⎦

+ ℑP

( n )T

(n) ⎤ Φ 02 i (n) ⎥ Φ12i ⎥ , Φ (22n )i ⎥⎦

(n) Φ 01 i

⎡ X 00( n ) ⎢ X ( n ) := ⎢ X 10( n ) ⎢ (n) ⎢⎣ X 20

⎡S S := ⎢ 1 ⎣0

SP

(n)

Θ1( n ) ⎤ ⎥, 0 ⎦

(n) ⎤ Ξ 02 i (n) ⎥ Ξ12i ⎥ , Ξ (22n )i ⎥⎦

S ℑP

(n)

⎡ Θ(00n )i ⎢ (n) Θi( n ) := ⎢ Θ10 i ( n ⎢Θ 20)i ⎣

P

(n)

⎡ Ξ1( n ) ℑSP ℑ = ⎢ ⎣ 0 (n)

(n) Θ01 i (n) Θ11 i

Θ(21n i)

0 ⎤ ⎥, Ξ (2n ) ⎦

(n) ⎤ Θ02 i (n) ⎥ Θ12i ⎥ , Θ(22n )i ⎥⎦

⎡ X ( n) 0 ⎤ := ⎢ ⎥, Y (n) ⎦ ⎣ 0

ε1 X 10( n )T

ε 2 X 20( n )T ⎤

X 11( n )

α X 21( n )T ⎥ ,

α X 21( n )

+P

( n )T

−1

X 22( n )

⎥ ⎥ ⎥⎦

⎡Y00( n ) ⎢ Y ( n ) := ⎢Y10( n ) ⎢ (n) ⎢⎣Y20

ε1Y10( n )T

ε 2Y20( n )T ⎤

Y11( n )

α Y21( n )T ⎥ ,

α Y21( n )

⎥

−1

Y22( n )

⎥ ⎥⎦

⎡Q 0 ⎤ Q := ⎢ 1 ⎥, ⎣ 0 Q2 ⎦ ⎡0 ℑ := ⎢ ⎣IN

IN ⎤ . 0 ⎥⎦

And the initial condition P (0) has the following form

⎡ X (0) P (0) = ⎢ ⎣ 0

T ⎡ X 00 ε1 X 10T ε 2 X 20 0 ⎢ X 11 0 0 ⎢ X 10 0 0 0 0 ⎤ ⎢ X 20 =⎢ (0) ⎥ 0 0 Y00 Y ⎦ ⎢ 0 ⎢ 0 0 0 Y10 ⎢ 0 0 Y20 ⎣⎢ 0

0 0 0

ε1Y10T 0 0

0 ⎤ ⎥ 0 ⎥ 0 ⎥ ⎥. ε 2Y20T ⎥ 0 ⎥ ⎥ Y22 ⎦⎥

(2.45)

Note that the considered algorithm (2.44) is original. The new algorithm (2.44) can be constructed by setting

P ( n +1) = P ( n ) + ∆P ( n ) and neglecting

O(∆P ( n )T ∆P ( n ) )

term.

Newton’s method is well known and is widely used to find a solution of the algebraic

2-18

equations, and its local convergence properties are well understood. First we show that the algorithm (2.43) is equivalent to the Newton’s method. Now, let us define the following matrix function: F ( P ) := AT P + PT A + Q − PT SP − ℑPT S ℑP − PT ℑSPℑ

(2.46)

where ⎡X P := ⎢ ⎣0

0⎤ . Y ⎥⎦

Taking the partial derivative of the function F ( P ) with respect to P yields ∇F ( P ) :=

∂vecF ( P) ∂ ( vecP )

T

(

P=P

)

(

)

T T = A − SP − ℑSPℑ ⊗ I 2 N + ⎡ I 2 N ⊗ A − SP − ℑSPℑ ⎤ U 2 N 2 N ⎢⎣ ⎥⎦

(

− S ℑP

)

T

(

)

(2.47)

T

⊗ ℑ − ⎡ ℑ ⊗ S ℑP ⎤ U 2 N 2 N . ⎣ ⎦

Taking the vec-operator transformation on both sides of (2.44) and (2.46), we obtain

⎡( Φ ( n )T ⊗ I 2 N ) U 2 N 2 N + I 2 N ⊗ Φ ( n )T ⎤ vecP ( n +1) ⎣ ⎦ − ⎡⎣( Θ( n )T ⊗ ℑ) U 2 N 2 N + ℑ ⊗ Θ( n )T ⎤⎦ vecP ( n +1) − vecΞ ( n ) = 0, vecF ( P ( n ) ) = ⎡⎣( Φ ( n )T ⊗ I 2 N ) U 2 N 2 N + I 2 N ⊗ Φ ( n )T ⎤⎦ vecP ( n )

− ⎡⎣( Θ( n )T ⊗ ℑ)U 2 N 2 N + ℑ ⊗ Θ( n )T ⎤⎦ vecP ( n ) − vecΞ ( n ) = 0.

(2.48a)

(2.48b)

Subtracting (2.48a) from (2.48b) and noting that ∇F ( P ( n ) ) := ⎡⎣( Φ ( n )T ⊗ I 2 N ) U 2 N 2 N + I 2 N ⊗ Φ ( n )T ⎤⎦ − ⎡⎣( Θ( n )T ⊗ ℑ) U 2 N 2 N + ℑ ⊗ Θ( n )T ⎤⎦

,

(2.49)

we have −1

vecP ( n +1) = vecP ( n ) − ⎡⎣∇F ( P ( n ) ) ⎤⎦ vecF ( P ( n ) ) ,

(2.50)

which is the desired result. We are concerned with good choices of the starting points which guarantee to find a required solution of a given GCMARE (2.31). Our new idea is to set the initial conditions

2-19

to the matrix (2.45). The fundamental idea is based on

P − P (0) = Ο ( µ

) which

is

derived from (2.43). Consequently, we will get the required solution with rate of the quadratic

convergence

via

the

Newton’s

method.

Moreover,

using

the

Newton-Kantorovich theorem [20], we will also prove the existence of the unique solution for the GCMARE (2.31). The main result of this section is as follows.

Lemma 5: by [14] under Assumptions 1–3, the new iterative algorithm (2.44) converges to the exact solution P* of the GCMARE (2.31) with the rate of quadratic convergence. Furthermore, the unique bounded solution P ( n ) of the GCMARE (2.31) is in the neighborhood of the exact solution P (*) . That is, the following conditions are satisfied

( ),

P ( n ) − P* ≤ Ο µ

P ( n ) − P* ≤

2n

(2.51a)

n =0,1,

1 ⎡ 1 − 1 − 2θ ⎤⎦ , β℘ ⎣

n =0,1,

(2.51b)

where

⎡X* 0 ⎤ P=P =⎢ , *⎥ ⎣0 Y ⎦ *

℘:= 6 S ,

β := ⎣⎡∇F ( P (0) ) ⎦⎤

η := ⎡⎣∇F ( P (0) ) ⎤⎦ F ( P ) := ( f x1 ,

−1

−1

,

⋅ F ( P (0) ) ,

, f x 6 , f y1 ,

θ := βη℘,

, f y6 ) ,

P (0) := P0 .

Remark 1: Using the Newton–Kantorovich theorem instead of the implicit function theorem, we can also prove (2.43). That is, the structure of the solutions X and Y are established by setting n = 0 in (2.51a).It should be noted that the asymptotic structure for the GCMARE (2.31) is established by using the Newton–Kantorovich theorem which is

2-20

different from implicit function theorem. As a result, the condition of the small parameter µ

for existence of the solutions of the GCMARE (2.31) would be clear.

Remark 2: In view of [17], we know that the solution of the GCMARE (2.31) is not unique and several nonnegative solutions exist. In this paper, it is very important to note that if the sufficient condition (2.42) holds, the new algorithm (2.44) converge to the desired positive semidefinite solution in the same way as the Lyapunov iterations [3]. Remark 3: In order to obtain the initial condition (2.45), we have to solve the CARE (2.39a) and (2.39b) which are independent of the perturbation parameters. In this situation, we can also apply the Newton’s method to these equations (2.39a) and (2.39b). The resulting algorithm is given by

(A −S s

X 00( k ) − S s 2Y00( k ) ) X 00( k +1) + X 00( k +1) ( As − S s1 X 00( k ) − S s 2Y00( k ) ) T

s1

− X 00( k ) S s 2Y00( k +1) − Y00( k +1) S s 2 X 00( k ) + X 00( k ) S s1 X 00( k )

(2.52a)

+ X 00( k ) S s 2Y00( k ) + Y00( k ) S s 2 X 00( k ) + Qs1 = 0,

(A −S s

X 00( k ) − S s 2Y00( k ) ) Y00( k +1) + Y00( k +1) ( As − S s1 X 00( k ) − S s 2Y00( k ) ) T

s1

− Y00( k ) S s1 X 00( k +1) − X 00( k +1) S s1Y00( k ) + Y00( k ) S s 2Y00( k )

(2.52b)

+ Y00( k ) S s1 X 00( k ) + X 00( k ) S s1Y00( k ) + Qs 2 = 0, with AsT X 00(0) + X 00(0) As − X 00(0) S s1 X 00(0) + Qs1 = 0

(2.52a)

AsT Y00(0) + Y00(0) As − Y00(0) S s 2Y00(0) + Qs 2 = 0 ,

(2.52b)

In the rest of this section, we explain the method for solve the generalized cross-coupled multiparameter algebraic Lyapunov equations (GCMALE) (44). So that we can completely solve the generalized cross-coupled multiparameter algebraic Riccati equations (GCMARE) (2.31) and get the controller as follow u (t ) = − K c x (t ) + E c r (t )

(2.54)

where

2-21

K c 1 = R 1−11 B 1Te X e ,

Ec1 = − R11−1 B1Te [( Ae1 )−1 ]T C1T Q1 ,

K c 2 = R 2−21 B 2T e Y e ,

Ec 2 = − R22−1 B2Te [( Ae 2 )−1 ]T C2T Q2 ,

⎡K ⎤ K c = ⎢ c1 ⎥ , ⎣ K c2 ⎦

⎡E Ec = ⎢ c1 ⎣0

0 ⎤ . Ec 2 ⎥⎦

Then we have digital controller as follows ud ( kT ) = − K d xd ( kT ) + Ed r * ( kT )

(2.55)

where K d = ( I m + K c H ) −1 K c G , Ed = ( I m + K c H ) −1 Ec

in which G = exp( AeT ), 3 ⎧ T2 2 T + Ae + …]Be , for singular Ae ⎪[ I T + Ae H =⎨ n . 2! 3! ⎪(G − I ) A−1 B , for nonsingular Ae n e e ⎩

In order to improve the performance of the tracker, tuning the tracker with Chaos-Evolutionary-Programming algorithm (CEPA) may be the optimum method. In next section, tuning the forward gain Ec via a CEPA tuning is applied to establish an effective tracker. Now, let’s study how the chaos-evolutionary programming can be applied to the tracker. First we should know what the chaos system selected for optimization.

2-22

2.2. The chaos-evolutionary-programming approach

Because of the (2.30), we find the forward gain Ec1 and Ec 2 . The forward gain Ec1 and Ec 2 that we find are just acceptive but not good enough. So we use the algorithm chaos evolutionary programming approach (CEPA) to optimal the parameter Ec1 and Ec 2 . We introduce the CEPA algorithm as following sections.

2.2.1. The chaos system selected for optimization

The chaotic equation for chaos optimization algorithm can be selected as the logistic mapping, namely

t k +1 = f ( µ , t k ) = µ t k (1 − t k )

k = 1, 2,

,N ,

(2.56)

where µ is a control parameter. It is easy to find that (2.56) becomes chaotic and stochastic as µ = 4 . When µ = 4 , the equation is described as

tk +1 = 4tk (1 − tk ) .

(2.57)

The chaotic space of (2.57) is [0 1]. The system is easy to be affected by initial conditions [18]. A slight change of initial conditions will case a violent distinction in the long-time behavior of chaotic system. Besides, the chaotic motion can also go non-repeatedly through every state in a certain domain.

2.2.2. The chaos-evolutionary-programming approach

If the optimization problems are continuous problems rather than discrete problems and the constraints of the variables are known, the optimization problems can be described as min f ( xi ) or max f ( xi ), i = 1, 2,

, n, for ai ≤ xi ≤ bi .

2-23

(2.58)

Suppose that the natural numbers are represented in the scale of notation with radix R, so that n = a0 + a1 R + a 2 R 2 + + a m R m , 0 ≤ a i ≤ R .

(2.59)

Write the digits of these numbers in the reverse order, preceded by a decimal point. This gives the number

φ R ( n ) = a 0 R −1 + a1 R −2 + + a m R − m−1 .

(2.60)

Halton [19] extended the two-dimensional result of Van Der Corput [21] to ρ -dimensions, when R1 ,R2,…,R ρ are mutually coprime. We show a binary scale and an illustration in Table 2.1 and Table 2.2.

Table 2.1 Natural numbers in binary scale (R=2) n (decimal)

n(binary)

φ 2 ( n ) (binary)

φ 2 ( n ) (decimal)

1

1

0.1

0.5

2

10

0.01

0.25

3

11

0.11

0.75

4

100

0.001

0.125

5

101

0.101

0.625

6

110

0.011

0.375

7

111

0.111

0.875

8

1000

0.0001

0.0625

2-24

Table 2.2 Quasi-random sequences

φ R ( n ) R=2

3

5

7

11

13

17

…

N=1

0.5000

0.333

0.2000

0.1429

0.0909

0.0769

0.0588

…

2

0.2500

0.6667

0.4000

0.2857

0.1818

0.1538

0.1176

…

3

0.7500

0.1111

0.6000

0.4286

0.2727

0.2308

0.1765

…

48

0.0469

0.1975

0.7680

0.9796

0.3967

0.7101

0.8304

…

49

0.5469

0.5309

0.9680

0.0029

0.4876

0.7870

0.8893

…

50

0.2969

0.8642

0.0160

0.1458

0.5785

0.8639

0.9481

…

Since φ R ( n ) < 1 , to satisfy this range, scaling any varying parameter (e.g., a real number η from its range ⎡⎣η η ⎤⎦ to [ 0 1] ) is required. Let the interval real ( ℑℜ ) matrix X ∈ ℑℜ n×m be a set of degenerate real matrices defined by X = [ L, U ] = {[ x ij ]| lij ≤ x ij ≤ uij ; 1 ≤ i ≤ n, 1 ≤ j ≤ m} ,

(2.61)

where L and U are constant real matrices. We introduce the variable ηij , 0 ≤ ηij ≤ 1 such that xij = lij + ηij (uij − lij ) ,

and use the notation η = [η11 ,

,η1m ,η 21 ,

(2.62) ,η2 m ,η n1 , ηnm ] . Then the interval matrix X can

be denoted as X (η ) . Let η11 = φ2 (n) 、η12 = φ3 (n) 、η13 = φ5 (n) …, and so on, to construct the desired initial population of size N (e.g., N = 50). Then CEPA can be used to optimize these problems. The developed CEPA for the minimal or maximal principle is described as follows: 1) Individual population:

2-25

Choose individual population based on the quasi-random sequence (QRS) to form an initial population P0 = [P1,P2 ,

,PN ] of size N by initializing each ρ -dimensional

solution vector Pi in S. 2) Objective function: Assign each Pi , i = 1,

, N an objective function score. Arrange Pi , i = 1,

, N in

the descending order, starting from the best one generated from the objective function score. 3) Fitness function: Assign each sorted Pi , i = 1,

, N a fitness function (FF) score to weigh those

high-quality individuals in the pool of individuals based on the obtained objective function scores. Search some P* in the solution Pi , i = 1,

, N , so that the objective

function (OF) value OF ( Pi ) is maximal, using

⎛ ⎞ β −β ⎟(OF ( P ) − OF ( P )) + β . FF (OF ( Pi )) = ⎜ i i ⎜ OF ( P ) − OF ( P ) ⎟ i i ⎝ ⎠ On the other hand, search some P* in the solution Pi , i = 1,

(2.63)

, N , so that the

objective function value OF ( Pi ) is minimal, using −1

⎡⎛ ⎤ ⎞ β −β ⎟(OF ( P ) − OF ( P )) + β ⎥ . FF (OF ( Pi )) = ⎢⎜ i i ⎢⎜⎝ OF ( Pi ) − OF ( Pi ) ⎟⎠ ⎥ ⎣ ⎦

2.64)

This function linearly maps the real-valued space [OF ( Pi ), OF ( Pi )] to any appropriate specified space, [ β , β ] , where β > 0 , for weighting the objective function scores. Hence, the better an individual is , the higher the objective function score that it will have. 4) Probability function:

2-26

Calculate the probability function (PF) score of each Pi , i = 1,

, N using the fitness

function score: PF ( FF ( Pi )) := PF ( Pi ) =

FF ( Pi )

∑i =1 FF ( Pi ) N

.

(2.65)

This equation guarantee that each individual a non-zero selection probability. It is an essential factor to determine the preservative or extinctive of the individual. 5) Mutation: Mutate each Pi , i = 1,

, N based on statistics to increase the population size from N

to 1.6N, assign Pi+ N the following value Pi + N , j := Pi , j (1 + sgn( Ν (0,1))γ (1 − PF ( Pi )))

for i = 1, 2,

, 0.6 N

(2.66a)

and values of 1.6 N to 2 N are produced from the best Pbest , j as follows: P1.6 N + k , j := Pbest , j (1 + k % × α1 )

for

k = 1, 2,

, 0.2 N ,

(2.66b)

P1.8 N + k , j := Pbest , j (1 − k % × α1 )

for

k = 1, 2,

, 0.2 N ,

(2.66c)

where α1 is a weighting factor, Pi, j is the j-th element in the i-th individual, N ( µ , σ 2 ) is the Gaussian random variable with mean µ and variance σ 2 , γ is a

weighting factor for the percentage change of Pi , j , and sgn(⋅) is the standard sign function. Whenever Pi + N , j ∈ / [ P j , P j ] , some modification is required ⎧P Pi + N , j := ⎨ j ⎩P j

if if

Pi + N , j < P j . Pi + N , j > P j

(2.67)

Properly adjusting the weighting factor γ can possibly avoid the undesired situation Pi + N , j ∈ / [ P j , P j ] . It is notable that γ heavily dominates the convergence rate of the

EP. 6) Select:

2-27

Calculate the objective function score of each Pi+ N , i = 1, function scores of Pi , i = 1,

,2 N . Record Pi , i = 1,

, N . Rank the objective

,2 N , in the descending order,

starting from the best individual in the pool of the population (proportional selection). The first N individuals are selected for the next generation, in which the top one of each * generation (elitist model), denoted Pg,i , always survives and is selected for the next

* generation. Whenever Pg,i is no longer the best during the evolutionary process, update

it by the newly generated best one. 7) Penalty: Tune γ in the following way to further avoid the search be trapped into a local extreme ⎧γ ⎪ γ := ⎨1.5γ ⎪0.5γ ⎩

if

| OF ( Pg*−1,i ) − OF ( Pg*,i ) |> ε

if

| OF ( Pg*−1,i ) − OF ( Pg*,i ) |≤ ε

if

* g −2,i

| OF ( P

,

) − OF ( P ) |≤ ε and * g ,i

* g −1,i

| OF ( P

(2.68)

) − OF ( P ) |≤ ε * g ,i

where ε is some tolerable error bound and g is the generation index. Then, go to Step 2) and continue until the desired extreme value OF ( Pg*,i ) cannot be further improved and/or the allowable generation is obtained. 8) EP termination condition: After some generations without improvement, stop the EP algorithm. 9) Then, carry the proposed chaotic search out: a) Generate N − 1 chaotic variables by (2.57). b) Change N − 1 chaotic variables to real variables in ⎡⎣ Pg*, j − θ , Pg*, j + θ ⎤⎦ . Furthermore, amplify the ergodic areas of the N − 1 chaotic variables to the variance ranges of real variables by (2.69). Set P1, j = p*g , j , which is the best population. Pk +1, j = ( p*g , j − θ ) + 2 × θ × tk , j

for k = 1, 2, …, N − 1 ,

2-28

(2.69a)

PN + k , j := Pg*, j (1 + k % × α 2 )

for k = 1, 2,

, 0.2 N ,

(2.69b)

P1.2 N + k , j := Pg*, j (1 − k % × α 2 )

for k = 1, 2,

, 0.2 N ,

(2.69c)

where θ = ξ Pg*, j , ξ denote some percentage of Pg*, j and α 2 is a weighting factor. c) Calculate the value of the objective function OF ( Pk ) . d) Sort the minimum or maximum value of OF ( Pk ) . Then update the P1, j from OF ( Pk ) .

e) Calculate tk +1, j = µtk , j (1 − tk , j ) for k := k + N − 1 , then go to (2.69) and repeat n1 times. f) Go to Step b) and repeat n2 times. When OF ( Pk ) holds, stop it.

The interval matrix EcI is set as a real matrix with an appropriate dimension defined as EcI = ⎡⎣ E

{

}

E ⎤⎦ = ⎡⎣eij ⎤⎦ | e ij ≤ eij ≤ eij ∈ IR n×n ,

(2.70)

where eij = eij + eij , e ij = eij − eij , and eij denotes the absolute value of eij . Then, the EcI is used as the CEPA’s initial interval.

2.3

An illustrative example

Example 2.1: In this section, to show the effect of the proposed the design to the linear quadratic Nash games based tracker for multiparameter singularly perturbed systems, let consider the following system

2-29

2

2

k =0

k =1

x0 (t ) = ∑ A0i xi (t ) + ∑ B0i ui (t ) , x0 (t ) = x00 ,

ε i xi (t ) = Ai 0 x0 (t ) + Aii xi (t ) + Bii ui (t ) ,

(2.71)

xi (t ) = xi0 ,

yi (t ) = Ci 0 x0 (t ) + Cii xi (t ) = Ci (t ) , ⎡ x0 (t ) ⎤ x(t ) = ⎢⎢ x1 (t ) ⎥⎥ , ⎢⎣ x2 (t ) ⎥⎦

i =1, 2 .

where 0 1 ⎤ ⎡ 0 0⎤ ⎡ 0 0⎤ ⎢ ⎥ ⎢ 0 0⎥ ⎥ −1 ⎥ 4.5 ⎢ 0 0⎥ ⎢ ⎥ −0.1⎥ , A01 = ⎢0.1 0 ⎥ , A02 = ⎢ 0 0 ⎥ , 0 ⎢ ⎥ ⎢ ⎥ ⎥ −0.05 0.1 ⎥ ⎢ 0 0⎥ ⎢ 0.1 0 ⎥ ⎢⎣ 0 0 ⎥⎦ ⎢⎣ 0 0 ⎥⎦ −32.7 0 ⎥⎦ 0 0⎤ 0 0⎤ ⎡0 0 0 , A20 = ⎢ ⎥ ⎥, 0 0⎦ ⎣0 0 0 −0.4 0 ⎦ ⎡0⎤ T B01 = B02 = [ 0 0 0 0 0] , B11 = B22 = ⎢ ⎥ , ⎣0.1⎦ ⎡0 ⎢0 ⎢ A00 = ⎢0 ⎢ ⎢0 ⎢⎣0 ⎡0 A10 = ⎢ ⎣0

0 4.5 0 0 0 −0.05 0 0 0 32.7 0 0 0 −0.4

x0 = [ 0 0 0 0 0 0.5 0.5 0 0] . T

⎡cos(4t ) ⎤ The small parameters are chosen as ε1 = ε 2 = 0.01 . The reference input r (t ) = ⎢ ⎥ . ⎣ sin(4t ) ⎦ In this simulation, we select the sampling period as T = 0.02s and the final simulation time as t = 1s . The initial values and parameters are R11 = R22 = 3.3 in (2.19), Q1 = Q2 = 105 in (2.19), Q1n = C1T Q1C1 = 105 diag (1 1 1 1 1 1 1 0 0 ) in (2.23a),

Q2 n = C2T Q2C2 = 105 diag (1 1 1 1 1 0 0 1 1) in (2.23b), Eigenvalue of Ae − B1e K c1 and Ae − B2 e K c 2 are shown as follows

( −1740.8

−10.35 −6.32 −4.29 −0.30 ± 2.79i −0.48 ± 0.32i 0 ) ,

2-30

( −1740.8

−10.35 −6.33 −4.29 −0.30 ± 2.79i −0.478.32i 0 ) ,

and ζ 1 = 10 −4 in (2.26). Then from (2.30) get the result are Ec1 = 174.02 = Ec11 , Ec 2 = 174.02 = Ec 22 , Ec12 = 0 , Ec 21 = 0 . So the initial value of CEPA are N=50, β = [1 ,10] ,

γ = 0.3 , α1 = 1 , ε = 0.01 , α 2 = 1 , n1 = 5 , n2 = 2 , generation=8. ⎡ ⎡ Ec11 − Ec11 EcI = ⎢ ⎣ ⎢⎣ ⎡⎣ Ec 21 − Ec 21 ⎡[ 0 348.04] =⎢ ⎣ [ 0 0]

Ec11 + Ec11 ⎤⎦

⎡⎣ Ec12 − Ec12 Ec 21 + Ec 21 ⎤⎦ ⎡⎣ Ec 22 − Ec 22

[0

0] ⎤ . [0 348.04]⎥⎦

Because of the interval ⎡⎣ Ec12 − Ec12

[0

Ec12 + Ec12 ⎤⎦ ⎤ ⎥ Ec 22 + Ec 22 ⎤⎦ ⎥⎦

Ec12 + Ec12 ⎤⎦ and ⎡⎣ Ec 21 − Ec 21

Ec 21 + Ec 21 ⎤⎦ are

0] and [ 0 0] , So by experience, it is chosen as ⎡ ⎡⎣ Ec11 − Ec11 Ec11 + Ec11 ⎤⎦ ⎡⎣ Ec12 − 0.5 Ec11 EcI = ⎢ ⎡⎣ Ec 22 − Ec 22 ⎢⎣ ⎡⎣ Ec 21 − 0.5 Ec 22 Ec 21 + 0.5 Ec 22 ⎤⎦ ⎡ [ 0 348.04] [ −87.01 87.01]⎤ . =⎢ [0 348.04] ⎥⎦ ⎣[ −87.01 87.01]

Ec12 + 0.5 Ec11 ⎤⎦ ⎤ ⎥ Ec 22 + Ec 22 ⎤⎦ ⎥⎦

The simulation results by CEPA are shown as follows ⎡153.94 −23.656 ⎤ Ec = ⎢ ⎥, ⎣ 18.95 142.79 ⎦ 84.71 173.32 −10.88 −0.03 ⎤ ⎡ 174.14 −2.69 2915.3 −1196.6 48.81 Kc = ⎢ ⎥, ⎣ −2.6948 174.14 −1195.5 2914.3 −48.755 −10.90 −0.03 84.74 173.32 ⎦ ⎡ 4.63 −0.69 ⎤ Ed = ⎢ ⎥, ⎣0.58 4.29 ⎦ 2.47 4.5 −0.35 −0.02 ⎤ ⎡ 5.24 −0.06 85.1 −36.62 1.32 . Kd = ⎢ 4.5 ⎥⎦ ⎣ −0.06 5.24 -36.59 85.07 −1.32 −0.36 −0.02 2.47 The simulations are shown on following figures

2-31

(a)

(b) Fig. 2.3. Time responses of ri (t ) , yci (t ) and ydi (t ) by proposed controllers without CEPA: (a) continuous-time system (b) sampled-data system

2-32

Fig.2.4. The time series of uci (t ) and udi (t ) without CEPA

(a)

2-33

(b) Fig. 2.5. Time responses of ri (t ) , yci (t ) and ydi (t ) by proposed controllers with CEPA: (a) continuous-time system (b) sampled-data system

Fig. 2.6. The time series of uci (t ) and udi (t ) with CEPA

2-34

Chapter 3 Conclusion

A tracker for the linear quadratic Nash game of multiparameter singularly perturbed sampled-data systems is newly proposed in this thesis. The proposed control methodology needs to solve generalized-cross-coupled multiparameter algebraic Riccati equation (GCMARE). First, the algorithm based on Newton’s method is used to improve the convergence rate of the Lyapunov iterations. The proposed algorithm has the property of the quadratic convergence. It shows that the Newton’s method can be used well to solve the GCMARE under the appropriate initial conditions. Then the low-gain sample-data controller with a high design performance is realized though the digital redesign method. Finally, in order to improve the tracking performance of the proposed methodology, the matrix sign function and the chaos evolutionary programming approach are proposed to turn the forward gain Ec . Simulation results are good and efficient to tracking performance.

3-1

Reference

[1] P. V Kokotovi´c, H. K. Khalil, and J. O’Reilly, Singular perturbation methods in control: Analysis and design. Philadelphia: SIAM. 1999. [2] Z. Gajic´, and M-. T. Lim, Optimal control of singularly perturbed linear systems and applications. New York: Marcel Dekker. 2001. [3] T.-Y. Li and Z. Gajic´, “Lyapunov iterations for solving coupled algebraic lyapunov equations of Nash differential games and algebraic Riccati equations of zero-sum games,” in New Trends in Dynamic Games and Applications. Boston, MA: Birkhauser, pp. 333–351, 1994 [4] J. R. Magnus and H. Neudecker, Matrix Differential Calculus with Applications in Statistics and Econometrics. New York: Wiley, 1999. [5] A. W. Starr and Y. C. Ho, “Nonzero-sum differential games,” J. Optim. Theory Appl., vol. 3, no. 3, pp. 184–206, 1969. [6] D. J. N. Limebeer, B. D. O. Anderson, and B. Hendel, “A Nash game approach to mixed H 2 / H ∞ control,” IEEE Trans. Autom. Contr., vol. 39, no. 1, pp. 69–82, Jan. 1994. [7] H. K. Khalil and P. V. Kokotoviæ, “Feedback and well-posedness of singularly perturbed Nash games,” IEEE Trans. Autom. Contr., vol. 24, no. 5, pp. 699–708, Oct. 1979. [8] H. K. Khalil and P. V. Kokotoviæ, “Control of linear systems with multiparameter singular perturbations,” Automatica, vol. 15, no. 2, pp. 197–207, 1979. [9] Z. Gajic´, D. Petkovski, and X. Shen, Singularly perturbed and weakly coupled linear system − a recursive approach. Lecture Notes in Control and Information Sciences. Berlin:Springer. 1990.

R-1

[10] H. Mukaidani, H. Xu, and K. Mizukami, “Recursive approach of H ∞ control problems. for singularly perturbed systems under perfect and imperfect state measurements.” International Journal of systems Science, vol. 30, 467–477, 1999. [11] H. Mukaidani, H. Xu, and K. Mizukami, “Numerical algorithm for solving cross-coupled algebraic Riccati equations of singularly perturbed systems,” Advances in Dynamic Games: Applications to Economics, Finance, Optimization and Stochastic Control, A. S. Nowak and K. Szajowski, Eds. Boston: Birkhauser, pp. 545–570, 2005. [12] H. Mukaidani, T. Shimomura, and H. Xu, “Numerical algorithm for solving cross-coupled algebraic Riccati equations related to Nash games of multimodeling systems,” in Proc. IEEE Conf. Decision and Control, pp. 4167–4172, 2002. [13] Z. Gajic, D. Petkovski, and X. Shen, “Singularly Perturbed and Weakly Coupled Linear System—A Recursive Approach,” Lecture Notes in Control and Information Sciences

Berlin, Germany: Springer-Verlag, vol. 140, 1990.

[14] H. Mukaidani, “A new design approach for solving linear quadratic Nash games of multiparameter singularly perturbed systems,” IEEE Trans. Circuits Syst. I, Regu. papers, vol. 52, no. 5, pp. 960-974, May 2005. [15] F. L.Lewis and V. L. Syrmos, Optimal Control, 2nd Edition, New York: Wiley, 1986 [16] S. M. Guo, L. S. Shieh, G.. Chen and C. F. Lin, “Effective chaotic orbit tracker: a prediction-based digital redesign approach,” IEEE Transaction on Circuits and Systems-I, Fundamental Theory and Application vol. 47 no.11, 1557-1570, 2000. [17] G. Jank and G. Kun et al., “Solutions of generalized Riccati differential equations and their approximation,” in Computational Methods and Function, Theory , N. Papamichael et al., Eds. Singapore: World Scientific, pp. 1–18, 1998. [18] K. Lu, J. H. Sun, R. B. Ou, and L. Y. Huang., Chaos dynamics, Shanghan: shanghai Translation Press Company, 1990

R-2

[19] J. H. Halton, “On the efficiency of certain quasi-random sequences of points in evaluating multidimensional integrals”, Numerische Mathematik, vol. 2, pp. 84-90, 1960 with Corrigenda on pp. 196. [20] T. Yamamoto, “A method for finding sharp error bounds for Newton’s method under the Kantorvich assumptions,” Num. Math., vol. 49, pp. 203–220, 1986. [21] J.C. Van Der Corput, “Verteilungsfunktionen,” Proc. Kon Akad. Wet., Amsterdam, vol. 38,

pp. A12-A21, pp. 1058-1066, 1935.

[22] S. M. Guo, K. T. Liu, J. S. H. Tsai, and L. S. Shieh, “An observer-based tracker for hybrid interval chaotic systems with saturating inputs: the chaos evolutionary programming approach,” Computers and mathematics with Applications (Master Thesis, Nation Cheng Kung University) [23] J. S. H. Tsai, L. S. Shieh and R. E. Yates, “Fast and stable algorithms for computing the principal nth root of a complex matrix and the matrix sector function,” Comput. Math. Application, vol. 15, no.11, pp. 903-913, 1988.

R-3

Appendix A. The principal nth root of a matrix and the associated matrix sector function22, 23 Definition A.1. Let a matrix A ∈ C m×m , eigenspectrum σ ( A) = {λi , i = 1,2, ", m} , eigenvalues

λi ≠ 0 and arg(λi ) ≠ π . The principal nth root of A defined as a positive integer and (a)

( A)

n

n

n

A ∈ C m×m , where n is

= A , (b) arg(σ ( n A )) ∈ ( −π / n, π / n) . The principal nth

root of a matrix is unique.

Definition A.2. Let

a

matrix

A ∈ C m× m

σ ( A) = {λi , i = 1,2,", m}

,

λi ≠ 0

,

and

arg(λi ) ≠ 2π (k + 1 / 2) / n for k ∈ [0, n − 1] . In addition, let M be a modal matrix of A , i.e. A = MJ A M −1 , where J A is a matrix containing Jordan blocks of A . Then, the matrix sector function of A , denoted by Sectorn ( A) or S n ( A) , is defined as " 0 0 ⎤ ⎡ S n (λ1 ) ⎢ 0 0 ⎥⎥ −1 S n (λ 2 ) M , Sectorn ( A) = S n ( A) = M ⎢ ⎢ # % # ⎥ ⎥ ⎢ 0 " S n (λ m ) ⎦ ⎣ 0

where

S n (λ i )

is the scalar n-sector function of

S n (λi ) = e j 2 kπ / n = λi

/

n

(A.1)

λi , which is defined as

λi n with k ∈ [0,1] where λi lies in the interior of the minor

sector in C bounded by the sector angles 2π ( k − 1 / 2) / n and 2π ( k + 1 / 2) / n , and n

λi n is the principal nth root of λi .When n = 2 , the scalar n-sector function of λi n

becomes

the

scalar

sign

Sector2 (λ i ) = S 2 (λi ) = e jkπ = λ i

function

/

2

λi ,

of

denoted

by

Sign(λi ) ,

i.e.

λi 2 ≡ Sign(λi ) with k ∈ [0,1] .

The matrix sector function S n ( A) defined in Definition A.2. can be expressed as

( )

S n ( A) = A n A n

−1

A-1

(A.2)

where

n

A n is the principal nth root of A n . Also, the associated matrix sign function,

denoted by Sign ( A) , becomes

( )

Sign( A) = S 2 ( A) = A 2 A 2

−1

(A.3)

The partitioned matrix sector function of A can be described as follows :

Definition A.3. Let

a

matrix

A ∈ C m× m

,

σ ( A) = {λi , i = 1,2,", m}

λi ≠ 0

,

and

arg(λi ) ≠ 2π ( p + 1 / 2) / n for p ∈ [0, n − 1] . Also, let M be a modal matrix of A . Then, the qth matrix n-sector function of A , denoted by S n( q ) ( A) , is defined as ⎤ ⎡ S n( q ) (λ1 ) " 0 0 ⎥ ⎢ (q) S n (λ 2 ) 0 0 ⎥ M −1 , S n( q ) ( A) = M ⎢ ⎥ ⎢ # % # ⎥ ⎢ 0 " S n( q ) (λ m )⎦⎥ ⎣⎢ 0

(A.4)

where the qth scalar n-sector function of λi , denoted by S n( q ) (λ i ) , is ⎧1, 2π (q − 1 / 2) / n < arg(λi ) < 2π (q + 1 / 2) / n for q ∈ [0, n − 1] ⎫ S n( q ) (λi ) = ⎨ ⎬ . ⎩0, otherwise ⎭ The qth matrix n-sector function of A can be obtained by the following equation

[

]

i −1 1 n S n ( A)e − j 2 qπ / n for q ∈ [0, n − 1] . (A.5) ∑ n i =1 When n = 2 and q = 1 , the qth matrix n-sector function of A becomes the matrix

S n( q ) ( A) =

sign minus function of A , denoted by Sign − ( A) or S 2(1) ( A) as Sign − ( A) = S 2(1) ( A) =

1 [I m − Sign( A)]. 2

(A.6a)

The matrix sign plus function of A is Sign + ( A) = I m − Sign − ( A) .

(A.6b)

A generalized fast and stable algorithm with rth-order convergence rate for computing the principal nth root of a given matrix A is listed below: When the order of the desired convergence rate r = 2 ,

{

G (l + 1) = G (l ) [2 I m + (n − 2)G (l )][I m + (n − 1)G (l )]

},

−1 n

A-2

(A.7a)

G (0) = A , lim G (l ) = I m ; l →∞

R(l + 1) = R(l )[2I m + (n − 2)G(l )] [I m + (n − 1)G(l )] , −1

(A.7b)

R(0) = I m , lim R (l ) = n A . l →∞

When the order of the desired convergence rate r = 3 , ⎧⎪⎡ ⎛ n 2 + 5n − 12 ⎞ ⎛ n 2 − 5n + 6 ⎞ 2 ⎤ ⎜ ⎟ ⎟⎟G (l )⎥ G (l + 1) = G (l )⎨⎢3I m + ⎜ ⎟G (l ) + ⎜⎜ 2 2 ⎪⎩⎣ ⎝ ⎠ ⎝ ⎠ ⎦

⎫⎪ ⎬ ⎪⎭

n

−1 ⎧⎪ ⎡ ⎛ n 2 − 3n + 2 ⎞ 2 ⎤ ⎫⎪ ⎛ n 2 + 3n − 4 ⎞ ⎟⎟G (l )⎥ ⎬ , ⎟⎟G (l ) + ⎜⎜ ⎨ × ⎢ I m + ⎜⎜ 2 2 ⎝ ⎠ ⎝ ⎠ ⎪⎩ ⎣ ⎦ ⎪⎭

(A.8a)

G (0) = A , lim G (l ) = I m ; l →∞

⎡ ⎛ n 2 − 5n + 6 ⎞ 2 ⎤ ⎛ n 2 + 5n − 12 ⎞ ⎟⎟G (l )⎥ ⎟⎟G (l ) + ⎜⎜ R (l + 1) = R (l ) ⎢3I m + ⎜⎜ 2 2 ⎝ ⎠ ⎝ ⎠ ⎣ ⎦ ⎡ ⎛ n 2 + 3n − 4 ⎞ ⎛ n 2 − 3n + 2 ⎞ 2 ⎤ ⎟⎟G (l ) + ⎜⎜ ⎟⎟G (l )⎥ , × ⎢ I m + ⎜⎜ 2 2 ⎝ ⎠ ⎝ ⎠ ⎣ ⎦

−1

(A.8b)

R(0) = I m , lim R (l ) = n A . l →∞

Also, a generalized fast and stable algorithm with rth-order convergence rate for computing the matrix sector function of a given matrix A is as follows : When the order of the desired convergence rate r = 2 ,

[

][

Q (l + 1) = Q (l ) 2 I m + (n − 2)Q n (l ) I m + (n − 1)Q n (l )

]

−1

,

(A.9)

Q (0) = A , lim Q (l ) = S n ( A) . l →∞

When the order of the desired convergence rate r = 3 ,

⎡ n 2 + 5n − 12 n n 2 − 5n + 6 2 n ⎤ Q(l + 1) = Q(l ) ⎢3I m + Q (l ) + Q (l )⎥ 2 2 ⎣ ⎦ −1

⎡ n 2 + 3n − 4 n n 2 − 3n + 2 2 n ⎤ × ⎢I m + Q (l ) + Q (l )⎥ , 2 2 ⎦ ⎣

A-3

(A.10)

Q (0) = A , lim Q (l ) = S n ( A) . l →∞

A-4

誌謝

感謝恩師蔡聖鴻教授在我碩士班求學過程中，在研究及課業上的指導與教誨，建 立正確的研究觀念與態度，讓學生受益匪淺。感謝張郁斌教授與陳嘉偉教授對本論文 的不吝指正與熱心建議，使本論文內容更趨完整。 最後感謝多變數實驗室的成員對於我在研究上及生活上的幫助及關懷，以及周邊 朋友的加油打氣，在此致上由衷的謝意。有了大家的支持，我才能順利完成此畢業論 文

在此致上最誠摯的謝意。

View publication stats