- Email: [email protected]
The robust root locus
0005-1098/90 $3.00 + 0.00 Pergamon Press pie ~) 1990 International Federation of Automatic Control
Autoraatica, Vol, 26, No. 2, pp, 283-292, 1990 Printed in Great Britain.
The Robust Root Locus* B. R. BARMISH~'§ and R. TEMPO~tll
A simple technique is described for generating the root loci of a feedback system which includes perturbations entering affine linearly into the coefficients of the plant.
KeyWordsmRoot locus; perturbations; robustness; value set. Abstract--In this paper, a simple technique is described for generating the root loci of a feedback system which includes perturbations q ~ Rt entering affine linearly into the coefficients of the plant. Denoting the perturbed plant by ~(s, q) and the compensator by C(s), we address the following problem: given a bounding set Q = R~ for q, find the locus of points z in the complex plane such that 1 + k~(z, q)C(z)= 0 for some k >-0 and some q a Q. Such points z are said to lie on the Robust Root Locus. One of the strengths of the technique presented in this paper is that it avoids the "combinatoric explosion" synonymous with gddding the/-dimensional set Q and plotting a large number of ordinary root loci associated with gfidding the points. Instead, the technique given here exploits only a 2-dimensional gridding of a bounded subset of the complex plane. For the more general case when perturbations enter multilinearly into the coefficients, a certain convex hull robust root locus is constructed. Finally, the robust root locus technique is illustrated using a model of a supersonic transport plane.
of the plant m a y be perturbed. Hence, the " t r u e " root locus may differ from the root locus generated using the so-called nominal system. In this regard, the main focal point of this p a p e r is the generation of root loci with respect to variations in m o r e than one parameter. This work can be viewed as a systematic approach to robustness problems which were alluded to in the fifties and sixties; e.g. see Truxal (1955). In the later work of Z e h e b and Walach (1977), a two-parameter root locus p r o b l e m is considered and rather specific assumptions (motivated by circuit theory) are m a d e about the class of perturbations. A n u m b e r of papers following Z e h e b and Walach's work deal with the so-called zero set concept; e.g. see Z e h e b and Walach (1981) and Fruchter, et al. (1987). It is seen that t~he zero set provides a rather general f r a m e w o r k for dealing with multi-parameter root loci. That is, for each value of the gain k, a function can be given whose zero set corresponds with the cross section of the root locus at k. Hence, in principle, a multi-parameter root locus can be generated by computing these cross sections as k varies. In practice, however, there are difficulties associated with numerical computation of the zero set. By strengthening the hypotheses on the class of allowable perturbations, we obtain an enormous simplification in the characterization of multi-parameter root loci. To this end, attention in this p a p e r is restricted to the case when the perturbations enter the plant coefficients either linearly or multilinearly. This leads to the following type of result: a point z ~ C can be " t e s t e d " for m e m b e r s h i p in the root locus by simply checking whether or not the point s = 0 in the complex plane lies inside or outside a specially constructed 2-dimensional convex polygon. Hence, one avoids the " c o m binatoric explosion" in computation associated with an /-dimensional gridding of the rectangle
1. INTRODUCTION AND FORMULATION T o MOTIVATE this paper, recall the classical R o o t Locus Problem: given a nominal plant transfer function ~ ( s ) and a fixed c o m p e n s a t o r C(s), find the locus of points z e C such that 1 + k~(z)C(z) = 0
(1)
for some k-> 0. That is, find the locus of the closed loop pole locations as a function of the gain k; e.g. see Krall (1970) for a survey of classical root locus techniques. The take-off point for this p a p e r is the following fact: when performing robustness analysis and design, various physical p a r a m e t e r s * Received 5 July 1988; revised 24 April 1989; received in final form 14 August 1989. The original version of this paper was not presented at any IFAC meeting. This paper was recommended for publication in revised form by Associate Editor V. Ku~era under the direction of Editor H. Kwakernaak. ~"Department of Electrical and Computer Engineering, University of Wisconsin-Madison, Madison, WI 53706, U.S.A. , CENS-CNR, Politecnico di Torino, Corso Duea degli Abruzzi 24, 10129 Torino, Italy. Author to whom all correspondence should be addressed. § Supported by the National Foundation under Grant No. ECS-8612948 and NATO-CNR Fellowship Program. II Supported by funds of CENS-CNR of Italy. PdlTO ~ : 2 - 6
283
284
B . R . BARMISH and R. TEMPO
Q. That is, instead of a potentially high dimensional gridding of Q, one can test candidate z points in a two-dimensional subset of the complex plane. The type of "zero inclusion" condition above goes back a long way in the robust stability literature (see, e.g., Zadeh and Desoer, 1963) and has been revived in the last few years (e.g. see Saeki, 1986; de Gaston and Safonov, 1988; Dasgupta, 1988; Minnichelli et al., 1989). Furthermore, the zero exclusion condition is easy to check because only two dimensions are involved and the set to be checked is a convex polygon. The first distinguishing feature of this paper is the use of the zero inclusion condition within the context of root locus rather than robust stability as in previous work. The second distinguishing feature is that the technique given here readily lends itself to practical application and computer-aided analysis and design. To make this argument more compelling, we include a rather detailed example involving the tuning of a rate gyroscope using a model of a supersonic transport plane. To study variations in the root loci under multiple parameter variations, we let q e R ~ denote a vector of perturbations with i-th component qi satisfying q/- -< qi -< q,+
(2)
where the q,.+ and q~- are prescribed bounds. Hence, the vector of perturbations q is confined to the/-dimensional rectangle
Q={q:qT<-qi<-q~;i=l
. . . . ,l}.
(3)
To denote the dependence of the plant transfer function on q, we write ~(s, q). We now present the formal definition of the Robust Root Locus (RRL): a point z ~ C is said to lie on the Robust Root Locus (RRL) if
1 + k~(z, q)C(z) = 0
(4)
for some k - 0 and some q e Q. Since it is also important to know the distribution of the closed loop poles for a fixed value of the gain k, we provide the following definition: for fixed k = k*, a point z e C is said to lie on the cross section of the R R L at k* if
1 + k*~(z, q)C(z) = 0
(5)
for some q ~ Q. The main objective of this paper is to provide a computationally tractable method to generate accurately the R R L for the case when the coefficients of the plant are affine linear functions of q. To this end, we do not simply overbound the distribution of the poles in the s-plane--we find the distribution precisely. This
lack of conservatism may seem somewhat surprising because of the complicated relationship between the zeros of a polynomial and its coefficients (e.g. see Marden (1966)). In a sense, this paper can be viewed as a complement to the growing body of literature on parameter space methods (see e.g. Ackermann (1980)). One reason why parameter space methods are powerful is because it is much easier to study the effect of feedback on the coefficients of the characteristic polynomial than on the closed loop poles. One of the earliest papers to fully recognize and exploit this fact is that of Fam and Meditch (1978). This is also recently recognized in Soh et al. (1987) and Biernacki, et al. (1987). Noting, however, that the closed loop poles may be highly sensitive to the coefficients of the characteristic polynomial, a design which "looks good" from a parameter space point of view may have poles which are badly "smeared" over a large region in the complex plane. Hence, having the R R L available, one has the opportunity to check whether or not this phenomenon is occurring; i.e. examine the cross section of the R R L at k = 1. Subsequently, by tuning the loop gain, there is an opportunity to obtain a more desirable distribution of closed loop poles: said another way, the R R L can be used for robustness enhancement. The remainder of this paper is organized as follows. In Section 2, we introduce the necessary notation for the parameterization of the closed loop system with respect to q e Q. Section 3 provides a simple technique to decide whether or not a given point z ~ C lies on the RRL. In Section 4, the R R L technique is reinterpreted in terms of the control system parameters and in Section 5, a crude bound for the cross sections of the R R L is computed. In Section 6, we describe two extensions of the theory. We note that the R R L can be generated with respect to design variables other than a simple loop gain and we extend the R R L technique to the case when the perturbations enter multilinearly into the plant coefficients. In this multilinear case, we show that it possible to generate a so-called Convex Hull Robust Root Locus (CH-RRL). Finally, in Section 7 we present a numerical example and in Section 8, we provide conclusions.
2. BASIC NOTATION 2.1. The extreme points of the rectangle Q As stated in Section 1, we assume that the vector of perturbations q is confined to the /-dimensional rectangle Q. This rectangle has at most L -" 2 t extreme points denoted
The Robust Root Locus
q~, q2 . . . . . qL. The i-th extreme point has j-th component given by qj=q~
or
i__ q~ - q ~÷.
2.2. The plant parameterization For each fixed q e Q, the plant is a strictly proper* single-input single-output system~represented by the transfer function Np(s, q) e(s, q ) = De(s, q-~)
(6)
where Np(s, q) and De(s, q) are polynomials in s of the form
We(s, q) = a,,(q)s"" + a,,_l(q)s"'-' + " . +a~(q)s+ao(q); De(s, q) = s a" + ba~_~(q)s a~-~
(7)
+ . . . + b~(q)s + bo(q). Furthermore, the coefficients a~(q) and b~(q) above are given affine linear functions of the perturbation vector q; i.e. 1
ai(q) = a~o + ~ aiiq~ ./=1
285
We introduce the vector di(q, k ) e R a'+a~ to represent the coefficients of A ( s , q , k ) ; i.e. 6(q,k) has i-th component 6i_l(q,k) and hence, we write d~,+dc-I
A(s, q, k)=Sac÷d'+
~
5i(q, k)s i.
(11)
i=0
In Section 4, we express the vector 6(q, k) more directly in terms of the coefficients ai~, bit, o:~ and fig describing the plant ~ ( s , q ) and the compensator C(s). This type of expression proves useful for numerical computation.
3. KEY IDEAS BEHIND THE CONSTRUCTION OF THE RRL
As emphasized in the Introduction, it is often not feasible to generate an RRL by gridding the /-dimensional set Q. In this section, we see that the computation of the RRL can be achieved using only a 2-dimensional gridding of a bounded subset of the complex plane. This reduction in complexity leads to computational tractability which will be demonstrated using the numerical example presented in Section 7.
for i = 0, 1. . . . . ne and I
bi(q) = bio + Z biiql i=~ for i =0, 1. . . . .
de-1.
2.3. The compensator As usual, the compensator to be applied to the plant ~(s, q) is a proper rational transfer function
C(s) = k Nc(s) Dc(s)
(8)
3.1. Key ideas for the special case of an interval plant with a static gain To motivate the general construction of the RRL, we consider the special case of interval polynomials for Np(s, q) and De(s, q) and a static gain compensator C(s)
k.
(12)
Hence, suppose that each plant coefficient ai(q) or bi(q) depends on one and only one perturbation parameter; say
where k -> 0 is an adjustable gain and
ai(q) = aio + qi+~
Nc(s) = o~cs"C + o~c_~s rI~-~ + • • • + oqs + O~o
for i = 0, 1. . . . .
Dc(s) = sa" + fldc-I S d c - I + " ' " + fll S + ~0
np and
bi(q) = bio + qne+i+2
(9) are prescribed polynomials.
2.4. Vector representation of the closed loop equation Using the notation introduced above, the closed loop characteristic equation is given by A(s, q, k) = kNe(s, q)Nc(s) + De(s, q)Dc(s).
=
for i = 0 , 1 . . . . . d e - 1 . Then, it is easily verified that the closed loop polynomial is given by A(s, q, k) rip
= s a~ + ~'~ (kaio + bio + kqi+~ + q,~÷i+z)s i i=0
dr,-- 1
(10)
*Strict properness is only assumed for notational convenience in the sequel. If instead, the plant is proper but not strictly proper, the results to follow only require minor modification; i.e. all the theory remains intact. ~"The development to follow readily extends to multi-input single-output systems and to single-input multi-output systems.
+
~
(bio+ q~+i+E)S i.
(13)
i=np+ l
Hence, we see that A(s, q, k) is an interval polynomial; i.e. we can write dr,--1
A(s, q, k) = s d~ + ~ i=0
t~i(q, k)s i
(14)
286
B . R . BARMISH and R. TEMPO
where
given by 6;(q, k) ~ [67(k), 67(k)]
for i = 0 , 1. . . . .
(15)
f2(to*, k*) = {z ~ C : R e Kl(/.to*, k*) -< Re z - Re Kz(/.~o*, k*);
n c - 1 and
Im g3(/.to*, k*) -.< Im z
kaio + bio + kq7+~ + q~e+i+2 I -~ f o r i = 0 , 1 . . . . . ne;
8- k ' ( ) - Jb~o + q~t.+i+2 / k fori=ne+l,n~,+2
-< Im K,(/.w*, k*)}. (16)
Observe that the sharpness of the bounds in (18) implies that all points in Q(to*,k*) are attainable by choice of q. Hence, it follows that s =jto* belongs to the cross section of the R R L at k* if and only if
..... de-1
and
=~
67(k)
kaio + bio + kqi++l + q.+~+i÷z . f o r i = 0 , 1 . . . . ,np; Ib,o+ q.+~+,+z
I
k
fori=ne+l,n~.+2
.....
0 ~ g2(to*, k*) (17)
Re Kl(jto*, k*) <- 0 <- Re Kz(jto*, k*)
Re K~(jto*, k*) <- Re A(j~o*, q, k*) -< Re K2(jto*, k*); Im K3(jto*, k*) -< Im A(jto*, q, k*) -< Im K4(jto*, k*)
(18)
where
K~(s, k ) = 6if(k) + 6?(k)s + , ~ ( k ) s ~ 3 -4- ~ ' ( k ) s
4 --~ ~ ( k ) s
5 -q-.. • ;
Kz(s, k) = 6~'(k) + 6 t ( k ) s + 6~(k)s z + 6~-(k)s 3 + 6 ; ( k ) s ~ + 6~-(k)s 5 + . . .
;
K3(s, k) = 6if(k) + 67(k)s + 6~'(k)s z + 6 f ( k ) s 3 + 6~'(k)s 4 + 6~-(k)s 5 + . . .
;
K4(s, k) = 6if(k) + 67(k)s + 6 f ( k ) s z + O~'(k)s 3 + b~'(k)s 4 + b~'(k)s 5 + . . . (19) are the four Kharitonov polynomials. Hence, defining the value set by fl(to*, k*) -" {A(/'to*, q, k*) :q;- _< q~ _< q~f o r i = l , 2 . . . . ,l},
(20)
it follows from (18) that ff2(to*, k*) is a rectangle
(22)
or, equivalently,
de-1.
We now motivate the construction of the cross section of the R R L at k = 1. In this regard, we use the well-known Kharitonov polynomials to determine which points on the imaginary axis lie on this cross section. The idea used to characterize these axis points is generalized to "non-axis points" and to "non-interval polynomials" in the next subsection. Indeed, for fixed frequency to*~ R and fixed gain k = k*, we want to decide if s =/'to* lies on the cross section of R R L at k*. To this end, first note that the theorem of Kharitonov (1978) leads to the sharpest possible bounds on A(jto*, q, k*). Namely, we have
-[- 6 ; ( k ) s
(21)
Im K3(jto*, k*) <--0 <- Im K4(/.to*, k*).
(23)
This leads to the following conclusion: For to*->0, the point s =/.to* belongs to the cross section of the R R L at k* if and only if (23) holds. It is easy to verify that for to*-<0, the identical condition holds with the labels reversed for Kl(/'to*, k*) and Kz(jto*, k*). Observation. An equivalent way to write (23) is H(to*, k*) ~ 0
(24)
where H(to, k) --+-max (Kl(jw, k), - K z ( j w , k),
K3(jto, k), - K 4 ( j w , k)}.
(25)
The way of checking the zero inclusion condition motivates a method for performing computations in the more general analysis to follow.
3.2. Key ideas for the general case Now we consider the more general case obtained when the plant coefficients depend affine linearly on q and we want to determine whether a given point s = z belongs to the RRL. Motivated by the interval polynomial analysis of the previous subsection, we fix k = k* and define the set f~(z, k*) "- {A(z, q, k * ) : q ~ Q}.
(26)
It now follows that z belongs to the cross section of the R R L at k* if and only if 0 ~ ~(z, k*).
(27)
Finally, to test zero inclusion condition above, we need a more concrete description of ~2(z, k*). To obtain such a description, we first note that an affine linear transformation T taking the rectangle Q into the complex plane has a convex polygonal image which is the convex hull
The Robust Root Locus of T applied to the extreme points qi; i.e.
T(Q) = c o n v {r(q~), T(q 2) . . . . , r(qZ)}.
(28)
Hence, if we consider T to be an affine linear mapping taking q into 6(q, k*), it follows that f~(z, k*) is the convex polygon described by £2(z, k*) = c o n v (A(z, qX, k*), A(z, q2, k*) . . . . .
A(z, qL, k*)}.
(29)
Book keeping simplification. Note that the zero inclusion condition (27) can be checked in many possible ways because fl(z, k*) is a convex polygon---even a brute force gfidding procedure will suffice. However, when performing digital computation, it is convenient to have a "compact formula" which can be used to check zero inclusion and can easily be coded. Recalling inequality (24) for interval polynomials with a static gain compensator, we now provide a similar result for this afl'ine linear case, i.e. we generate a function H(z, k) having the following property (proved in Appendix A): a point z belongs to the cross section of the R R L at k* if and only if H(z, k*) <- O.
function involves a sweep of the scalar variable p over the interval [0, 1]. 4. A M O R E
A simple example of an acceptable path is described by the boundary of the unit circle given by ~ r ( P ) = cos 2rrp + j sin 2~rp.
(32)
Another possibility for F is given in Barmish (1989): There, it is argued that computations are facilitated when F is a polyhedral set. For example, one reasonable choice is the unit "diamond" described by IResl + [Imsl = 1.
tr,. -" 0 for i > n o Finally, we define the compensator parameter matrices as follows: 0
0
0
0
0
0
0-
0
a~o 0
0
0
0
0
0~2 0
~q 0
0
0
0
O~ =" 0~3 0
a~2 0
0
Or"x
0
0
0
0
0
0
0
0
~ R (dc+d'')×2(d"+l) (37)
-0/300 0 0 /3x 0 /30 0 /32 0 /31 0~--" 0 /33 0 /32
numerical
generation
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
/3dc-1
0
0
0
0
0
1
0
/3d~--2
0 /3d¢-x
E R (dc+d~)x2(ae+l).
(38)
(34)
for i > ne and j = 0, 1 . . . . , l. Subsequently, we define the ~erturbation coefficient matrix as
Oab "~"
that
0
0
Step 3. Let
Notice
O~dc_1 0 Ogdc_2 0 O:d¢ 0 tree--1 0
and
(zx, z2) -" Re zx. Re z2 + Im Zl • Im z2. (35)
0
0
Similarly, to describe the plant perturbations, we let %'-0
where for zl, z2 e C, we use the inner product
H(z, k*) "- max h(z, k*, p).
F O R 6(q, k)
FORMULA
(33)
Step 2. Define the function h(z, k*, p) "- mi2n L ( ~ r ( P ) , iX(Z, qi, k*))
CONCRETE
The computations required in Section 3 can be further simplified by taking advantage of the underlying control system structure. In this section, we derive an expression for the closed loop vector 6(q, k) in terms of the control system parameters. The computation below essentially amounts to a minor modification of the idea used in Biernacki et al. (1987) to isolate the effect of the gain k on the vector 6(q, k). Note that this leads to computational savings because various matrices need only be computed once as k varies along the R R L . Indeed, for notational convenience, we define
(30)
We now provide a recipe for H(z, k*). Step 1. Pick any continuous path F in the complex plane which encircles the origin and let ~ r ( P ) denote a parameterization of this path which is obtained by varying p between 0 and 1; i.e. F -" ((I)r(p) : p ~ [0, 11). (31)
287
(36) of this
aol
ao2
• • • ao~-
box
bo~
-..
:
:
adel
ade2
0
0
bo~ :
R2(a~+Oxt
• . . ad~t
...
0 (39)
288
B . R . BARMISH and R. TEMPO
and the nominal plant coefficient vector as -aoo
and the desired bounding region is ~k--' {Z ~ C : ( R e z) 2 + ( I m z ) 2 < R ~ } .
"
boo O0 ~"
, . .
~ l
(40)
(2e+l).
ad~O -
1
-
The following lemma is obtained by a lengthy but straightforward computation. 4.L Lemma (See Appendix B for proof). Let O~,, O~, O,~ and 00 be defined as in (37), (38), (39) and (40), respectively. Then 6(q, k) = (kO~ + O~)(00 + O~bq).
(41)
(46)
Note that even greater computational savings can be achieved by using tighter bounds for the zeros of A(s, q, k); e.g. see Marden (1966). 6. DISCUSSION OF EXTENSIONS The assumption that the coefficients ai(q) of Nv(s, q) and bi(q) of De(s, q) depend afffine linearly on q clearly does not capture all "uncertainty structures" of interest. More generally, when ai(q) and bi(q) are polynomic functions of q, application of the R R L technique of this paper requires generating "overbounds" for the coefficients; e.g. suppose
Remark. To exploit the function H(z, k) in
a3(q) = q2t + q2q3 -- 4qi + 5
checking if a point z belongs to the R R L , one simply substitutes (41) into (11) which is then used in (34).
then to use the techniques in this paper, one can introduce new perturbations
fh = q~; 42 = q2~; 5. SIMPLIFICATION OF RRL COMPUTATION USING A CRUDE BOUND Note that if k = k* is fixed, we have a crude bound ~k. for the cross section of the R R L at k*. Hence, one can avoid testing zero inclusion condition for an unbounded set of z; i.e. computationally we need only generate a 2-dimensional gridding of @k" and check if zero belongs to ff~(z, k*) when z is a grid point. Indeed, for k-> 0, we seek a bound ~k on the zeroes of A(s, q, k); i.e.
43 =
q2q3
with bounds
14d-<4t; 1421-<42~; 1431<4~43 and replace a3(q) by d3(41, 42, 43) = 42 + 43 - 441 + 5. The key point to note, however, is that such a procedure leads to conservatism in the analysis; i.e. the resulting R R L obtained using an overbounding polynomial dp -
G(s, q) "-- ~, a,(q)s'
{z e C:A(z, q, k) = 0 for some q e Q} ~_ ~ .
(47)
i=1
(42) To generate an acceptable bound @~, we use the well known fact (see Marden, 1966) that for a monic n-th order polynomial p ( s ) = s " + ~'d~o~a~s~, all its zeros are interior to the circle of radius
R=l+max{la,[:i=O,
1. . . . .
n-l}.
(43)
Applying this result to equation (11), it follows that all zeros of A(s, q, k) lie interior to a circle of radius
Rk = 1 + max 16,(q j, k)l.
will be an overestimate of the true RRL. This motivates the following question: is there a class of "physically meaningful" uncertainty structures for the ai(q) and bi(q) for which the overbounding procedure is "tight" in a certain sense? The answer to this question is "yes." If the ai(q) and hi(q) are multilinear affine* in q, then stronger results can be obtained. Note that many so-called physically based parameterizations lead to multilinear affine functions for the ai(q) and bi(q) (see Dasguspta and Anderson (1987).
(44)
i.j
To complete the bounding process, we substitute for 6i(q j, k) in terms of the given plant and compensator coefficient functions. Indeed, letting e,- denote a unit vector in the i-th coordinate direction and substituting (41) into (44) we obtain
g k ----- 1 q- max lef(kO~, + Ot~)(0o + OabqJ)l i,j
(45)
6.1. The convex hull robust root locus for multilinear affine structures Continuing on the issue of tight R R L approximations for multilinear affine uncertainty * A function f : Q --* R is said to be multilinear affine if the following condition holds: if all c o m p o n e n t s ql . . . . . qt except one are fixed, then f is affine linear. F o r example, the function equation f (q) = 3qlq2q3 - 6 q t q 3 + 4q2q3 + 2q~ 2q2 + q3 - 1 is multilinear aliine.
The Robust Root Locus
6.2. Compensator tuning using the RRL and CH-RRL Analogous to classical root locus analysis, the R R L and C H - R R L can also be used to study variations in the closed loop poles with respect to control parameters other than the loop gain k. For example, consider the problem of "optimally locating" the pole a~ of a fixed compensator
structures, we introduce the following two definitions: a point z e 12 is said to lie on the Convex Hull Robust Root Locus (CH-RRL) if 0 e conv ffa(z, k)
(48)
for some k -> 0. For fixed k = k*, a point z e 12 is said to lie on the cross section of the C H - R R L at k* if (48) is satisfied with k = k*; i.e. 0 e conv ~(z, k*). Note that the R R L is a subset of the CH-RRL; i.e. if 0 e f2(s*, k), it follows that 0 e conv ff2(s*, k). Moreover, the C H - R R L degenerates to the R R L when the a,(q) and b~(q) are affine linear. For this special case, it is easy to show that conv fa(z, k) = fa(z, k)
C(s) = - - .
S--~
7. NUMERICAL EXAMPLE
In this section, we illustrate the application of the R R L technique using a model for the Concorde supersonic transport plane. The data for this system is found in Dorf (1980) and given in Fig. 1. In this configuration, the gain k of the rate gyroscope is used for stability augmentation. To this end, the specification on the closed loop system is to have the "dominant poles" having natural frequency
for all z e C and all k -> 0. To complete this subsection, it remains to characterize conv ~(z, k). To this end we exploit the so-called Mapping Theorem; e.g. see Zadeh and Desoer (1963). In a number of recent papers the Mapping Theorem has also proven to be quite useful in a robust stability context (e.g. see Saeki (1987); Sideris and de Gaston (1987) and de Gaston and Safonov (1988)). Stated in terms of the notation given here, the Mapping Theorem indicates that for fixed k = k* and fixed z e C, the convex hull of f~(z, k*) is a convex polygon in the complex plane generated from the extreme points q". More precisely,
(50)
(52)
at medium-weight cruise and y = 0.20 at light-weight cruise. For the parameters o),, and ¢,,, we study the effect of variations up to 10% around the nominal values ~o~, = 2.5;
~,, = 0.30.
Reformulation using q parameter notation. We make the following identifications between the
A i r c r a f t
Pit ch rate
,
lFilter
¢ = 0.707.
y = 0.02
Actuator
t
(Si)
Robustness issue. The issue of robustness arises because the aircraft parameters y, ton and ¢ may vary. In this example, we have
As in Section 3, to carry out C H - R R L computations, it is often convenient to use the function H(z, k) as defined in (36): for fixed k = k*, a point z e C belongs to the cross section of the R R L at k* if and only if H(z, k*) <_0. As was the case with the R R L , the computational burden associated with the C H - R R L is decreased by exploiting a bound ~ for the zeros of A(z, q, k). In this regard, it can readily be shown that the radial bound R~ given by (44) still holds.
r v o
a,n = 2 . S and damping
conv fa(z, k*) = c o n v {A(z, ql, k*),A(z, q2, k*),
S •
1
By generating an R R L (or a CH-RRL) with respect to or, we obtain a better handle on robust compensator design.
(49)
. . . . A(z, qZ., k*)}
289
I Rate
-
J gyro
FtG. 1. Block diagram for supersonic transport plane.
290
B . R . BARMISH and R. TEMPO
aircraft parameters and the uncertainties q~
Imaglnerg
y~0.11+qt;
5,'
2¢,,~o,, ~ 1.5 + q2; 2
o~,, -~ 6.25 + q3.
5.6
(53) 4.8
Given the ranges of variation above, we consider
4.0
3,2
q~ e [-0.09, 0.09]; q2 e [-0.285, 0.315];
2,4
(54)
Real
q3 E
[-1.1875, 1.3125].
10.2 -II 0
Next, we describe the plant in the qg notation. A lengthy but straightforward computation leads
9.0
JB.~
J..~ "
°2.4
-3.2
tO -4.0
3~(s, q)
-4.9
= 39.06 (0.176 + 1.6ql)s + (0.0176 + 0.16q0 s 2 + (1.5 + q2)s + (6.25 + q3) ;
-5,6
the compensator transfer function is C(s) = looo
25s 2 + 100s + 100 s ~ + 80. ls* + 3808s ~ + 56380s 2" + 255600s + 25000
Computational results. The computational results to follow were generated on an Apollo mainframe computer. First, by running an "ordinary root locus" for the nominal system (set q t = q2 = q3 = 0), we obtained the following: (i) The relevant range of rate gyro gains is 10-< 39.06. k -< 50. Equivalently, letting /~ -" k • 39.06,
-6.4
FIG.
2.
;ross
section of RRL for s = - 9 . 5 +j4.0.
~=390.6
around
"insensitive poles" in (ii) and (iii) above are not shown in the plot.
Choice of gyro gain k. Recalling the specifications for closed loop, it remains to select the gyro gain k. This is accomplished as follows: on the R R L , we superimpose lines emanating from zero at angles of 135° and 225 ° (representing the desired damping ratio) and a circle of radius 2.5 (representing the desired natural frequency). From the cross sections of the RRL, a reasonable compromise between these two specifications is achieved with/~ = 1953; i.e.
we investigate the behavior of the closed loop poles for
/~
k =
39.06
= 50.
390.6 -
(ii) Of the seven closed loop poles, one of the poles is at
2,0
2,6
s ~ -0.1 2,4
and is insensitive to changes in the gain/~. (iii) There is also a pair of closed loop poles at s --- - 3 0 . 9 +j39.5
and
2.2
2.0
s ~- - 3 0 . 9 - j 3 9 . 5 .
Real
As the gain/~ increases, these poles move further into the left half plane. It was also observed in some preliminary R R L runs that the poles at s~-0.1, s-~-30.9+j39.5 and s ~ - - 3 0 . 9 j39.5 were almost totally insensitive to variations in the q~. In the Figs. 2-5, we provide some highlights of the graphical output generated by the R R L computer program. These figures show a number of cross sections of the RRL. Note that the
-1.4
i
i
i
.1,3
-1.2
.11
i -1.o
i -o.9
~ -o8
t
~
-o.7
.05
~
-:~.2 .2.4 -2.$ .2.8
FIG.
3. Cross
section ol RRL for s = - 1 . 0 +]2.3.
~=390.6
around
The Robust Root Locus ~ Imaginary
4
~ 9.2
-~.6
I -8.8
I -8.4
I -8.0
I -7,S
T -72
!
, ~6.8
Real ~
.~
~0
-~2
section of RRL for /~=1953 s = -7.9 +j9.6.
around
~Imaginary 2,8
2.4
2.0
1.6
1.2
Real 2.2
I -20
I -1.8
I -1.6
I -1.4
I -1.2
I -1.0
I -0.8
I -0,6
B "
-12
-1 6
"2.0
-2,4 ,
-2~:
FIG. 5. Cross
one obtains (within the limits of numerical round-off) the exact distribution of the closed loop poles. At the level of multilinear uncertainties, however, some overbounding is required which is tight in the sense that we obtain the CH-RRL as an approximation for the true RRL. This suggests a fruitful area for future research--the generation of an RRL for more complicated uncertainty structure.
-6.4
.g
FIG. 4. Cross
291
section of RRL for /~=1953 s = - 1 . 6 ±jl.9.
around
8. CONCLUSION
The main technical novelty of this paper is the use of the zero inclusion condition to obtain the robust root locus. To date, such exclusion condition has been used largely in a robust stability context. As illustrated by the numerical example, the resulting computational technique is easy to use and readily lends itself to useful graphical display. Moreover, gridding of the Q set is avoided---computationally only a two dimensional gridding is used. A second important point to note is that for affine linear uncertainty structures, we obtain the "true" RRL. That is, at each value of the gain k,
REFERENCES Ackermann, J. E. (1980). Parameter space design of robust control systems. IEEE Trans. Aut. Control, AC-25, 1058-1072. Barmish, B. R. (1989). A generalization of Kharitonov's four polynomial concept for robust stability problems with linearly dependent coefficient perturbations. IEEE Trans. Aut. Control, AC-34, 157-165. Biemacki, R. M., H. Hwang and S. P. Bhattacharyya (1987). Robust stability with structured real parameter perturbations. IEEE Trans. on Aut. Control, AC-32, 495-506. Dasgupta, S. and B. D. O. Anderson (1988). Physically based parameterizations for designing adaptive controllers. Automatica, 23, 469-477. Dasgupta, S. (1988). Kharitonov's theorem revisited. Syst. Control Lett., 11, 381-384. de Gaston, R. R. E. and M. G. Safonov (1988). Exact calculation of the multiloop stability margin. IEEE Trans. Aut. Control, AC-33, 156-171. Doff, R. C. (1980). Modern Control Systems. AddisonWesley, Reading, MA. Faro, A. T. and J. S. Meditch (1978). A canonical parameter space for linear systems. IEEE Trans. Aut. Control, AC-23, 454-458. Fruchter, G., U. Srebro and E. Zeheb (1987). On several variable zero set and application to M I M e robust feedback stabilization. IEEE Trans. Circuits Syst. CAS-34, 1210-1220. Kharitonov, V. L. (1978). Asymptotic stability of an equilibrium position of a family of systems of linear differential equations. Differential'nye Uraveniya, 14, 1483-1485. Krali, A. M. (1970). The root locus method: a survey. SIAM Rev., 12, 64-72. Marden, M. (1966). Geometry of polynomials. Mathematical Surveys, no. 3, American Mathematical Society, Providence, RI. Minnichelli, R. J., J. J. Anagnost and C. A. Desoer (1989). An elementary proof of Kharitonov's theorem and extensions. IEEE Trans. Aut. Control, AC-34, 995-998. Saeki, M. (1986). A method of robust stability analysis with highly structured uncertainties. IEEE Trans. Aut. Control, AC-31, 935-940. Sideris, A. and R. R. E. de Gaston (1986). Multivariable stability margin calculation with uncertainty correlated parameters. Prec. 1EEE Conf. Decis. Control, Athens, Greece. Soh, Y. C., R. J. Evans, I. R. Petcrsen and R. E. Betz (1987). Robust pole assignment. Automatica, 23, pp. 601-610. Truxal, J. G. (1955). Automatic Feedback Control System Synthesis. McGraw-Hill, New York. Zadeh, L. A. and C. A. Desocr (1963). Linear System Theory. McGraw-Hill, New York. Zcheb, E. and E. Walach (1977). 2-parameter root-loci concept and some applications. IEEE Int. J. Circuit Theory Applic., 5, 305-315. Zchcb, E. and E. Walch (1981). Zero sets of multiparametcr functions and stability of multidimensional systems. IEEE Trans. Acoustic, Speech, Signal Processing, ASSP-29, 197-206.
292
B . R . BARMISH and R. TEMPO kO,, + O~ =
A P P E N D I X A Justification of inequality (30) Since ~(z,k*) is a convex polygon, zero inclusion condition holds if and only if there is no straight line separating the origin from the extreme points A(z, qq k*) of f~(z, k*). That is, given any z • C, there exists at least one i -< L such that Re z • Re A(s*, qi, k*) + Im z • Im A(s*, qZ, k*) _-<0. (55) Since the left-hand side above depends affine linearly on z, we do not have to check (55) for all z • C. We need only consider z along some continuous path in the complex plane which encircles the origin; i.e. condition zero inclusion condition holds if and only if for each p • [0, 1], there exists some i • {1, 2 . . . . . L} such that <~r(P), A(s*, ql, k*)) - 0 .
-
koq~ flo
0
0
koq fl~ koto [30
0
0
0
0
0
0
0
0
0
0
kot2 f12 koq [3~ ko~3 [33 k~2 f12 0
0
0
0
0
0
0
0
ketac-~ flac-~ kO~ac
1
kctac-'. [3dc-, kaac-~
flac-~ (60)
(56) and
Equivalently, letting
h(s*, p, k*) = rain ( ~ r ( P ) , A(s*, q~, k*)), i:~L
F a°(q) 1
(57)
Oo+O,,oq=~b°(q)[
(61)
it follows that 0 e f~(s*, k*) if and only if max h(s*, p, k*) <-O.
0~p:~l
(58) where
Hence, it follows that z belongs to the cross section of the RRL at k* if and only if
a,( q ) ± 0
H(z, k*) <-O.
for i > he. The proof is complete by simply computing the row by column product of (60) and (61).
(59)
Acknowledgements--The authors express gratitude to Mr A P P E N D I X B Proof of Lemma 4.1 Straightforward computations indicate that
Andres Takach for his work on the numerical computations. They are also grateful to Mr Richard Do~phus for his detailed comments on the manuscript.
Autoraatica, Vol, 26, No. 2, pp, 283-292, 1990 Printed in Great Britain.
The Robust Root Locus* B. R. BARMISH~'§ and R. TEMPO~tll
A simple technique is described for generating the root loci of a feedback system which includes perturbations entering affine linearly into the coefficients of the plant.
KeyWordsmRoot locus; perturbations; robustness; value set. Abstract--In this paper, a simple technique is described for generating the root loci of a feedback system which includes perturbations q ~ Rt entering affine linearly into the coefficients of the plant. Denoting the perturbed plant by ~(s, q) and the compensator by C(s), we address the following problem: given a bounding set Q = R~ for q, find the locus of points z in the complex plane such that 1 + k~(z, q)C(z)= 0 for some k >-0 and some q a Q. Such points z are said to lie on the Robust Root Locus. One of the strengths of the technique presented in this paper is that it avoids the "combinatoric explosion" synonymous with gddding the/-dimensional set Q and plotting a large number of ordinary root loci associated with gfidding the points. Instead, the technique given here exploits only a 2-dimensional gridding of a bounded subset of the complex plane. For the more general case when perturbations enter multilinearly into the coefficients, a certain convex hull robust root locus is constructed. Finally, the robust root locus technique is illustrated using a model of a supersonic transport plane.
of the plant m a y be perturbed. Hence, the " t r u e " root locus may differ from the root locus generated using the so-called nominal system. In this regard, the main focal point of this p a p e r is the generation of root loci with respect to variations in m o r e than one parameter. This work can be viewed as a systematic approach to robustness problems which were alluded to in the fifties and sixties; e.g. see Truxal (1955). In the later work of Z e h e b and Walach (1977), a two-parameter root locus p r o b l e m is considered and rather specific assumptions (motivated by circuit theory) are m a d e about the class of perturbations. A n u m b e r of papers following Z e h e b and Walach's work deal with the so-called zero set concept; e.g. see Z e h e b and Walach (1981) and Fruchter, et al. (1987). It is seen that t~he zero set provides a rather general f r a m e w o r k for dealing with multi-parameter root loci. That is, for each value of the gain k, a function can be given whose zero set corresponds with the cross section of the root locus at k. Hence, in principle, a multi-parameter root locus can be generated by computing these cross sections as k varies. In practice, however, there are difficulties associated with numerical computation of the zero set. By strengthening the hypotheses on the class of allowable perturbations, we obtain an enormous simplification in the characterization of multi-parameter root loci. To this end, attention in this p a p e r is restricted to the case when the perturbations enter the plant coefficients either linearly or multilinearly. This leads to the following type of result: a point z ~ C can be " t e s t e d " for m e m b e r s h i p in the root locus by simply checking whether or not the point s = 0 in the complex plane lies inside or outside a specially constructed 2-dimensional convex polygon. Hence, one avoids the " c o m binatoric explosion" in computation associated with an /-dimensional gridding of the rectangle
1. INTRODUCTION AND FORMULATION T o MOTIVATE this paper, recall the classical R o o t Locus Problem: given a nominal plant transfer function ~ ( s ) and a fixed c o m p e n s a t o r C(s), find the locus of points z e C such that 1 + k~(z)C(z) = 0
(1)
for some k-> 0. That is, find the locus of the closed loop pole locations as a function of the gain k; e.g. see Krall (1970) for a survey of classical root locus techniques. The take-off point for this p a p e r is the following fact: when performing robustness analysis and design, various physical p a r a m e t e r s * Received 5 July 1988; revised 24 April 1989; received in final form 14 August 1989. The original version of this paper was not presented at any IFAC meeting. This paper was recommended for publication in revised form by Associate Editor V. Ku~era under the direction of Editor H. Kwakernaak. ~"Department of Electrical and Computer Engineering, University of Wisconsin-Madison, Madison, WI 53706, U.S.A. , CENS-CNR, Politecnico di Torino, Corso Duea degli Abruzzi 24, 10129 Torino, Italy. Author to whom all correspondence should be addressed. § Supported by the National Foundation under Grant No. ECS-8612948 and NATO-CNR Fellowship Program. II Supported by funds of CENS-CNR of Italy. PdlTO ~ : 2 - 6
283
284
B . R . BARMISH and R. TEMPO
Q. That is, instead of a potentially high dimensional gridding of Q, one can test candidate z points in a two-dimensional subset of the complex plane. The type of "zero inclusion" condition above goes back a long way in the robust stability literature (see, e.g., Zadeh and Desoer, 1963) and has been revived in the last few years (e.g. see Saeki, 1986; de Gaston and Safonov, 1988; Dasgupta, 1988; Minnichelli et al., 1989). Furthermore, the zero exclusion condition is easy to check because only two dimensions are involved and the set to be checked is a convex polygon. The first distinguishing feature of this paper is the use of the zero inclusion condition within the context of root locus rather than robust stability as in previous work. The second distinguishing feature is that the technique given here readily lends itself to practical application and computer-aided analysis and design. To make this argument more compelling, we include a rather detailed example involving the tuning of a rate gyroscope using a model of a supersonic transport plane. To study variations in the root loci under multiple parameter variations, we let q e R ~ denote a vector of perturbations with i-th component qi satisfying q/- -< qi -< q,+
(2)
where the q,.+ and q~- are prescribed bounds. Hence, the vector of perturbations q is confined to the/-dimensional rectangle
Q={q:qT<-qi<-q~;i=l
. . . . ,l}.
(3)
To denote the dependence of the plant transfer function on q, we write ~(s, q). We now present the formal definition of the Robust Root Locus (RRL): a point z ~ C is said to lie on the Robust Root Locus (RRL) if
1 + k~(z, q)C(z) = 0
(4)
for some k - 0 and some q e Q. Since it is also important to know the distribution of the closed loop poles for a fixed value of the gain k, we provide the following definition: for fixed k = k*, a point z e C is said to lie on the cross section of the R R L at k* if
1 + k*~(z, q)C(z) = 0
(5)
for some q ~ Q. The main objective of this paper is to provide a computationally tractable method to generate accurately the R R L for the case when the coefficients of the plant are affine linear functions of q. To this end, we do not simply overbound the distribution of the poles in the s-plane--we find the distribution precisely. This
lack of conservatism may seem somewhat surprising because of the complicated relationship between the zeros of a polynomial and its coefficients (e.g. see Marden (1966)). In a sense, this paper can be viewed as a complement to the growing body of literature on parameter space methods (see e.g. Ackermann (1980)). One reason why parameter space methods are powerful is because it is much easier to study the effect of feedback on the coefficients of the characteristic polynomial than on the closed loop poles. One of the earliest papers to fully recognize and exploit this fact is that of Fam and Meditch (1978). This is also recently recognized in Soh et al. (1987) and Biernacki, et al. (1987). Noting, however, that the closed loop poles may be highly sensitive to the coefficients of the characteristic polynomial, a design which "looks good" from a parameter space point of view may have poles which are badly "smeared" over a large region in the complex plane. Hence, having the R R L available, one has the opportunity to check whether or not this phenomenon is occurring; i.e. examine the cross section of the R R L at k = 1. Subsequently, by tuning the loop gain, there is an opportunity to obtain a more desirable distribution of closed loop poles: said another way, the R R L can be used for robustness enhancement. The remainder of this paper is organized as follows. In Section 2, we introduce the necessary notation for the parameterization of the closed loop system with respect to q e Q. Section 3 provides a simple technique to decide whether or not a given point z ~ C lies on the RRL. In Section 4, the R R L technique is reinterpreted in terms of the control system parameters and in Section 5, a crude bound for the cross sections of the R R L is computed. In Section 6, we describe two extensions of the theory. We note that the R R L can be generated with respect to design variables other than a simple loop gain and we extend the R R L technique to the case when the perturbations enter multilinearly into the plant coefficients. In this multilinear case, we show that it possible to generate a so-called Convex Hull Robust Root Locus (CH-RRL). Finally, in Section 7 we present a numerical example and in Section 8, we provide conclusions.
2. BASIC NOTATION 2.1. The extreme points of the rectangle Q As stated in Section 1, we assume that the vector of perturbations q is confined to the /-dimensional rectangle Q. This rectangle has at most L -" 2 t extreme points denoted
The Robust Root Locus
q~, q2 . . . . . qL. The i-th extreme point has j-th component given by qj=q~
or
i__ q~ - q ~÷.
2.2. The plant parameterization For each fixed q e Q, the plant is a strictly proper* single-input single-output system~represented by the transfer function Np(s, q) e(s, q ) = De(s, q-~)
(6)
where Np(s, q) and De(s, q) are polynomials in s of the form
We(s, q) = a,,(q)s"" + a,,_l(q)s"'-' + " . +a~(q)s+ao(q); De(s, q) = s a" + ba~_~(q)s a~-~
(7)
+ . . . + b~(q)s + bo(q). Furthermore, the coefficients a~(q) and b~(q) above are given affine linear functions of the perturbation vector q; i.e. 1
ai(q) = a~o + ~ aiiq~ ./=1
285
We introduce the vector di(q, k ) e R a'+a~ to represent the coefficients of A ( s , q , k ) ; i.e. 6(q,k) has i-th component 6i_l(q,k) and hence, we write d~,+dc-I
A(s, q, k)=Sac÷d'+
~
5i(q, k)s i.
(11)
i=0
In Section 4, we express the vector 6(q, k) more directly in terms of the coefficients ai~, bit, o:~ and fig describing the plant ~ ( s , q ) and the compensator C(s). This type of expression proves useful for numerical computation.
3. KEY IDEAS BEHIND THE CONSTRUCTION OF THE RRL
As emphasized in the Introduction, it is often not feasible to generate an RRL by gridding the /-dimensional set Q. In this section, we see that the computation of the RRL can be achieved using only a 2-dimensional gridding of a bounded subset of the complex plane. This reduction in complexity leads to computational tractability which will be demonstrated using the numerical example presented in Section 7.
for i = 0, 1. . . . . ne and I
bi(q) = bio + Z biiql i=~ for i =0, 1. . . . .
de-1.
2.3. The compensator As usual, the compensator to be applied to the plant ~(s, q) is a proper rational transfer function
C(s) = k Nc(s) Dc(s)
(8)
3.1. Key ideas for the special case of an interval plant with a static gain To motivate the general construction of the RRL, we consider the special case of interval polynomials for Np(s, q) and De(s, q) and a static gain compensator C(s)
k.
(12)
Hence, suppose that each plant coefficient ai(q) or bi(q) depends on one and only one perturbation parameter; say
where k -> 0 is an adjustable gain and
ai(q) = aio + qi+~
Nc(s) = o~cs"C + o~c_~s rI~-~ + • • • + oqs + O~o
for i = 0, 1. . . . .
Dc(s) = sa" + fldc-I S d c - I + " ' " + fll S + ~0
np and
bi(q) = bio + qne+i+2
(9) are prescribed polynomials.
2.4. Vector representation of the closed loop equation Using the notation introduced above, the closed loop characteristic equation is given by A(s, q, k) = kNe(s, q)Nc(s) + De(s, q)Dc(s).
=
for i = 0 , 1 . . . . . d e - 1 . Then, it is easily verified that the closed loop polynomial is given by A(s, q, k) rip
= s a~ + ~'~ (kaio + bio + kqi+~ + q,~÷i+z)s i i=0
dr,-- 1
(10)
*Strict properness is only assumed for notational convenience in the sequel. If instead, the plant is proper but not strictly proper, the results to follow only require minor modification; i.e. all the theory remains intact. ~"The development to follow readily extends to multi-input single-output systems and to single-input multi-output systems.
+
~
(bio+ q~+i+E)S i.
(13)
i=np+ l
Hence, we see that A(s, q, k) is an interval polynomial; i.e. we can write dr,--1
A(s, q, k) = s d~ + ~ i=0
t~i(q, k)s i
(14)
286
B . R . BARMISH and R. TEMPO
where
given by 6;(q, k) ~ [67(k), 67(k)]
for i = 0 , 1. . . . .
(15)
f2(to*, k*) = {z ~ C : R e Kl(/.to*, k*) -< Re z - Re Kz(/.~o*, k*);
n c - 1 and
Im g3(/.to*, k*) -.< Im z
kaio + bio + kq7+~ + q~e+i+2 I -~ f o r i = 0 , 1 . . . . . ne;
8- k ' ( ) - Jb~o + q~t.+i+2 / k fori=ne+l,n~,+2
-< Im K,(/.w*, k*)}. (16)
Observe that the sharpness of the bounds in (18) implies that all points in Q(to*,k*) are attainable by choice of q. Hence, it follows that s =jto* belongs to the cross section of the R R L at k* if and only if
..... de-1
and
=~
67(k)
kaio + bio + kqi++l + q.+~+i÷z . f o r i = 0 , 1 . . . . ,np; Ib,o+ q.+~+,+z
I
k
fori=ne+l,n~.+2
.....
0 ~ g2(to*, k*) (17)
Re Kl(jto*, k*) <- 0 <- Re Kz(jto*, k*)
Re K~(jto*, k*) <- Re A(j~o*, q, k*) -< Re K2(jto*, k*); Im K3(jto*, k*) -< Im A(jto*, q, k*) -< Im K4(jto*, k*)
(18)
where
K~(s, k ) = 6if(k) + 6?(k)s + , ~ ( k ) s ~ 3 -4- ~ ' ( k ) s
4 --~ ~ ( k ) s
5 -q-.. • ;
Kz(s, k) = 6~'(k) + 6 t ( k ) s + 6~(k)s z + 6~-(k)s 3 + 6 ; ( k ) s ~ + 6~-(k)s 5 + . . .
;
K3(s, k) = 6if(k) + 67(k)s + 6~'(k)s z + 6 f ( k ) s 3 + 6~'(k)s 4 + 6~-(k)s 5 + . . .
;
K4(s, k) = 6if(k) + 67(k)s + 6 f ( k ) s z + O~'(k)s 3 + b~'(k)s 4 + b~'(k)s 5 + . . . (19) are the four Kharitonov polynomials. Hence, defining the value set by fl(to*, k*) -" {A(/'to*, q, k*) :q;- _< q~ _< q~f o r i = l , 2 . . . . ,l},
(20)
it follows from (18) that ff2(to*, k*) is a rectangle
(22)
or, equivalently,
de-1.
We now motivate the construction of the cross section of the R R L at k = 1. In this regard, we use the well-known Kharitonov polynomials to determine which points on the imaginary axis lie on this cross section. The idea used to characterize these axis points is generalized to "non-axis points" and to "non-interval polynomials" in the next subsection. Indeed, for fixed frequency to*~ R and fixed gain k = k*, we want to decide if s =/'to* lies on the cross section of R R L at k*. To this end, first note that the theorem of Kharitonov (1978) leads to the sharpest possible bounds on A(jto*, q, k*). Namely, we have
-[- 6 ; ( k ) s
(21)
Im K3(jto*, k*) <--0 <- Im K4(/.to*, k*).
(23)
This leads to the following conclusion: For to*->0, the point s =/.to* belongs to the cross section of the R R L at k* if and only if (23) holds. It is easy to verify that for to*-<0, the identical condition holds with the labels reversed for Kl(/'to*, k*) and Kz(jto*, k*). Observation. An equivalent way to write (23) is H(to*, k*) ~ 0
(24)
where H(to, k) --+-max (Kl(jw, k), - K z ( j w , k),
K3(jto, k), - K 4 ( j w , k)}.
(25)
The way of checking the zero inclusion condition motivates a method for performing computations in the more general analysis to follow.
3.2. Key ideas for the general case Now we consider the more general case obtained when the plant coefficients depend affine linearly on q and we want to determine whether a given point s = z belongs to the RRL. Motivated by the interval polynomial analysis of the previous subsection, we fix k = k* and define the set f~(z, k*) "- {A(z, q, k * ) : q ~ Q}.
(26)
It now follows that z belongs to the cross section of the R R L at k* if and only if 0 ~ ~(z, k*).
(27)
Finally, to test zero inclusion condition above, we need a more concrete description of ~2(z, k*). To obtain such a description, we first note that an affine linear transformation T taking the rectangle Q into the complex plane has a convex polygonal image which is the convex hull
The Robust Root Locus of T applied to the extreme points qi; i.e.
T(Q) = c o n v {r(q~), T(q 2) . . . . , r(qZ)}.
(28)
Hence, if we consider T to be an affine linear mapping taking q into 6(q, k*), it follows that f~(z, k*) is the convex polygon described by £2(z, k*) = c o n v (A(z, qX, k*), A(z, q2, k*) . . . . .
A(z, qL, k*)}.
(29)
Book keeping simplification. Note that the zero inclusion condition (27) can be checked in many possible ways because fl(z, k*) is a convex polygon---even a brute force gfidding procedure will suffice. However, when performing digital computation, it is convenient to have a "compact formula" which can be used to check zero inclusion and can easily be coded. Recalling inequality (24) for interval polynomials with a static gain compensator, we now provide a similar result for this afl'ine linear case, i.e. we generate a function H(z, k) having the following property (proved in Appendix A): a point z belongs to the cross section of the R R L at k* if and only if H(z, k*) <- O.
function involves a sweep of the scalar variable p over the interval [0, 1]. 4. A M O R E
A simple example of an acceptable path is described by the boundary of the unit circle given by ~ r ( P ) = cos 2rrp + j sin 2~rp.
(32)
Another possibility for F is given in Barmish (1989): There, it is argued that computations are facilitated when F is a polyhedral set. For example, one reasonable choice is the unit "diamond" described by IResl + [Imsl = 1.
tr,. -" 0 for i > n o Finally, we define the compensator parameter matrices as follows: 0
0
0
0
0
0
0-
0
a~o 0
0
0
0
0
0~2 0
~q 0
0
0
0
O~ =" 0~3 0
a~2 0
0
Or"x
0
0
0
0
0
0
0
0
~ R (dc+d'')×2(d"+l) (37)
-0/300 0 0 /3x 0 /30 0 /32 0 /31 0~--" 0 /33 0 /32
numerical
generation
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
/3dc-1
0
0
0
0
0
1
0
/3d~--2
0 /3d¢-x
E R (dc+d~)x2(ae+l).
(38)
(34)
for i > ne and j = 0, 1 . . . . , l. Subsequently, we define the ~erturbation coefficient matrix as
Oab "~"
that
0
0
Step 3. Let
Notice
O~dc_1 0 Ogdc_2 0 O:d¢ 0 tree--1 0
and
(zx, z2) -" Re zx. Re z2 + Im Zl • Im z2. (35)
0
0
Similarly, to describe the plant perturbations, we let %'-0
where for zl, z2 e C, we use the inner product
H(z, k*) "- max h(z, k*, p).
F O R 6(q, k)
FORMULA
(33)
Step 2. Define the function h(z, k*, p) "- mi2n L ( ~ r ( P ) , iX(Z, qi, k*))
CONCRETE
The computations required in Section 3 can be further simplified by taking advantage of the underlying control system structure. In this section, we derive an expression for the closed loop vector 6(q, k) in terms of the control system parameters. The computation below essentially amounts to a minor modification of the idea used in Biernacki et al. (1987) to isolate the effect of the gain k on the vector 6(q, k). Note that this leads to computational savings because various matrices need only be computed once as k varies along the R R L . Indeed, for notational convenience, we define
(30)
We now provide a recipe for H(z, k*). Step 1. Pick any continuous path F in the complex plane which encircles the origin and let ~ r ( P ) denote a parameterization of this path which is obtained by varying p between 0 and 1; i.e. F -" ((I)r(p) : p ~ [0, 11). (31)
287
(36) of this
aol
ao2
• • • ao~-
box
bo~
-..
:
:
adel
ade2
0
0
bo~ :
R2(a~+Oxt
• . . ad~t
...
0 (39)
288
B . R . BARMISH and R. TEMPO
and the nominal plant coefficient vector as -aoo
and the desired bounding region is ~k--' {Z ~ C : ( R e z) 2 + ( I m z ) 2 < R ~ } .
"
boo O0 ~"
, . .
~ l
(40)
(2e+l).
ad~O -
1
-
The following lemma is obtained by a lengthy but straightforward computation. 4.L Lemma (See Appendix B for proof). Let O~,, O~, O,~ and 00 be defined as in (37), (38), (39) and (40), respectively. Then 6(q, k) = (kO~ + O~)(00 + O~bq).
(41)
(46)
Note that even greater computational savings can be achieved by using tighter bounds for the zeros of A(s, q, k); e.g. see Marden (1966). 6. DISCUSSION OF EXTENSIONS The assumption that the coefficients ai(q) of Nv(s, q) and bi(q) of De(s, q) depend afffine linearly on q clearly does not capture all "uncertainty structures" of interest. More generally, when ai(q) and bi(q) are polynomic functions of q, application of the R R L technique of this paper requires generating "overbounds" for the coefficients; e.g. suppose
Remark. To exploit the function H(z, k) in
a3(q) = q2t + q2q3 -- 4qi + 5
checking if a point z belongs to the R R L , one simply substitutes (41) into (11) which is then used in (34).
then to use the techniques in this paper, one can introduce new perturbations
fh = q~; 42 = q2~; 5. SIMPLIFICATION OF RRL COMPUTATION USING A CRUDE BOUND Note that if k = k* is fixed, we have a crude bound ~k. for the cross section of the R R L at k*. Hence, one can avoid testing zero inclusion condition for an unbounded set of z; i.e. computationally we need only generate a 2-dimensional gridding of @k" and check if zero belongs to ff~(z, k*) when z is a grid point. Indeed, for k-> 0, we seek a bound ~k on the zeroes of A(s, q, k); i.e.
43 =
q2q3
with bounds
14d-<4t; 1421-<42~; 1431<4~43 and replace a3(q) by d3(41, 42, 43) = 42 + 43 - 441 + 5. The key point to note, however, is that such a procedure leads to conservatism in the analysis; i.e. the resulting R R L obtained using an overbounding polynomial dp -
G(s, q) "-- ~, a,(q)s'
{z e C:A(z, q, k) = 0 for some q e Q} ~_ ~ .
(47)
i=1
(42) To generate an acceptable bound @~, we use the well known fact (see Marden, 1966) that for a monic n-th order polynomial p ( s ) = s " + ~'d~o~a~s~, all its zeros are interior to the circle of radius
R=l+max{la,[:i=O,
1. . . . .
n-l}.
(43)
Applying this result to equation (11), it follows that all zeros of A(s, q, k) lie interior to a circle of radius
Rk = 1 + max 16,(q j, k)l.
will be an overestimate of the true RRL. This motivates the following question: is there a class of "physically meaningful" uncertainty structures for the ai(q) and bi(q) for which the overbounding procedure is "tight" in a certain sense? The answer to this question is "yes." If the ai(q) and hi(q) are multilinear affine* in q, then stronger results can be obtained. Note that many so-called physically based parameterizations lead to multilinear affine functions for the ai(q) and bi(q) (see Dasguspta and Anderson (1987).
(44)
i.j
To complete the bounding process, we substitute for 6i(q j, k) in terms of the given plant and compensator coefficient functions. Indeed, letting e,- denote a unit vector in the i-th coordinate direction and substituting (41) into (44) we obtain
g k ----- 1 q- max lef(kO~, + Ot~)(0o + OabqJ)l i,j
(45)
6.1. The convex hull robust root locus for multilinear affine structures Continuing on the issue of tight R R L approximations for multilinear affine uncertainty * A function f : Q --* R is said to be multilinear affine if the following condition holds: if all c o m p o n e n t s ql . . . . . qt except one are fixed, then f is affine linear. F o r example, the function equation f (q) = 3qlq2q3 - 6 q t q 3 + 4q2q3 + 2q~ 2q2 + q3 - 1 is multilinear aliine.
The Robust Root Locus
6.2. Compensator tuning using the RRL and CH-RRL Analogous to classical root locus analysis, the R R L and C H - R R L can also be used to study variations in the closed loop poles with respect to control parameters other than the loop gain k. For example, consider the problem of "optimally locating" the pole a~ of a fixed compensator
structures, we introduce the following two definitions: a point z e 12 is said to lie on the Convex Hull Robust Root Locus (CH-RRL) if 0 e conv ffa(z, k)
(48)
for some k -> 0. For fixed k = k*, a point z e 12 is said to lie on the cross section of the C H - R R L at k* if (48) is satisfied with k = k*; i.e. 0 e conv ~(z, k*). Note that the R R L is a subset of the CH-RRL; i.e. if 0 e f2(s*, k), it follows that 0 e conv ff2(s*, k). Moreover, the C H - R R L degenerates to the R R L when the a,(q) and b~(q) are affine linear. For this special case, it is easy to show that conv fa(z, k) = fa(z, k)
C(s) = - - .
S--~
7. NUMERICAL EXAMPLE
In this section, we illustrate the application of the R R L technique using a model for the Concorde supersonic transport plane. The data for this system is found in Dorf (1980) and given in Fig. 1. In this configuration, the gain k of the rate gyroscope is used for stability augmentation. To this end, the specification on the closed loop system is to have the "dominant poles" having natural frequency
for all z e C and all k -> 0. To complete this subsection, it remains to characterize conv ~(z, k). To this end we exploit the so-called Mapping Theorem; e.g. see Zadeh and Desoer (1963). In a number of recent papers the Mapping Theorem has also proven to be quite useful in a robust stability context (e.g. see Saeki (1987); Sideris and de Gaston (1987) and de Gaston and Safonov (1988)). Stated in terms of the notation given here, the Mapping Theorem indicates that for fixed k = k* and fixed z e C, the convex hull of f~(z, k*) is a convex polygon in the complex plane generated from the extreme points q". More precisely,
(50)
(52)
at medium-weight cruise and y = 0.20 at light-weight cruise. For the parameters o),, and ¢,,, we study the effect of variations up to 10% around the nominal values ~o~, = 2.5;
~,, = 0.30.
Reformulation using q parameter notation. We make the following identifications between the
A i r c r a f t
Pit ch rate
,
lFilter
¢ = 0.707.
y = 0.02
Actuator
t
(Si)
Robustness issue. The issue of robustness arises because the aircraft parameters y, ton and ¢ may vary. In this example, we have
As in Section 3, to carry out C H - R R L computations, it is often convenient to use the function H(z, k) as defined in (36): for fixed k = k*, a point z e C belongs to the cross section of the R R L at k* if and only if H(z, k*) <_0. As was the case with the R R L , the computational burden associated with the C H - R R L is decreased by exploiting a bound ~ for the zeros of A(z, q, k). In this regard, it can readily be shown that the radial bound R~ given by (44) still holds.
r v o
a,n = 2 . S and damping
conv fa(z, k*) = c o n v {A(z, ql, k*),A(z, q2, k*),
S •
1
By generating an R R L (or a CH-RRL) with respect to or, we obtain a better handle on robust compensator design.
(49)
. . . . A(z, qZ., k*)}
289
I Rate
-
J gyro
FtG. 1. Block diagram for supersonic transport plane.
290
B . R . BARMISH and R. TEMPO
aircraft parameters and the uncertainties q~
Imaglnerg
y~0.11+qt;
5,'
2¢,,~o,, ~ 1.5 + q2; 2
o~,, -~ 6.25 + q3.
5.6
(53) 4.8
Given the ranges of variation above, we consider
4.0
3,2
q~ e [-0.09, 0.09]; q2 e [-0.285, 0.315];
2,4
(54)
Real
q3 E
[-1.1875, 1.3125].
10.2 -II 0
Next, we describe the plant in the qg notation. A lengthy but straightforward computation leads
9.0
JB.~
J..~ "
°2.4
-3.2
tO -4.0
3~(s, q)
-4.9
= 39.06 (0.176 + 1.6ql)s + (0.0176 + 0.16q0 s 2 + (1.5 + q2)s + (6.25 + q3) ;
-5,6
the compensator transfer function is C(s) = looo
25s 2 + 100s + 100 s ~ + 80. ls* + 3808s ~ + 56380s 2" + 255600s + 25000
Computational results. The computational results to follow were generated on an Apollo mainframe computer. First, by running an "ordinary root locus" for the nominal system (set q t = q2 = q3 = 0), we obtained the following: (i) The relevant range of rate gyro gains is 10-< 39.06. k -< 50. Equivalently, letting /~ -" k • 39.06,
-6.4
FIG.
2.
;ross
section of RRL for s = - 9 . 5 +j4.0.
~=390.6
around
"insensitive poles" in (ii) and (iii) above are not shown in the plot.
Choice of gyro gain k. Recalling the specifications for closed loop, it remains to select the gyro gain k. This is accomplished as follows: on the R R L , we superimpose lines emanating from zero at angles of 135° and 225 ° (representing the desired damping ratio) and a circle of radius 2.5 (representing the desired natural frequency). From the cross sections of the RRL, a reasonable compromise between these two specifications is achieved with/~ = 1953; i.e.
we investigate the behavior of the closed loop poles for
/~
k =
39.06
= 50.
390.6 -
(ii) Of the seven closed loop poles, one of the poles is at
2,0
2,6
s ~ -0.1 2,4
and is insensitive to changes in the gain/~. (iii) There is also a pair of closed loop poles at s --- - 3 0 . 9 +j39.5
and
2.2
2.0
s ~- - 3 0 . 9 - j 3 9 . 5 .
Real
As the gain/~ increases, these poles move further into the left half plane. It was also observed in some preliminary R R L runs that the poles at s~-0.1, s-~-30.9+j39.5 and s ~ - - 3 0 . 9 j39.5 were almost totally insensitive to variations in the q~. In the Figs. 2-5, we provide some highlights of the graphical output generated by the R R L computer program. These figures show a number of cross sections of the RRL. Note that the
-1.4
i
i
i
.1,3
-1.2
.11
i -1.o
i -o.9
~ -o8
t
~
-o.7
.05
~
-:~.2 .2.4 -2.$ .2.8
FIG.
3. Cross
section ol RRL for s = - 1 . 0 +]2.3.
~=390.6
around
The Robust Root Locus ~ Imaginary
4
~ 9.2
-~.6
I -8.8
I -8.4
I -8.0
I -7,S
T -72
!
, ~6.8
Real ~
.~
~0
-~2
section of RRL for /~=1953 s = -7.9 +j9.6.
around
~Imaginary 2,8
2.4
2.0
1.6
1.2
Real 2.2
I -20
I -1.8
I -1.6
I -1.4
I -1.2
I -1.0
I -0.8
I -0,6
B "
-12
-1 6
"2.0
-2,4 ,
-2~:
FIG. 5. Cross
one obtains (within the limits of numerical round-off) the exact distribution of the closed loop poles. At the level of multilinear uncertainties, however, some overbounding is required which is tight in the sense that we obtain the CH-RRL as an approximation for the true RRL. This suggests a fruitful area for future research--the generation of an RRL for more complicated uncertainty structure.
-6.4
.g
FIG. 4. Cross
291
section of RRL for /~=1953 s = - 1 . 6 ±jl.9.
around
8. CONCLUSION
The main technical novelty of this paper is the use of the zero inclusion condition to obtain the robust root locus. To date, such exclusion condition has been used largely in a robust stability context. As illustrated by the numerical example, the resulting computational technique is easy to use and readily lends itself to useful graphical display. Moreover, gridding of the Q set is avoided---computationally only a two dimensional gridding is used. A second important point to note is that for affine linear uncertainty structures, we obtain the "true" RRL. That is, at each value of the gain k,
REFERENCES Ackermann, J. E. (1980). Parameter space design of robust control systems. IEEE Trans. Aut. Control, AC-25, 1058-1072. Barmish, B. R. (1989). A generalization of Kharitonov's four polynomial concept for robust stability problems with linearly dependent coefficient perturbations. IEEE Trans. Aut. Control, AC-34, 157-165. Biemacki, R. M., H. Hwang and S. P. Bhattacharyya (1987). Robust stability with structured real parameter perturbations. IEEE Trans. on Aut. Control, AC-32, 495-506. Dasgupta, S. and B. D. O. Anderson (1988). Physically based parameterizations for designing adaptive controllers. Automatica, 23, 469-477. Dasgupta, S. (1988). Kharitonov's theorem revisited. Syst. Control Lett., 11, 381-384. de Gaston, R. R. E. and M. G. Safonov (1988). Exact calculation of the multiloop stability margin. IEEE Trans. Aut. Control, AC-33, 156-171. Doff, R. C. (1980). Modern Control Systems. AddisonWesley, Reading, MA. Faro, A. T. and J. S. Meditch (1978). A canonical parameter space for linear systems. IEEE Trans. Aut. Control, AC-23, 454-458. Fruchter, G., U. Srebro and E. Zeheb (1987). On several variable zero set and application to M I M e robust feedback stabilization. IEEE Trans. Circuits Syst. CAS-34, 1210-1220. Kharitonov, V. L. (1978). Asymptotic stability of an equilibrium position of a family of systems of linear differential equations. Differential'nye Uraveniya, 14, 1483-1485. Krali, A. M. (1970). The root locus method: a survey. SIAM Rev., 12, 64-72. Marden, M. (1966). Geometry of polynomials. Mathematical Surveys, no. 3, American Mathematical Society, Providence, RI. Minnichelli, R. J., J. J. Anagnost and C. A. Desoer (1989). An elementary proof of Kharitonov's theorem and extensions. IEEE Trans. Aut. Control, AC-34, 995-998. Saeki, M. (1986). A method of robust stability analysis with highly structured uncertainties. IEEE Trans. Aut. Control, AC-31, 935-940. Sideris, A. and R. R. E. de Gaston (1986). Multivariable stability margin calculation with uncertainty correlated parameters. Prec. 1EEE Conf. Decis. Control, Athens, Greece. Soh, Y. C., R. J. Evans, I. R. Petcrsen and R. E. Betz (1987). Robust pole assignment. Automatica, 23, pp. 601-610. Truxal, J. G. (1955). Automatic Feedback Control System Synthesis. McGraw-Hill, New York. Zadeh, L. A. and C. A. Desocr (1963). Linear System Theory. McGraw-Hill, New York. Zcheb, E. and E. Walach (1977). 2-parameter root-loci concept and some applications. IEEE Int. J. Circuit Theory Applic., 5, 305-315. Zchcb, E. and E. Walch (1981). Zero sets of multiparametcr functions and stability of multidimensional systems. IEEE Trans. Acoustic, Speech, Signal Processing, ASSP-29, 197-206.
292
B . R . BARMISH and R. TEMPO kO,, + O~ =
A P P E N D I X A Justification of inequality (30) Since ~(z,k*) is a convex polygon, zero inclusion condition holds if and only if there is no straight line separating the origin from the extreme points A(z, qq k*) of f~(z, k*). That is, given any z • C, there exists at least one i -< L such that Re z • Re A(s*, qi, k*) + Im z • Im A(s*, qZ, k*) _-<0. (55) Since the left-hand side above depends affine linearly on z, we do not have to check (55) for all z • C. We need only consider z along some continuous path in the complex plane which encircles the origin; i.e. condition zero inclusion condition holds if and only if for each p • [0, 1], there exists some i • {1, 2 . . . . . L} such that <~r(P), A(s*, ql, k*)) - 0 .
-
koq~ flo
0
0
koq fl~ koto [30
0
0
0
0
0
0
0
0
0
0
kot2 f12 koq [3~ ko~3 [33 k~2 f12 0
0
0
0
0
0
0
0
ketac-~ flac-~ kO~ac
1
kctac-'. [3dc-, kaac-~
flac-~ (60)
(56) and
Equivalently, letting
h(s*, p, k*) = rain ( ~ r ( P ) , A(s*, q~, k*)), i:~L
F a°(q) 1
(57)
Oo+O,,oq=~b°(q)[
(61)
it follows that 0 e f~(s*, k*) if and only if max h(s*, p, k*) <-O.
0~p:~l
(58) where
Hence, it follows that z belongs to the cross section of the RRL at k* if and only if
a,( q ) ± 0
H(z, k*) <-O.
for i > he. The proof is complete by simply computing the row by column product of (60) and (61).
(59)
Acknowledgements--The authors express gratitude to Mr A P P E N D I X B Proof of Lemma 4.1 Straightforward computations indicate that
Andres Takach for his work on the numerical computations. They are also grateful to Mr Richard Do~phus for his detailed comments on the manuscript.