An Algorithm for Factoring Hermitian Matrix Polynomials Relative to the Imaginary Axis

2a-08 6

Copyright © 1996 IFAC 13th Triennial World Congress, San Francisco, USA

An Algorithm for Factoring Hermitian Matrix Polynomials Relative to the Imaginary Axis Chyi Hwang and Bo-Win Lin DCl'u,r·tmcnt of Chemical Engineering iVa.hona[ Chnng Chcnq University Chia- Yi 621 TA TW4N

and Tong-Yi Guo lJepartment of Chemical Engineering Na.tional Kaohsiunq Institute of Technology Kaoshiung 807 TAIW4N

Abstract: In this paper we apply Cauchy's residue theory to convert the problem of factoring a Hermitian matrix polynomial relative to the imaginary axiH into tha.t of solving et block- Toeplitz system of linear equations in which the coefficients are indefinite integrals of rational functions. The involved indefinite integrals can be transformed into definite ones by the variable transformation 8 = jtan£, which can then be ac.curately and effectively evaluated by a numerical integration scheme capable of acc.urac.v control. Since a block- Toeplitz system of linear equations can be efficiently solved by it fast algorithm, the proposed "Igorit.hrn is nOllit.erative in nature and therefore doc~ not have t.ll(~ problem of convergence and inHial guess assoc.iated \l,.'ith all iterative a.pproach, Keywords: Polynomials Factorization methods, Linear systems, Numerical methods, Wiener filters, LQG control 1

1. INTRODUCTION

The problem of factoriTl!'; a Hermitian matrix poly· nomial relative to the imaginary axis has a close cOIlnection with t.he studies of optimal LQG and JI= controller synthesis (nOIl!,;iorno, 1969; Jonckheere and SilvNmaIl, 1981: JOllckheere and Verma, 1986; K ucera, 1979; Park and Youia, 1992; Youla ct aI., 197G), optimal prediction and filtering (Roberts and "Iewmann, 1988; Wiener, 1949; Wong and Thomas, 19(1), slochast.ic realization th,,ory (Caines, 1988), and network synthesis (Anderson and Vongpanitlerd. 197:l). Several authors (Aliyev

ct al., 1990a, 1990b, 1991; Aliyev and Larin, 1993; Anderson et al., 1974; Anderson and Vongpanitlerd, 1973; Callicr, 1985; Davis, 1963; Jezek and Kucera, 1985; Larin, 1992; Rissen, 1973; Tuel, 1968; Wilson, 1972) have devoted their effort to develop algorithms for the factorization of Hcrmitia.n matrix polynomials. Except for the algorithms due to Davis (1963) and Callicr (1%5), most factoriz"tion algorithms that have been proposed up to now are iterative ones. The factorization proc:edures are based on solving it matrix algebraic Riccati equation (Aliycv et al., 1990a, 1990b, 1991; Anderson et al., 1974; Tuel, 1968) or a suc:cessive of Lyapunov equations

1263

(Aliyev a.nd Larin, 1993; Larin, 1992), on performing a triangular fartorization of Toeplitz matriees (Rissen, 1973), or on the Newton-Raphson method (Jewk awl Kacera, 1985; Wilson, (972). The construc.tion of solution to an algebraic matrix Riccati equation or Lyapullov equation is llsnally aCCOIIlplished by using the matrix sign functiOIl approach (Roberts, 1980) or the Newt.oJl-Raphson approach (AJlderson, 1978). It is well-known hot.h the matrix sign function and NewtoTl- Raphson approach are iterative algorithms tnat have the problems of convergence and initial gllesses. On the other hand, the recursive nat.ure of the algorithms hinders them III adapting to a paraUcI computing pTlvironment,

Although the algorithms proposed by Davis (1963) and Calli"r (1985) fi)r factoring a matrix polynomial are finite-step procE~dures~ they requ.ire the zeros of a high order polynomial, whirh is usually a numerically ill-conditioned computation. In this pap"r, we extend the polynomial factorization approach proposed by Arbenz (1979) and Hwang et al. (1995) to factoriy;e a. Hermitjan nlatrix polynomial relative to the imaginary axis. The algorithm presented here is also based on using Cauchy's residue theory to convert the problem of matrix polynomial factorization along the imaginary axis into t.hat of eva.1uating a set of indefinite integrals and tha.l of solving a block-To~plitz system or linear equations. It is note that a indefinite intpgral along t.he imaginary axis can he transformed into a definite one by the variable transforrnaticJT] ,<; = tan~ (Hwang and Chuang, 1994), and that a definite integral ran be partitioned into several ones \\li1,h a smaller integration interval. Consequently, the values of resultant indefinite integrals can be obtained accurately with a numerical integration scheme capable of accuracy control, and the computation time eat) be reduced in a parallel computing environment. This parallel computation capability of ev,duating integrals along with the availability of fast algorithms (Akaike, 1975; Glentis and Kalouptsidis, 199·1: .Iou ct al., 1986) for a block-Toeplitz system of linear equations makes the proposed algorithm of matrix polynomial factorization attractive to real-time applications.

Given a regular m

m lIermitian matrix polynomial

X

pes) = (-lrIms2n+P2n_182n-l+ ... +P,S+Po (1) with P(s) - pT(_s) and P(jw) > 0 (positive definite) for all real UI~ the problem considered is to determine an m x m matrix polynomial G( s) satisfying the relations hi p

(2) under the conditioTl that G- (.s) has no poles in the right half-plane. Here and in the sequel, Im is the rn X rn j(h~lltity matrix, the s11persr.ript T denotes transpose. 1

The matrix polynomial G(s) that satisfies (2) can be written in the form

G(.-)

= I",8 + G , s,,-l + ... + G n - 1", + G n n

(3)

where the Tn >< rn coefficient matrices Go, G 1 ,···, G n - 1 are to be determined under the aSSll rnption that G - I (.,,) has no poles inside the right half-plane, i.e., G- T ( -.s) = (G( _s)T)-1 is analytic inside the closed left ha.lf-plane. Let

{jw; {Re

-00

j8 . IT

.

<

W

< oo} a7r

U

2 <:: (I <:: "2' R

-->

oo}

(4)

be the closed path that encloses the entire left halfplane. Then, it follows rrom Callchy's residue theorem (Conway, 1978) that the contour integrals of "kG-'1"( _,,), k = 0,1, ... , n along the curve C vanish, i,e.

] ,,-

+ A~=JE' (Rc,O)kG-T( -Re,o)d(Re j o) , o (5) Since skG- T (_8) can be expanded as nH ,kG- T ( -.,) = Ims-n+k +

afs-

1

+...

(6)

the second integral on the right-hand side of (5) is evaluated to be

2, PROBLEM STATEMENT AND FORMULATION

k = 0, l. ... n - 2 k=n-I

1264

(7)

A combination of (c») and (7) leads to

+.,=

k

Q(c')

T

. -.1,:X) s G-- (-,<;)d,<;

/

{

O,

k=O,I. ... ,n-2

lflml

k-:::;::n-(

(8)

Now subst.itut.ing J

G(S)p-l(S) n-l

.,np-l(S)

+L

GI-<,p-l(.,) (9)

= (-I)"Q2ns2n + Q2n_l,,2n-l + ... + Q1S + Qo

(13) with the matrix Q2n being singular can be converted into the regular one by constructing a matrix polynomial T(.s) such that TTr-s)Q(s)T(s) =' P(.,), where Pis) is given in the form of (1). After finding the factor G(s) for pes), the factor H(s) that satisfies HT( - 8 )H( s) = Q( s) is then obtained as H(s) = G(S)T-l (s). A procedure for constructing such a transformation matrix T(s) is described in (Aliycv et ai., 1990a, 1991) .

/::::0

3. COMPUTATIONAL ALGORITHMS

into (8), we obtain

r:

In view of the foregoing derivations, it is obvious that the presentpd approach to the factorization of Hers"+'P-'(s)ds+ G, s'+'P-'(s)ds mitian mat.rix polynomials with respect to the imaginary axis depends heavily on the evaluation of ink=(I,l, ... ,n-2 0, r10) { 1!"I"" k = 11 - 1 definite integrals VI in (I]). As indicated in (H wang and Chuang, 1994), the indefinite integral along the Letting imaginary axis can be converted into a definite onc 1= { -1 by the varia.ble transformation s = jtan(B/2). With 1=0,1, .. ·,2n-l VI = . 8 P (.s)ds, this varia.ble change, the indefinite integrals Vl be,/ -J= (11 ) come and noting = we put the equations in (10) e B B dB into the following; matrix fann: VI = 2 Re{(jtan(,,))'P-' (jtan(....-)}-.----;; o L £ I + cosO (14 ) V n- 1 Vn-:l VD Since the accuracy of the matrix polynomial factor Vu V",_l VI G(s) depends on Vl, it is essential to ensure the accuracy of evaluating the above definite integrals. An effective approach to accurately evaluating dPiiuite integrals (1.1) is to solve the initial-value problems V V 2n-2 V 2n-:~

l,:,

~

vi v"

l

(:D

n

-V"

)

-V n +1 (

- V 21,-l +-

l

(12 ) KIm

This is a block- Toeplitz system of linear equations for the unknown coefficient matrices of the solution G(s). Hence, we have transformed the problem of factoring a Hermitian matrix polynomial relative to the imaginary axis to that of pvaluating the indefinite integrals V, in (11) and solving the blockToeplitz system of equations in (12). It is noted in the above deriva.tion. wc have assumed that the mat.rix polynomial P( 8) is regular. However a general matrix polynomial of form

e

from = 0 to 8 = 11", with " numerical integration scheme capable of accuracy control. With this dif~ ferential equation approach, the values of integrals VI are given by VI = V I(,,). Once the matdx cntdes VI, I = 0: 1,·· ., 2n - 1 of the block-Toeplitz system of linear equations (12) are accurately evaluated, one can apply one of the available fast block-'l'oeplitz solvers (Akaike, 1975; Glentis and Kaloll ptsidis, 11194) to find the. coefficient matrices Gi, i = 1,2"", n. It is noted that in a parallel computing environment, the computation t.ime for eva,luating the definit.e integrals Vl can he reduced if each of these integrals is partitioned into

1265

At, i.e.,

Besides, the computation efficiency can be enhanced

Table 1 shows the obtained coefficient matrices G , and G 2 of G( s) for various values of tolerance parameter E. l'he number of funct.ion calls, Nj,k, performed in solving t.he set of differential equatioIls are also included. It can be verified that the e"ch obtaiTJ(,d G(,,) in t.he last. t.hree rows of Table 1 satisfies G-1(-s)G(s) = P(s) and G- 1 (s) has all its poles in the Jpft. half-plane.

if the explicit expression for inverse of the matrix polynomia.l P(s) is determined a priori by the the

5. CONCLlISIONS

114 ones \Vitll a smaller interval

~"'I

f

7r /

rr

()dB

(16)

. (A1-1)7r/A1

algorithm given by lnouye (Inouve, 1979).

4. AN ILLUSTRATIVE EXAMPLE Consider the Hermitian Hlfll.rix polynomial Pis) =

9]

+[~

C,

l7

It is desired to find the matrix polynomial G( s)

sllch that P(.s) = G"/"(-S)G(8) and G-'(s) is analytic in tl", closed right half-plane. The inverse of P( s) is given hy

We have presented a noniterative method for the factorization of Hermitian Inatrix polynomials \vith respect t.o t.he imaginary axis. The method is based on using the Ca.u('hy~s residue theory to transform the factorization problem to a block- Toeplitz system of linear equations in which the coefficients are integrals of rational functions
evaluate the the indefinite integrals. Since a definite integral GHl be partitioTlPd into several ones with a smaller integration interval, and the block- Toepllt7. system of equations can he efficiently solved by a rast algorithm, the proposed method is pa,rt.iClJlarly at.tractiv(~ 1.0 h(-, adapted 1.0 a. parallel computing environment.

p-'(.,)

ACKNOWLEDGEMENT

[~ ~]s4+[~~a s~

- 12s 6

'"(,) ,,,(,I I [ 1'2ds) 1'22(8)

+ 408 4

-

3:~s2

+4

This work was support.ed h.v the National Sf'.ienf'.e COllTlcil of Repllblic: of Chill,. under Grants KSC842214-E-194-001 and KSC85-nI4-E-151-UOI.

For computing VI = [VI.],k], I = 0,1,2,3, we solv(' the following 4 set.s of different.ial equations dV{,j.k

dB

2Re{(jtan(

~ ))lrJ,k(jta.n(~))} 1 +~:osO'

1= 0, 1,2,3 tor j,k = 1,2 by invoking t.he subroutine DIVPAG of the popular nrSL mathematical library (IMSL, 1989). In calling DIVPAG, the initial step length is set to 10-:3 while all other parameters are set to their default values. The accuraev of integration is controlled by the tolerance parameter E. The smaller value of E: is scL tIll' more accurate int.egrals values are obtained\ and consequently the more llUlnber of function calls are p(:·rformed.

REFERENCES Akaike, H. (1!J75). SlOll.

Block Toeplitz matrix inver-

HIMv! J. Appl. Math., 24,234-241.

Aliy"v, F.A., B.A. BonlYllg and V.B. Larin (19903). Factorization of polynomial matrices respect to the ima.ginary and t.he unit circle. Soviet J. of Automation and Information Sciences, 20(4), .il-59 Aliyev, 1o'.A., B.A. Bordyug and V.B. Larin (1990h). All algorithm for the Factorization of polynomial matrices. Soviet Math. Vokl., 40, 1:37-140. Aliyev, F.A., II.A. Bordyug and V.B. Larin (1991). Fact.orization of polynomial matrices

1266

and Hcparation of ra,tional matrices. Soviet J. COTTllmler and S'ystems S'cit;ncc.s~ 28(2),17-58. Aliyev, F.A. a.nd v.n. Larin (1993). Generalized Lyapunov eCjlJ<'lt.ion and factorization of matrix polynomiah. ,'ystem8 1'1 Control Letters, 21, 485-49J. Anderson, B.U.a. (1978). Second order convergent algorithms for t.he st.ea(ly-state Riccati equation. 1nl. j. Control, 28, 295-306. Anderson, n.D.O., K.L. Hitz and N.D. Diem (1974). Recursive algorithm for spectral factorization, {FP}} Tmns, Circuit,., and Systems, 21, 742-750. Anderson, I3.Jl.(). ami S. Vongpanitlerd (197,1). Network Analysis and Synthr-sis -- A lV[ordf:T'Tl Systems 'l'heory ilpproQ(:h, r~Tlglewood CliffH, N ..J.: Prentice Hall. Arbenz, K. (1979). Spektrale faktorzerlegung eines poly noms mi! hilfe de LaplaceTransformatioll. Arch. Elf'k. Ubertr., 33,42142:t Bongiorno, J.J. (1969). \1inirnurn sensitivity design of linear rnultivariable feedback control systems by matrix factorizat.ion, IEEE Trans, Automat. Cont1'., 1-1. 665-Gn. Caines, r.E. (19~8). ['i"ell'r Stoehastic Processes, Wiley, New York. Callier, F.M. (HJH5). On polynomial matrix spectral fadorization by sym metric extradion. JEEF Tmns. iI "tomat. Contr., 30, 453-464. Conway, .LB. (197H). Functions of One Complex, 2nd Ed .. S)lrin~er- Veriag, New York. Davis, M.C. (19G3). On faci.OT'ing the spe<:t.ral matrix. IFFF Tr(}.U8. it ulo'flwt. C~ont1'.) 8,296:l05. Glentis, G.O. and N. Kaloup1.sidis (1994). Efficient solution of block linear systems with Toeplitz entdeh using a channel decomposition technique. Signal Frocc8"ing, 37, 15-60. Tlwang, C. and Y.H. Chuang (1994). Computation of optimal reduced-order models with time delay. Chem. Eng. Sri., 49, :1219-3296. Hwang, C., .l.H. llwang, n.w. tin and T.Y. Guo (199S). Polynomial factorizaJion via Holving differential equations and Toeplitz system of equations. subrnited to Control a.nd Computers, [MSL (1989). Use,-'" AIrman! for iMSL Mathl foibrw'Y: Fortran Subroutines for Math-

emaiieal A pplical,ions, Version 1.1. In01ryc, Y. (1979). An ;tlgorithm for inverting polynomial matrices. 1nl. J. C()ntrol~ 30 , 989999. Jezek, J. and V. Kueera (1985). Efficient algorithm for ma.t.rix spcrtral f'l.ctorization. A utomat£ca, 21, 663-669. Jonckheere, L.A. and L.).1. SilverrnaIl (1981). Spectral theory of the linear-quadratic optimal control problem: analyt ical factorization of rational matrix-valued functions. SIAM.1. Control and Optimization, 19, 262-281. Jonckheere, E.A. and lvI. Vcrrna (1986). A spectral characterization of H=-optimal feedback pel'fonna,Tlce a.nd it.s dficient computatio~. Systems {1 Cont'rol Letters; 8, 1:3-2'2. Jou, I.C., Y.U. Ru and \V.S. Feng (1986). Novel implementa.t.ion of pipelined Tocplitz system solver. ProI'. IEEE, 74. 1463-1464. Kucera, V. (1979). iJiscrete Linear Con/ml: The Polynomial iI pproach, Wiley, Chichester. Larin, v.n. (1992). TI,e generalized Lyapunov equat.ion aml fftctorization of matrix polynomia.1s, J. of Automation awl Information Sciences, 25(6), 16. Park, K. and D.C. Youla (1992). Numerical calculation of the optimal Ihree-degree-of-freedom Wiener-Hopf controller;. Int. J. Control, 56, 227-244. Rissen, ,I. (1973). Algorithm for triangular decomposition of block Hankel and Toeplitz matrices \vith applkation to factoring positive matrix polynomials. Mathematics of Computation, 27,

117-15,1. Roberts, .I.D. (1980). Lineal' model reduction and solution of the algebraic Hjccati equation by the use of th{~ ~ign function. 1nl. J. Control, 32, 677-678. Roberts, A.P., and M.M. Newmann (1988). Polynomial approach to Wiener filtering. into J. Contml, 47, 681-696. Tuel, W,G. (1968). Computer algorithm for spectral factorization. IBM.1, Res, Dev., 12, 163170. WieJ"'I', N. (1949). t;xtr'opoiation, interpolation and Smoothing 0] Stationa.ry Time Series, \ViIcy, New York. Wilson, G.T. (1!l72). The faetorization of mat ri-

1267

cial spectral densities. SIA!If J. Appl. Math. , 23, 420-426. Wong, E. and fc. Thomas (1(J(iI). On the multidimensional predictioTl alld fil t. f.:! ring problem and the factori zation of spectral matrices. J. Frank/in [nsl. , 272 , 87-!l9 .

Table f

10- 3 1O - ~

10- 5 10- 6 10- 7

10- " 10- 9

G, -1.184331 2.299154 1.19662 1 2.599120 2.295817 1.179272 -1.178139 3.710149 2.294221 -1.176698 1. 17697 1 3.706694 2. 294127 1.176500 -1.l76515 3.705096 2.294117 - 1.l 76471 1.176470 3.705886 2,294 11 8 1.176471 3.705882 - 1.l76470 1.1 764 71 2.204118 :3.705882 1.176471

Youla, f) .C. and J .J . l:Iongiorno (1985). A fecdback theo ry of two- d~gTee-of-freedom optimal Wiener-Hopf design. IPPE Tram. Automat. Cont /'., 30, 652- 6G5 . Youla, D.C. , H. Jabr an (I J.J. ilongiorno (1976). Modern Wiener-Hopf design of optimal controllers. Part 11: T he mult ivariablc case. [ E1:)1:) Tmns. Automat. Contr. , 21 , 319-338.

1

Gz 1.730327 2.599120 1.830844 2.777825 1.822133 2.762096 1.823432 2.764535 1.823518 2.764686 1.823.530 2.761707 1.823539 2.764706

1208720 2,908331 1.298185 3.067142 1.292465 3,055851 1.293968 3. 058559 1.29410" 3.058801 1.294118 3.058825 - 1.294118 3.058824

1268

lVit

N 1Z

N21

Nn

52

49

49

49

71

93

84

84

100

116

104

104

106

152

129

129

143

151

160

160

182

202

186

186

226

242

245

24 3