- Email: [email protected]
Robustness Analysis Using the Mapping Theorem without Frequency Gridding
2d-09 2
Copyright © 1996 IFAC 13th Triennial World Cong ress. San Francisco. USA
REM ROBU STNE SS ANAL YSIS USING THE MAPP ING THEO DING GRID Y WITH OUT FREQ UENC Gilles Ferrer es " J ean-Fran~ois Magni ., 1 Edouarri • Automatic Control Departm ent (DERA ), CERT- ONERA , 2 Avenue France Cedex, e B elin, BP 4025, F31055 Toulous
comput in g a. lower Abstra ct. This pa per proposes a new numerical techniq ue for The te(~hnique is g. griddin cy frequen bound of th e robust stabilit y margin , withou t t hp. numeric al that way a such in ed organiz is based on the Ma pping Theore m and unlikc in edges of number the of instead vertices of comple xity is related to the number a ban k ring conside trated demons is cy efficien al numeric The o ther known techniq ues. le for tradab he to of unce rtai n system s. T he stabilit y margin comput a tion is shown damped }' lightl several h t more than twelve uncerta in parame ters even for system s wi modes. tric varia.ti ons K e yword s. Sta bility margin , stabilit y test, robustDess, parame 1. INTRO DUCTI ON
Checkin g the sta bility of a family of polynom ials pes, 6), wh ere the vector 8 of Illlcerta in parame ters belongs t o an hype rrectan gle D, is a classical problem in ro bustnes s analysis. Exact solut ions a re availa ble when the coefficients of P(s, 6) belon g to interval s or arc affine functions of 6 (n amely the Kharit onov and Edge theorem s sec (Dasgu pta , 1988), (Bartle t! et al. 1987)). The problem is howeve r mu ch mo re difficult when the coefficie nts of P(s, c5 ) a re multilin e;u- or polynom ial functio ns of b. In the context of a multilin ear depend ency, an ou t.erboundi ng set of the family of polyno mials P{ s, h) ~ na mely the polytop e conv (P(s , D ;)) where D ; is a vertex of the hyperre ctanglc D, can be conside red. It ~uffices t hen to check the stabili ty of the (expose d) edges so aB to check t he stabili ty of t he poly tape of polynomia.ls. There is howeve r a combin atoria l explosi on with the size /11 1 of the problem , because of the 2 - (2'" -1 ) edges of the t.un para meters). An uncer polytop e (m. is the number of in (Bouguerra et al, cl propo~e heuristi c method has been burden: conside rtional a comput the 1990) for redu cing is perfor med, test logical a e, polytop he t of ing an edge for the staion t condi nt sufficie a to onds corresp ch whi hility of th e edge, further comput ation is thus necessa ry only if t his logical test is not verified . It can be howeve r noted that t he compn ta tio nal comple xity is still more 7 11 or less related to 4HI ::::: 2/11 "'\ (2 - 1), since it is still necessa ry to con ~i cler a ll edges of {,he polytop e. Morees I Suppor ted by Di rection d es Rcd lert:hes et Etudes Te<::hniqu
over, the comput ational require ment for verifyin g this logical t est may reveal non negligib le, depend in g on the im plemen tation. An a priori attracti ve solution is to combin e a frequen cy griddi ng wit h the zero exclusio n principle. Let conside r the con vex hull of t he points P(jw, D;) in the com plex j w (remem ber D j is one the 2Tf1 vertices of plane at s the hyperre clangle D ). The poly to pe of poly nomia ls is sta hle if a nd only if one elemen t of the polytop c is stable a.nd if the convex hull does not contain the origin of th e comple x plane a t all frequencies w. Noting N the number of frequen cy points, the comput ational comple xi ty is now N 2"1 . Neverth eles!-S, even if a frequen cy griddin g is rather usual ill robustn ess a nalysis pro blems (especia lly in the context of I.£~a.naly sis) . t his techniq ue may reveal unreliable in the specific context of lightly da mp ed b e nd~ ing mo nes, since it appe.ars possibl e t o miss the critica l frequen cy in t his case (FreuIlfle.nberg and Morto n, 1992).
=
The key idea of this paper is t o directly manipu late the continu um of convex hulls corresp onding to w E [O, oo{ and t o check whethe r the origin of the co mplex plane belongs t o this continu um. As a ma tter of fact , this freque ncy a pproacl, gi ves U H further insig ht inside our problem , and it is especially possible to reinterp ret the suffi cient conditio n in (Bougu crra et ai, 1990)) in this new co ntext . A sufficien t condit ion for t he sta bility of th e polyt ope of polynom iah, is then proposc d: its (".omputatio nal comple xity is 2111 . To check the stabilit y of the Volytop e of polyno mials (or equival ently to check wh et her t he continu um of c.onvex hulls contain s the origin of t he comple x pla. ne), it is necessary in the worst
(DRET)
3404
case t.o apply 2H1 times this sufficient condition. Even if the worst case complexity is still 4"\ numerous examples treated in §4 tend to prove tha.t the mean complexity is closer to 2fl1. As an important consequence, it seems pot>Siblc to deal with problems involving up to 12 uncertain parameters with a recwonable computational burden (about 5 minutes for polynomia.ls with degree 20) instead of 7 parameters with a. hrute application of the Edge Theorem. The paper is organized as follows. In §2 are sketched existing methods which are not based 011 frequency griddingo §3 presents the proposed method for computing the lower bound of the robustness margin, while the examples of §4 emphasize the practical interest of the method.
In view of this elementary result, two existing t echniques are briefly ske tched. M e thod 1. In order to check t he stability of a polytope of polynomials, (Bialas, 1985) f:onsiders a test which enables to check the ,tability of each edge . Bialas test requires the computation of thi~ eigenvalues of a n x n matrix deriveu from Hurwitz matrices. Here is used an alternative t est (see Bouguerra. et al, see also (ChapeUat a.nd Battacharya , 1989)) which is an immediate COll f)P-qllence of LeUlma 1. Let us consider two vertices DJ: aJld Dj of the pa.rameter box D. We have to check whether the segm ents [0, P(jw , D d ] and [0, P(jw, D d ] are collinear at some (n'! al) frequency w, 01' similarly have the same slope: 'S(P(jw, D,)) 'S(P(jw, Dd) fR(P(jw, Dd) - 1R(P(jw, Dc))
2. PROBLEM STATEMENT, REVIEW OF EXISTl,' W METHODS
Thiti is equivalent to lo ok for positive real roots of th e w-polynomial
Let us consider a polynomial
!R(P(jw, DtlP(P(jw, Dd) -l?( P(jw, D, ))"(P (jw, Dtl) (3)
bED where D is an hyper rectangle of Rill i. c. a set of points defined by D = {6= (61 ,62 , ... ,6",) , 6,- ::;6i ::;6t}
6t .
for given bounds 6; and The coefficients ado), i = 1, ... , n are assumed to be multiltnwT' junctions of 01, ... , li,rl' The nominal poin t in the parameter box will be denoted Do, it corresponds to the vector 0.5((01 + lit) •... , (D~ + 6;t)). The 2'" vertices will be denoted Dj. The polynomial of (1) is generally given under the following {orm:
P(s,b) =det (s I-A - B diag(6" ... , hm )C) (2 ) where A E lR":<" , B E R" Xrt' and C E
JRJ/lXTI.
The value set at EL given point .~ is the domain of the complex plane defined as pes, D). The following result, which is a direct consequence of the Boundary Crossing and Mapping Theorems, will be used. Lemma 1. Let us a.c;sume that the nominal polynomial p es , Do) is stable. The system is stable if at the frequency w = 0 the value se t. P (O, D) does not contain the origin and if one never encounters, for all real values of w > 0 , thp. case where two verti ce~ of the convex hull and the origill are on the same straight line, with the origin in b etween the two vertice.').
(Note that if we consider this polynomial with respect to w2 its degree is equal to n -I , so the involved computation is of the same order than with the original test of Bialas.) If t he polynomial of Equation (3) admits a positive real root Wo . it remains to check whether the origin of the complex plane is b etwl~en the points P(jwo, D{) anu P(jwo. Dd. The answer to this qu estion is affirmative if and only if: ~(P(jwo , D,))fR(P(jwo,Dk)) < 0 or fR(P(jwo, D,))'R(P(jwo , D .)) 0 and "'(P (jwo , D,))'S(P(jwo, Dk)) ::; 0
=
(4)
It comes out that the complf·xity of Bia las test and of the altern ative test is similar: computation of the roots of 2(11 -1)(2" -1) polynomials of degree n is required in both cases. Method 2. A similar idea. is used in (Bouguerra et at, 1990), still considering 211- 1(21< -1) edges but the computat ion of the roots iD replaced in most cases by a logical test which is assumed to be faster. I nstead of considering the explanations given in (Bouguerra et al, 1990) let us giVf~ o ur point of view concerning this technique: in that way the proposed technique of §3 will appear (1.<; a 'tep furth er beyond Methods 1 and 2. Clearly edgeR of the convex hull, say the segm.mt [P(jw , Dd , P(jw , Dd] , v:m not touch the origin for w > 0 (so that Lemma 1 enables to conc.1 ucle) if these points are in a same half plane for all values of w. Thi~ propert.y is satisfied if for example
3405
!R(P(jw, Dd) and !R(P(jw, D.» ,ame sign \/w or
(S(P(jw, Dd) alld (S(P(jw, D.» ,ame sign \/w
• its 4oI:t _axis rl is collinear CO the oriented segment OMI_', the "y_axis" is such that the frame is direct. It is clear that the coordinates of the point M which were (X, Y) her.ome (X', V')
Therefore, Bouguerra td al propose to proceed as follows. If for some interval oC frequencies we have :
" _ XpX+YpV . Y'_ -YpX+XrV X X',.. +" y2 ' X 'p+p V'
!R(P( jw, D,)!R(P(jw , Dd) < 0 aDd (S(P(jw , DJ)\J(P(jw, Dd) < 0
As only the signs of X ' and Y' will be needed} we shall consider the following simplified coordin a tes:
(5)
X'=XrX+YtY; Y'=-Y,X+Xr Y (6) check conditions (3) a nd (4), otherwise the edge is known to be stable. The numerical complexity corresponds to t he computation of the positive real root~ of 2"+1 polynomials with degree d ose to n / 2 (th~e polynomials are !R(P(jw, D;» and U(P(jw , D;» , i = 1, . .. , 2"). From the knowledge of these roots , tlw logica l test of Equation (5) is considered for the 2"-'(2" - 1) edges. For edges such that (5) is satisfied in some intervals of frequencies, StnI'm Theorem enables to check the existence of real roots inside each suc.h interval. Sturm Theorem will be replaced here by the test based 011 Equations (3) and (4) . III the procedure proposed in §3 (Method 3), we shall try to avoid the logical t est concerning all edges because this test is not as fast as it might appear at first sight (depending on the implementation). We shall also improve the test of Equation (5) by considering a frame which turns with th e value set so that chcr:king the fact t ha.t P(jw , D,) and P(jw, D. ) belong to a same half plane becomes much more efficient. In fact instead of eliminating edges} wc eliminate intervals of frequencies for which we know that. the nmtinnum of convex hulls does not contain the origin. If we can eliminate all frequencies, it means that the system is robustly stable. Otherwise the first frequency at which th e (:ontinuum tou ches the origin is characterized.
• its "y-axis" is collinear to the oriented segment OM F , the "x-axis" is such that the frame is direct. It is dear that the coordina.tes of the p oint M become here (X' , V' ) (with simplification as ahovc):
X'=Y"X-X,.,Y; Y'=X,..X+YpY (7) First step. It is necessary to check that, at frequency
w _ 0, the value set pen, D) does not contain th e origin . Tt is well known th at P(O,D) = [ Min ;P (O, D, ) Maxi P(O , D;)I. SO, if the o rigin belongs to this set., P( .. , 6) is not robustly stable} otherwise go to next ste ps. Second ste p. In order to use Lemma 1 in an efficient way} it seems to be worth c( !O~idcring a frame which turns with the image of the lw minal parameter vector Do. This idea can be justified by considering the locus corresponding to the nominal ('ase. It is a point which remains (of course) on th e "xaaxis'l • When th e box around Do is no t too large, the valu e set will remain in the half plane X'( w) > 0 (see (6». So for s uch small boxes a sufficient c.ondition for stability is X'(w) > 0 for all positive fr~ucncies. Ivlore precisdy we have the following !)ufficicnt condition for sta bility:
= =
Lemma 2. Let XFCw) !R(P(jw , Do» and Vdw) If for i 1 ... 2" , the polynomials X,..(w) R(P(jw, D;» + Y" '(w) ~~(P(jw, D;» have 110 positive real root, the polynomial P( .• .I) is robustly stable.
~(P(jw, Do».
3. A NEW ROBUST STABILITY TEST Our main objective in this section is to derive a robust Rtability test based on value set idea.~ but without frequency gridding, also we intent to avoid numerical complexity related to the number of edges 2 11 - l (2" -1) (m denotes the number of uncertain parameters) by keeping it as close as possible to the number of vertices 21n. Change of coordinate in the complex plane. Before giving a description of the proposed stability test, we shall consider a change of variable. Let us consider a point M in the complex plane with affix X + JY. An other point M r with affix XF + .fYv i~ used to define a frame in the c:omplex plane £U:i follows:
So to check stability} if Lemma 2 is not !:iufficicnt to conclude, it means that the w 2 -polynomial X'(w) ~ xp(w)31(P(jw , D . JJ
+ YF(W)"'(P(jw, V i ))
(8)
admits some positive real roo ts, let WI be such that wf is the smallest of these roots. \Vhen w = w. th e image of one vertex touches the ll.y" -axis. Two pOSl;)ibilities are depicted in Figures 1 and 2. In the first case} we must find which vertex of D CO rTl ~8ponrls to the point (denoted. M 1 in Figure 1) with maximal slope in th e moving frame. More precisely WI ! look for the vertex Dil (Mc =P(jwJ,D;,)) such thaj SeD;,) defined by :
3406
-Yr(w,)'R(P(jw" Di, )) + XF(W, )'>(P(jW" Di, )) X r(w,)lI1(P(jw" Di,)) + Yr(w,)'>(P(jwJ, Di, ))
!'i
(9)
I
is maximal. At the next step , this point will be used to define the moving frame. In the second case (see Figure 2) the point 111, is ob tained directly. Additional steps. In the additional steps, the moving frame turns now with the image of Dil defined at the previous step. But DOW, it is its ""y'.a.xis" as in (7) which conta ins the image of Dil i.e. AIl P(jw , D i \) . As previously we are interested in the first frequency (,.1)2 at which the value set touches the new "y' -(tXis" . If such a frequency does not exist, the polynomial p es, a) is robustly stable. The algorit.hm ends.
=
Fig. 3. ith iicrntion: The ve1tex M i + 1 touches the axis on the negative side. therefon~ (t he boundary of) the convex hull contains the origin. So the algorithm ends . Generally we cannot conclude about instability except for some pairs of" vertices defining the segment [Ml . W2"1, see comment 3.2.
Second, the image of D j 2 is OIL the other side. In which case we define a new frame which will turn with M2 011 its Uyl -axis" i:LDd proceed t·xactly as above until the algorithm ends.
,I Fig. \. Initialization: the frame (3:' , y') turns with the nominal point. Mo on the .T' -axis until the convex
hull touches the y' -axis here on the negative side.
Fig. 4. ith iteration: The vertex Mi+l touches the axis on the positive sid~ .
Fig. 2. As above, but the convex h'ull touches th e yl -aris
on the.
positiv~
side.
Other wise, let us denote D'l the vertex of D the image of which touche~ the axis at w = W2. To compute this frequency let us define Xr(w) R(P(jw, Dil» and YF(W) = <;}( P(jw, D i , From Equation (7) tile frequency W2 is the first real root (> wt) for whi ch we ha.ve X'(w) 0:
».
=
=
X'(w) = Yp(w) !lI( P(jw, DJ) - Xr(w),>(P(jw, D;))
(10)
considering all vertices D j excep t D ,1 . The two possibilities that might be encountered are illustrated in Figuret; 3 and 4. First, the image M2 (= P(jW2, D j;))) of Di 2 is on th e negative side of the yl-axis,
Comment. ,"1. 1. Computation of the stability maryin. A lower bound of the stability ra.dius can be computed by di chotomy using at each step the robust stability test of the proposed procedure. Comment. 3.2. Conse1"Vativm. Very often , t he dichotomy described in Comment 3.1 gives the exact value of the stab ility margin. The exact ma rgin is obtained when the edge which tonches t he origin is such that the two vertices of the parameter box defining this edge are adjacent (only onc coordinate differs). Comment. 3.S. An additional intermediate step. To keep the presentation of the algorit.hm as short as possible wc have omitted a step that we use between Step 2 and 3. In the same spirit of Steps 2 a.nd 3 we propose to coo-
3407
sider after Step 2 a frame which turns with its ';y" -axis being parallel to the segment [M; Md. Up to 50 % computing time can be sa.ved in that way.
Comment. ,1.4, A ~ufficie.nt cnndition for' Ilunstability" Often, the frequen cy at. which the convex hull touches the origin belongs to a set of a few freq uencies known a priori. A fast test for uIlstability (dichotomy: see Commp-nt 3.1) COILsits of checking if at Ofi(: of these frequencies, each quadrant of t he complex plane contains at least the image of one vertex.
4. NUMERlCAL JUSTIFICATION The organization of the proposed algorithm being heuristic, many numerical examples Inllst be considered to evaluate its efficiency. All t he rC'$ ults given here are obtained by using interpreted language (Matlab) .
A simple num erical example. The first numerical example considered for making co mparisons was generated as follows: t en eigenvalues l/ Jij + jVl i are chosen at random with VX i and Vlj being random number between zero and one. T he matrix A E R 20x:W is the corresponding Schur diagonal matrix. It comes out that t he matrix A has ten complex modes with damping ratio from 0.13 to 0.99. B and C (see Equation (2) are matrice, ofrandotD numbers between zero and o nc. The number of unknown paramet ers ut is considered form one to fourteen i. e. we consider the first m ::: 1) .. . , 14 columns (resp. rows) of B (resp . C). The com puting tim e necessary for checking stability is given in Table 1. m
rvtethou 1
I
0. 185
2 3 4 5 G
7 8 9 10 Ij 12 U
14
ra.tio
Method 3
0.41s
2.3
0 .265 0.38s
1.285 4.4 18
;U
O.67.~
3.4 3 .7
1.285 2.47."1 [•.05.'1 9.75.'1
1 6.,') .~
Imn 03s -1mn 09.'J 16 mn 255 lit OGmn
3.9
3.0 4.0 4 ."
19.~
.18....
Imn 17.)'
2mn :J3s 5mn 10" lOmn 168 20mrr.
As expected this ratio become.,; equal to 4 for method 1 and equal to 2 for Method 3. For this class of c:onsidereu systems, only o ne change of frame is required to test stability. In other words Lemma 2 suffices. This property will no t be true for systems with low damped modes , at lea..<;t when stability is checked in the vicin ity of the stability margin. Note that using standard in tt~ rpreted software, no numerical difficulties were enconntered at least for 20th degree polynomials and 14 unknown parameters. Examples with low damped modes. A bank of random systems generated as follows is considered. Five pairs of complex eigenvalues a re chosen as Vlli + jVl i where V/li h:i a. random number between - 1 and 0, V{i is a random number between 0 and 20 . An additional large negative real eig-envaiue is alf:;O considered. From these eigen values A E 1R 11 xlI is defi ned as the corresponding Schur diago nal form. The matrices Band C of (2) are generated at ran dom. For eadl system , the stabili ty margin was computed . The analysis of compu tin g tinle given in Table 2. It corresponds to parameter boxes diameter equal to 90% of the stability margin. m
Method 1
Mth 2 wct Meth o,i 1
Method 3
3
0.80 s
2
l!
5
4. 08 11.6
D.5{ ) O.6t-
O.55s
4
.~
O. 3 :~
5 5
6
55 .6
.~
7 8 9 10
3mn 50
,~
0.4~'
].5 .... 3 .3 s 7.88
O.l t)
54.'1
38 .~ 18 s '0,
11mn51s
D .2~.
47,.,m 10 8
0.0<1
3h 50rftn
0.5 ';
7 30 9 I
10
ra.tio
Table 2. Computing time 1.5 1.8 1.. 1.9 2.0
1.9 2.0 2.0 2.0 2.0 2.0
:l.O 2.0
Table 1. Com puting time For example, if 5 minutes can be accepted for testing stability, Method 1 is efficient up to 7 parameters while Method J can be used up to 12 parameters. Colu mns two and four give the ratio of the req1lired computing time each time one ad ditional pa rameter is considered,
The third colum n of this arra y concerns Method 2. It is the ratio of the uumber of edges th at must be tested (as if Method 1 were used) out of the tot al number of edges. If logical test may really be neglected, the computing time nee
=
In a ny case, Method 3 is mu ch more efficient. The efficiency of Ivlet hod 3 is even better whfm the same analysis is made using a pararnder box diameter equal to
3408
half the stabilit y margin (instea
5. CONCLUSION We have propose d in tbis paper a method for verifyin g the stabilit y of a polytop e of polynom ials. Our frequen cy a.pproach, which has given us a. further insight inside the prohlem , is based on two ideas: - frequen cy depend ent cllanges of coordin ates are applied in the comple x plane. - one verifies that the continu um of convex hulls does not contain the origin of the comple x plane for w E [O,oo( by verifyin g this propert y for w E [0, w,j U [Wl' W2] u .... The examples of §4 ha.ve iIIustrateci the reducti on of the comput ational burden when compar ing our method to other existing method s. rvlorcov er I the results reveal general ly ouly slightly conserv ative. The initial aim of the paper was to evaluat e rohllstne5~ of a closed loop without. frequenc.y griddin g, beca use this techniq ue may reveal unreliab le when conside ring lightly damped bending modes. III the CO[ltt~xt of J.t-analysis, a first method is tu directly comput e t he maxima l structured singula r value (s.s.v .) OVE~r th e frequen cy range: this maxima l s.s. v. can be obtaine d a.'i the solution of an augmen ted skewed p.-probl em corresp onding to a statespac.e approa ch (see (Ferrere s and M 'Saad , 1994) and include d referenc es). \Ve have conside red in thLe; paper an alterna tive mE:!thod ~ ba!:ierl on the Edge and Mappin g Theore ms.
6. REFER ENCES A.C. Bartlet t , C.V. Hallot, and Hua ng Lin (1988) . Root location s of an entire polytop e I)f polyno mials: it suffices to check the edges. In Mathem atic" of Control, Signals "nd System s, 1:61-71. S. Bialas (1985). A necessa ry and sufficien t conditi on for the stabilit y of convex combin ations of stable polynom ials of matrices. B ull. Polish. Acad. Sci., 33(9-10):474480. H. Bougue rra, B .C. Chang, H.H. Yeh, and S.S. Banda (1990). Fast sta.bility checkin g for the convex combin ation of stable polyno mials. IEEE Tmnsactions on Automatic Control, AC-35( 5) :58fi-588. H. Chapcl lat ant! S.P. I3hattacharyya (1989). An alterna tive proof of kha.ritonovls theorem . IEEE Transac tions on Automatic Control, AC-34( 4):448- 450. R.L. Dailey (1990). A new algorith m for the real structured singula r value, in ProCJ.:ed ings American Control
Conference.
s. Dasgup ta (1988). Kharito m,v's t.h eorem revisited. System and Control Letter, 11:381 - 384. G . Ferrere s and M. M'Saad (J 994) . Direct comput ation
of the maxima.l s.s,v. over tile fretluen cy ra nge using the 1) tool. In Pmc. of the 3Srd Con! on Decision and
ContT.. J . Freund enberg and B. Morton (1992). Robust control
of a booster vehicle using h ')O and ssv t.echniques. In Proc. Conference on Decision and Con trol, pages 24482453. A. Packard and J.C . Doyle (1993). The comple x structured singula r value. Automa tica, 2~(1).
In the context of .u-ana1ysis, the problem is that the J-Lupper bound corresp onding to the state-sp ace approac h is general ly nwrc conserv ative than the J-L·upper bound obtaine d using a freque ncy griddin g (Pacbr d and Doyle, 1993). This problem does not, exist iD the context of our method , since the manipu lation of the continu um of convex hulls is equival ent to an infinite dimensio nal frequen cy griddilJg. On the other hand, the (:omput ation of the It-uppe r bound is polynom ial time, wherea s o ur method is exponenti a l time. As a matte r of fa.d, both method s appear rather comple mentar y. H the method involvin g the .u-uppe r bound is applica ble for larger amount s of paramete rs , our method appea.rs ne verthele ss well adapted to problemf:i with reasona ble amount s of parame ters.
3409
Copyright © 1996 IFAC 13th Triennial World Cong ress. San Francisco. USA
REM ROBU STNE SS ANAL YSIS USING THE MAPP ING THEO DING GRID Y WITH OUT FREQ UENC Gilles Ferrer es " J ean-Fran~ois Magni ., 1 Edouarri • Automatic Control Departm ent (DERA ), CERT- ONERA , 2 Avenue France Cedex, e B elin, BP 4025, F31055 Toulous
comput in g a. lower Abstra ct. This pa per proposes a new numerical techniq ue for The te(~hnique is g. griddin cy frequen bound of th e robust stabilit y margin , withou t t hp. numeric al that way a such in ed organiz is based on the Ma pping Theore m and unlikc in edges of number the of instead vertices of comple xity is related to the number a ban k ring conside trated demons is cy efficien al numeric The o ther known techniq ues. le for tradab he to of unce rtai n system s. T he stabilit y margin comput a tion is shown damped }' lightl several h t more than twelve uncerta in parame ters even for system s wi modes. tric varia.ti ons K e yword s. Sta bility margin , stabilit y test, robustDess, parame 1. INTRO DUCTI ON
Checkin g the sta bility of a family of polynom ials pes, 6), wh ere the vector 8 of Illlcerta in parame ters belongs t o an hype rrectan gle D, is a classical problem in ro bustnes s analysis. Exact solut ions a re availa ble when the coefficients of P(s, 6) belon g to interval s or arc affine functions of 6 (n amely the Kharit onov and Edge theorem s sec (Dasgu pta , 1988), (Bartle t! et al. 1987)). The problem is howeve r mu ch mo re difficult when the coefficie nts of P(s, c5 ) a re multilin e;u- or polynom ial functio ns of b. In the context of a multilin ear depend ency, an ou t.erboundi ng set of the family of polyno mials P{ s, h) ~ na mely the polytop e conv (P(s , D ;)) where D ; is a vertex of the hyperre ctanglc D, can be conside red. It ~uffices t hen to check the stabili ty of the (expose d) edges so aB to check t he stabili ty of t he poly tape of polynomia.ls. There is howeve r a combin atoria l explosi on with the size /11 1 of the problem , because of the 2 - (2'" -1 ) edges of the t.un para meters). An uncer polytop e (m. is the number of in (Bouguerra et al, cl propo~e heuristi c method has been burden: conside rtional a comput the 1990) for redu cing is perfor med, test logical a e, polytop he t of ing an edge for the staion t condi nt sufficie a to onds corresp ch whi hility of th e edge, further comput ation is thus necessa ry only if t his logical test is not verified . It can be howeve r noted that t he compn ta tio nal comple xity is still more 7 11 or less related to 4HI ::::: 2/11 "'\ (2 - 1), since it is still necessa ry to con ~i cler a ll edges of {,he polytop e. Morees I Suppor ted by Di rection d es Rcd lert:hes et Etudes Te<::hniqu
over, the comput ational require ment for verifyin g this logical t est may reveal non negligib le, depend in g on the im plemen tation. An a priori attracti ve solution is to combin e a frequen cy griddi ng wit h the zero exclusio n principle. Let conside r the con vex hull of t he points P(jw, D;) in the com plex j w (remem ber D j is one the 2Tf1 vertices of plane at s the hyperre clangle D ). The poly to pe of poly nomia ls is sta hle if a nd only if one elemen t of the polytop c is stable a.nd if the convex hull does not contain the origin of th e comple x plane a t all frequencies w. Noting N the number of frequen cy points, the comput ational comple xi ty is now N 2"1 . Neverth eles!-S, even if a frequen cy griddin g is rather usual ill robustn ess a nalysis pro blems (especia lly in the context of I.£~a.naly sis) . t his techniq ue may reveal unreliable in the specific context of lightly da mp ed b e nd~ ing mo nes, since it appe.ars possibl e t o miss the critica l frequen cy in t his case (FreuIlfle.nberg and Morto n, 1992).
=
The key idea of this paper is t o directly manipu late the continu um of convex hulls corresp onding to w E [O, oo{ and t o check whethe r the origin of the co mplex plane belongs t o this continu um. As a ma tter of fact , this freque ncy a pproacl, gi ves U H further insig ht inside our problem , and it is especially possible to reinterp ret the suffi cient conditio n in (Bougu crra et ai, 1990)) in this new co ntext . A sufficien t condit ion for t he sta bility of th e polyt ope of polynom iah, is then proposc d: its (".omputatio nal comple xity is 2111 . To check the stabilit y of the Volytop e of polyno mials (or equival ently to check wh et her t he continu um of c.onvex hulls contain s the origin of t he comple x pla. ne), it is necessary in the worst
(DRET)
3404
case t.o apply 2H1 times this sufficient condition. Even if the worst case complexity is still 4"\ numerous examples treated in §4 tend to prove tha.t the mean complexity is closer to 2fl1. As an important consequence, it seems pot>Siblc to deal with problems involving up to 12 uncertain parameters with a recwonable computational burden (about 5 minutes for polynomia.ls with degree 20) instead of 7 parameters with a. hrute application of the Edge Theorem. The paper is organized as follows. In §2 are sketched existing methods which are not based 011 frequency griddingo §3 presents the proposed method for computing the lower bound of the robustness margin, while the examples of §4 emphasize the practical interest of the method.
In view of this elementary result, two existing t echniques are briefly ske tched. M e thod 1. In order to check t he stability of a polytope of polynomials, (Bialas, 1985) f:onsiders a test which enables to check the ,tability of each edge . Bialas test requires the computation of thi~ eigenvalues of a n x n matrix deriveu from Hurwitz matrices. Here is used an alternative t est (see Bouguerra. et al, see also (ChapeUat a.nd Battacharya , 1989)) which is an immediate COll f)P-qllence of LeUlma 1. Let us consider two vertices DJ: aJld Dj of the pa.rameter box D. We have to check whether the segm ents [0, P(jw , D d ] and [0, P(jw, D d ] are collinear at some (n'! al) frequency w, 01' similarly have the same slope: 'S(P(jw, D,)) 'S(P(jw, Dd) fR(P(jw, Dd) - 1R(P(jw, Dc))
2. PROBLEM STATEMENT, REVIEW OF EXISTl,' W METHODS
Thiti is equivalent to lo ok for positive real roots of th e w-polynomial
Let us consider a polynomial
!R(P(jw, DtlP(P(jw, Dd) -l?( P(jw, D, ))"(P (jw, Dtl) (3)
bED where D is an hyper rectangle of Rill i. c. a set of points defined by D = {6= (61 ,62 , ... ,6",) , 6,- ::;6i ::;6t}
6t .
for given bounds 6; and The coefficients ado), i = 1, ... , n are assumed to be multiltnwT' junctions of 01, ... , li,rl' The nominal poin t in the parameter box will be denoted Do, it corresponds to the vector 0.5((01 + lit) •... , (D~ + 6;t)). The 2'" vertices will be denoted Dj. The polynomial of (1) is generally given under the following {orm:
P(s,b) =det (s I-A - B diag(6" ... , hm )C) (2 ) where A E lR":<" , B E R" Xrt' and C E
JRJ/lXTI.
The value set at EL given point .~ is the domain of the complex plane defined as pes, D). The following result, which is a direct consequence of the Boundary Crossing and Mapping Theorems, will be used. Lemma 1. Let us a.c;sume that the nominal polynomial p es , Do) is stable. The system is stable if at the frequency w = 0 the value se t. P (O, D) does not contain the origin and if one never encounters, for all real values of w > 0 , thp. case where two verti ce~ of the convex hull and the origill are on the same straight line, with the origin in b etween the two vertice.').
(Note that if we consider this polynomial with respect to w2 its degree is equal to n -I , so the involved computation is of the same order than with the original test of Bialas.) If t he polynomial of Equation (3) admits a positive real root Wo . it remains to check whether the origin of the complex plane is b etwl~en the points P(jwo, D{) anu P(jwo. Dd. The answer to this qu estion is affirmative if and only if: ~(P(jwo , D,))fR(P(jwo,Dk)) < 0 or fR(P(jwo, D,))'R(P(jwo , D .)) 0 and "'(P (jwo , D,))'S(P(jwo, Dk)) ::; 0
=
(4)
It comes out that the complf·xity of Bia las test and of the altern ative test is similar: computation of the roots of 2(11 -1)(2" -1) polynomials of degree n is required in both cases. Method 2. A similar idea. is used in (Bouguerra et at, 1990), still considering 211- 1(21< -1) edges but the computat ion of the roots iD replaced in most cases by a logical test which is assumed to be faster. I nstead of considering the explanations given in (Bouguerra et al, 1990) let us giVf~ o ur point of view concerning this technique: in that way the proposed technique of §3 will appear (1.<; a 'tep furth er beyond Methods 1 and 2. Clearly edgeR of the convex hull, say the segm.mt [P(jw , Dd , P(jw , Dd] , v:m not touch the origin for w > 0 (so that Lemma 1 enables to conc.1 ucle) if these points are in a same half plane for all values of w. Thi~ propert.y is satisfied if for example
3405
!R(P(jw, Dd) and !R(P(jw, D.» ,ame sign \/w or
(S(P(jw, Dd) alld (S(P(jw, D.» ,ame sign \/w
• its 4oI:t _axis rl is collinear CO the oriented segment OMI_', the "y_axis" is such that the frame is direct. It is clear that the coordinates of the point M which were (X, Y) her.ome (X', V')
Therefore, Bouguerra td al propose to proceed as follows. If for some interval oC frequencies we have :
" _ XpX+YpV . Y'_ -YpX+XrV X X',.. +" y2 ' X 'p+p V'
!R(P( jw, D,)!R(P(jw , Dd) < 0 aDd (S(P(jw , DJ)\J(P(jw, Dd) < 0
As only the signs of X ' and Y' will be needed} we shall consider the following simplified coordin a tes:
(5)
X'=XrX+YtY; Y'=-Y,X+Xr Y (6) check conditions (3) a nd (4), otherwise the edge is known to be stable. The numerical complexity corresponds to t he computation of the positive real root~ of 2"+1 polynomials with degree d ose to n / 2 (th~e polynomials are !R(P(jw, D;» and U(P(jw , D;» , i = 1, . .. , 2"). From the knowledge of these roots , tlw logica l test of Equation (5) is considered for the 2"-'(2" - 1) edges. For edges such that (5) is satisfied in some intervals of frequencies, StnI'm Theorem enables to check the existence of real roots inside each suc.h interval. Sturm Theorem will be replaced here by the test based 011 Equations (3) and (4) . III the procedure proposed in §3 (Method 3), we shall try to avoid the logical t est concerning all edges because this test is not as fast as it might appear at first sight (depending on the implementation). We shall also improve the test of Equation (5) by considering a frame which turns with th e value set so that chcr:king the fact t ha.t P(jw , D,) and P(jw, D. ) belong to a same half plane becomes much more efficient. In fact instead of eliminating edges} wc eliminate intervals of frequencies for which we know that. the nmtinnum of convex hulls does not contain the origin. If we can eliminate all frequencies, it means that the system is robustly stable. Otherwise the first frequency at which th e (:ontinuum tou ches the origin is characterized.
• its "y-axis" is collinear to the oriented segment OM F , the "x-axis" is such that the frame is direct. It is dear that the coordina.tes of the p oint M become here (X' , V' ) (with simplification as ahovc):
X'=Y"X-X,.,Y; Y'=X,..X+YpY (7) First step. It is necessary to check that, at frequency
w _ 0, the value set pen, D) does not contain th e origin . Tt is well known th at P(O,D) = [ Min ;P (O, D, ) Maxi P(O , D;)I. SO, if the o rigin belongs to this set., P( .. , 6) is not robustly stable} otherwise go to next ste ps. Second ste p. In order to use Lemma 1 in an efficient way} it seems to be worth c( !O~idcring a frame which turns with the image of the lw minal parameter vector Do. This idea can be justified by considering the locus corresponding to the nominal ('ase. It is a point which remains (of course) on th e "xaaxis'l • When th e box around Do is no t too large, the valu e set will remain in the half plane X'( w) > 0 (see (6». So for s uch small boxes a sufficient c.ondition for stability is X'(w) > 0 for all positive fr~ucncies. Ivlore precisdy we have the following !)ufficicnt condition for sta bility:
= =
Lemma 2. Let XFCw) !R(P(jw , Do» and Vdw) If for i 1 ... 2" , the polynomials X,..(w) R(P(jw, D;» + Y" '(w) ~~(P(jw, D;» have 110 positive real root, the polynomial P( .• .I) is robustly stable.
~(P(jw, Do».
3. A NEW ROBUST STABILITY TEST Our main objective in this section is to derive a robust Rtability test based on value set idea.~ but without frequency gridding, also we intent to avoid numerical complexity related to the number of edges 2 11 - l (2" -1) (m denotes the number of uncertain parameters) by keeping it as close as possible to the number of vertices 21n. Change of coordinate in the complex plane. Before giving a description of the proposed stability test, we shall consider a change of variable. Let us consider a point M in the complex plane with affix X + JY. An other point M r with affix XF + .fYv i~ used to define a frame in the c:omplex plane £U:i follows:
So to check stability} if Lemma 2 is not !:iufficicnt to conclude, it means that the w 2 -polynomial X'(w) ~ xp(w)31(P(jw , D . JJ
+ YF(W)"'(P(jw, V i ))
(8)
admits some positive real roo ts, let WI be such that wf is the smallest of these roots. \Vhen w = w. th e image of one vertex touches the ll.y" -axis. Two pOSl;)ibilities are depicted in Figures 1 and 2. In the first case} we must find which vertex of D CO rTl ~8ponrls to the point (denoted. M 1 in Figure 1) with maximal slope in th e moving frame. More precisely WI ! look for the vertex Dil (Mc =P(jwJ,D;,)) such thaj SeD;,) defined by :
3406
-Yr(w,)'R(P(jw" Di, )) + XF(W, )'>(P(jW" Di, )) X r(w,)lI1(P(jw" Di,)) + Yr(w,)'>(P(jwJ, Di, ))
!'i
(9)
I
is maximal. At the next step , this point will be used to define the moving frame. In the second case (see Figure 2) the point 111, is ob tained directly. Additional steps. In the additional steps, the moving frame turns now with the image of Dil defined at the previous step. But DOW, it is its ""y'.a.xis" as in (7) which conta ins the image of Dil i.e. AIl P(jw , D i \) . As previously we are interested in the first frequency (,.1)2 at which the value set touches the new "y' -(tXis" . If such a frequency does not exist, the polynomial p es, a) is robustly stable. The algorit.hm ends.
=
Fig. 3. ith iicrntion: The ve1tex M i + 1 touches the axis on the negative side. therefon~ (t he boundary of) the convex hull contains the origin. So the algorithm ends . Generally we cannot conclude about instability except for some pairs of" vertices defining the segment [Ml . W2"1, see comment 3.2.
Second, the image of D j 2 is OIL the other side. In which case we define a new frame which will turn with M2 011 its Uyl -axis" i:LDd proceed t·xactly as above until the algorithm ends.
,I Fig. \. Initialization: the frame (3:' , y') turns with the nominal point. Mo on the .T' -axis until the convex
hull touches the y' -axis here on the negative side.
Fig. 4. ith iteration: The vertex Mi+l touches the axis on the positive sid~ .
Fig. 2. As above, but the convex h'ull touches th e yl -aris
on the.
positiv~
side.
Other wise, let us denote D'l the vertex of D the image of which touche~ the axis at w = W2. To compute this frequency let us define Xr(w) R(P(jw, Dil» and YF(W) = <;}( P(jw, D i , From Equation (7) tile frequency W2 is the first real root (> wt) for whi ch we ha.ve X'(w) 0:
».
=
=
X'(w) = Yp(w) !lI( P(jw, DJ) - Xr(w),>(P(jw, D;))
(10)
considering all vertices D j excep t D ,1 . The two possibilities that might be encountered are illustrated in Figuret; 3 and 4. First, the image M2 (= P(jW2, D j;))) of Di 2 is on th e negative side of the yl-axis,
Comment. ,"1. 1. Computation of the stability maryin. A lower bound of the stability ra.dius can be computed by di chotomy using at each step the robust stability test of the proposed procedure. Comment. 3.2. Conse1"Vativm. Very often , t he dichotomy described in Comment 3.1 gives the exact value of the stab ility margin. The exact ma rgin is obtained when the edge which tonches t he origin is such that the two vertices of the parameter box defining this edge are adjacent (only onc coordinate differs). Comment. 3.S. An additional intermediate step. To keep the presentation of the algorit.hm as short as possible wc have omitted a step that we use between Step 2 and 3. In the same spirit of Steps 2 a.nd 3 we propose to coo-
3407
sider after Step 2 a frame which turns with its ';y" -axis being parallel to the segment [M; Md. Up to 50 % computing time can be sa.ved in that way.
Comment. ,1.4, A ~ufficie.nt cnndition for' Ilunstability" Often, the frequen cy at. which the convex hull touches the origin belongs to a set of a few freq uencies known a priori. A fast test for uIlstability (dichotomy: see Commp-nt 3.1) COILsits of checking if at Ofi(: of these frequencies, each quadrant of t he complex plane contains at least the image of one vertex.
4. NUMERlCAL JUSTIFICATION The organization of the proposed algorithm being heuristic, many numerical examples Inllst be considered to evaluate its efficiency. All t he rC'$ ults given here are obtained by using interpreted language (Matlab) .
A simple num erical example. The first numerical example considered for making co mparisons was generated as follows: t en eigenvalues l/ Jij + jVl i are chosen at random with VX i and Vlj being random number between zero and one. T he matrix A E R 20x:W is the corresponding Schur diagonal matrix. It comes out that t he matrix A has ten complex modes with damping ratio from 0.13 to 0.99. B and C (see Equation (2) are matrice, ofrandotD numbers between zero and o nc. The number of unknown paramet ers ut is considered form one to fourteen i. e. we consider the first m ::: 1) .. . , 14 columns (resp. rows) of B (resp . C). The com puting tim e necessary for checking stability is given in Table 1. m
rvtethou 1
I
0. 185
2 3 4 5 G
7 8 9 10 Ij 12 U
14
ra.tio
Method 3
0.41s
2.3
0 .265 0.38s
1.285 4.4 18
;U
O.67.~
3.4 3 .7
1.285 2.47."1 [•.05.'1 9.75.'1
1 6.,') .~
Imn 03s -1mn 09.'J 16 mn 255 lit OGmn
3.9
3.0 4.0 4 ."
19.~
.18....
Imn 17.)'
2mn :J3s 5mn 10" lOmn 168 20mrr.
As expected this ratio become.,; equal to 4 for method 1 and equal to 2 for Method 3. For this class of c:onsidereu systems, only o ne change of frame is required to test stability. In other words Lemma 2 suffices. This property will no t be true for systems with low damped modes , at lea..<;t when stability is checked in the vicin ity of the stability margin. Note that using standard in tt~ rpreted software, no numerical difficulties were enconntered at least for 20th degree polynomials and 14 unknown parameters. Examples with low damped modes. A bank of random systems generated as follows is considered. Five pairs of complex eigenvalues a re chosen as Vlli + jVl i where V/li h:i a. random number between - 1 and 0, V{i is a random number between 0 and 20 . An additional large negative real eig-envaiue is alf:;O considered. From these eigen values A E 1R 11 xlI is defi ned as the corresponding Schur diago nal form. The matrices Band C of (2) are generated at ran dom. For eadl system , the stabili ty margin was computed . The analysis of compu tin g tinle given in Table 2. It corresponds to parameter boxes diameter equal to 90% of the stability margin. m
Method 1
Mth 2 wct Meth o,i 1
Method 3
3
0.80 s
2
l!
5
4. 08 11.6
D.5{ ) O.6t-
O.55s
4
.~
O. 3 :~
5 5
6
55 .6
.~
7 8 9 10
3mn 50
,~
0.4~'
].5 .... 3 .3 s 7.88
O.l t)
54.'1
38 .~ 18 s '0,
11mn51s
D .2~.
47,.,m 10 8
0.0<1
3h 50rftn
0.5 ';
7 30 9 I
10
ra.tio
Table 2. Computing time 1.5 1.8 1.. 1.9 2.0
1.9 2.0 2.0 2.0 2.0 2.0
:l.O 2.0
Table 1. Com puting time For example, if 5 minutes can be accepted for testing stability, Method 1 is efficient up to 7 parameters while Method J can be used up to 12 parameters. Colu mns two and four give the ratio of the req1lired computing time each time one ad ditional pa rameter is considered,
The third colum n of this arra y concerns Method 2. It is the ratio of the uumber of edges th at must be tested (as if Method 1 were used) out of the tot al number of edges. If logical test may really be neglected, the computing time nee
=
In a ny case, Method 3 is mu ch more efficient. The efficiency of Ivlet hod 3 is even better whfm the same analysis is made using a pararnder box diameter equal to
3408
half the stabilit y margin (instea
5. CONCLUSION We have propose d in tbis paper a method for verifyin g the stabilit y of a polytop e of polynom ials. Our frequen cy a.pproach, which has given us a. further insight inside the prohlem , is based on two ideas: - frequen cy depend ent cllanges of coordin ates are applied in the comple x plane. - one verifies that the continu um of convex hulls does not contain the origin of the comple x plane for w E [O,oo( by verifyin g this propert y for w E [0, w,j U [Wl' W2] u .... The examples of §4 ha.ve iIIustrateci the reducti on of the comput ational burden when compar ing our method to other existing method s. rvlorcov er I the results reveal general ly ouly slightly conserv ative. The initial aim of the paper was to evaluat e rohllstne5~ of a closed loop without. frequenc.y griddin g, beca use this techniq ue may reveal unreliab le when conside ring lightly damped bending modes. III the CO[ltt~xt of J.t-analysis, a first method is tu directly comput e t he maxima l structured singula r value (s.s.v .) OVE~r th e frequen cy range: this maxima l s.s. v. can be obtaine d a.'i the solution of an augmen ted skewed p.-probl em corresp onding to a statespac.e approa ch (see (Ferrere s and M 'Saad , 1994) and include d referenc es). \Ve have conside red in thLe; paper an alterna tive mE:!thod ~ ba!:ierl on the Edge and Mappin g Theore ms.
6. REFER ENCES A.C. Bartlet t , C.V. Hallot, and Hua ng Lin (1988) . Root location s of an entire polytop e I)f polyno mials: it suffices to check the edges. In Mathem atic" of Control, Signals "nd System s, 1:61-71. S. Bialas (1985). A necessa ry and sufficien t conditi on for the stabilit y of convex combin ations of stable polynom ials of matrices. B ull. Polish. Acad. Sci., 33(9-10):474480. H. Bougue rra, B .C. Chang, H.H. Yeh, and S.S. Banda (1990). Fast sta.bility checkin g for the convex combin ation of stable polyno mials. IEEE Tmnsactions on Automatic Control, AC-35( 5) :58fi-588. H. Chapcl lat ant! S.P. I3hattacharyya (1989). An alterna tive proof of kha.ritonovls theorem . IEEE Transac tions on Automatic Control, AC-34( 4):448- 450. R.L. Dailey (1990). A new algorith m for the real structured singula r value, in ProCJ.:ed ings American Control
Conference.
s. Dasgup ta (1988). Kharito m,v's t.h eorem revisited. System and Control Letter, 11:381 - 384. G . Ferrere s and M. M'Saad (J 994) . Direct comput ation
of the maxima.l s.s,v. over tile fretluen cy ra nge using the 1) tool. In Pmc. of the 3Srd Con! on Decision and
ContT.. J . Freund enberg and B. Morton (1992). Robust control
of a booster vehicle using h ')O and ssv t.echniques. In Proc. Conference on Decision and Con trol, pages 24482453. A. Packard and J.C . Doyle (1993). The comple x structured singula r value. Automa tica, 2~(1).
In the context of .u-ana1ysis, the problem is that the J-Lupper bound corresp onding to the state-sp ace approac h is general ly nwrc conserv ative than the J-L·upper bound obtaine d using a freque ncy griddin g (Pacbr d and Doyle, 1993). This problem does not, exist iD the context of our method , since the manipu lation of the continu um of convex hulls is equival ent to an infinite dimensio nal frequen cy griddilJg. On the other hand, the (:omput ation of the It-uppe r bound is polynom ial time, wherea s o ur method is exponenti a l time. As a matte r of fa.d, both method s appear rather comple mentar y. H the method involvin g the .u-uppe r bound is applica ble for larger amount s of paramete rs , our method appea.rs ne verthele ss well adapted to problemf:i with reasona ble amount s of parame ters.
3409