Short Communication
An interval method for computing the stability margin
of real uncertainty problems
Zheng-Ming Ge* R and Li-Wei ChuS
Department of Mechanical Engineering, National Chiao Tung University, 1001 Ta Hsueh Road, Hsinchu, Taiwan 30010, R.O.C.
SUMMARY
Frequently, in practical control system design, some designing parameters are uncertain. These uncertain parameters may vary with temperature, humidity or other environmental variable, and these variations will have an impact on the stability of the system. In this paper, we use a global optimal method}interval method, by which stability margins of these uncertain parameters can be computed. Copyright 2000 John Wiley & Sons, Ltd.
KEY WORDS: Stability margin;k-analysis; real uncertainty; interval arithmetic; global optimization
1. INTRODUCTION
The original idea of interval analysis was to bound rounding errors. However, interval mathema-tics can be said to have begun with the appearance of Moore's book [1] Interval Analysis in 1966. Moore's work transformed this simple idea into a useful tool for error analysis. Since the appearance of Moore's book, several persons have used interval analysis to solve the global optimization problem and systems of non-linear equations [2, 3]. Thus, the interval analysis has tended gradually to become an important mathematical tool for solving the problems of global optimization and systems of nonlinear equations. And we can also use interval analysis in robust control-system design.
Consider a robustness stability problem associated with real uncertain parameters. A funda-mental problem addressed in a large number of papers [4, 5] is: Determine the maximum uncertainty bound at which a system is stable. Our main technical objective in this paper is to
* Correspondence to: Z.-M. Ge, Department of Mechanical Engineering, National Chiao Tung University, 1001 Ta
Hsuch Road, Hsinchu, Taiwan 30010, R.O.C. R Professor.
S Graduate student.
Figure 1. A closed-loop feedback system.
show how interval analysis can be used to determine the stability margin of a robustness stability problem associated with real uncertain parameters. Although structured singular valuek-analysis can solve the type of problem, the stability margin yielded byk-analysis is more conservative than that obtained by using the interval method; see Section 2 for a detailed description.
2. PROBLEM FORMULATION
Consider a robustness stability problem as shown in Figure 1, where M(s)3RHL;L
, and the real uncertain part*31L;L. We know that this is a structured singular value problem. If we claim the closed-loop system is internally stable then the necessary and su$cient condition is
det(I!M( jw)*)O0 for w3[0,R).
Hence, if the uncertain part, #*#)(k*(M))\, then the system is always stable. But the condition#*#)(k*(M))\ is stringent. Here, we will use the interval method to obtain a loosen condition for the uncertain part* under which the system is still stable.
For convenience, we suppose M3"L;L and * is a diagonal matrix with real uncertain elements (if the block structure of* is not diagonal, we can transform the non-diagonal-structure problem into diagonal-structure problem by using certain techniques from Cheng and DeMoor [6]). Let *"+diag(dI,2, dLIL)"dG31,, and k*(M) is de"ned as
k*(M) :"
1
min+pN(*)"*3*; det(I!M*)"0,. Let*P is a subset of * de"ned as follows:
*P:"+*:*3*; "dG")r, i"1,2, n, where r31>
hence, we can say*P is a &&square'' in which "dG")r, i"1,2, n, and the length of &&square'' is r. At "rst we de"ned r:"sup+r: det(I!M*)O0,∀*3*P,. Then we will prove k*
(M)"1/r which mean the calculation of k-norm is to determine the maximum &square' within which det(I!M*)O0.
Let** is a solution of det(I!M*)"0, where**3* and**"diag(d* I,2, d*LIL). And let
pN(**)"min+pN(*)"*3*; det(I!M*)"0,. We know that pN(**)"max+"d*G": i"1,2, n,"r*, it follows thatk*(M)"1/r*.
Figure 2. The plot of f0"f'"0.
If r'r* then "d*G"(r(i"1,2, n). Hence,**3*P
and det(I!M**)O0, we have the
contradiction det(I!M**)"0. If r(r* then there exist r and* where r(r (r*, and * 3*PY such that det(I!M* )"0 and pN(* )"r (r*. It is contradicte with min+pN(*)"*3*; det(I!M*)"0,"pN(**)"r*. Therefore, r"r* and k*(M)"1/r.
From the above description, we know that the signi"cance ofk-analysis is to determine the maximum &square' of uncertain parameters within which det(I!M*)O0 (i.e. the system is stable). Let det(I!M*)"f0(d,2, dL)#f'(d,2, dL) j, where f0, f' are the real and imaginary
parts of det(I!M*) dependent on d,2, dL. If det(I!M*)O0, then f0(d,2, dL)O0 or
f'(d,2, dL)O0.
We can illustrate the geometrical signi"cance ofk*(M) by n"2, the geometrical signi"cance of
k*(M) is shown in Figure 2(a), which means f0O0 or f'O0 within the largest square on a plane
consisting ofd, d. Therefore, the structured singular value of k*"1/r.
From Figure 2(a), we know that if the uncertain part#*#)(k*(M))\"r, i.e. "d")r and
"d")r, then the system is always stable. But this claim is stringent. In fact, i.e. we let d and d fall within the rectangle shown in Figure 2(b), then the system is still stable. Therefore, the robustness stability problem can be transformed into a problem of systems of equations:
f0(d,2, dL)"0
f'(d,2, dL)"0 (M3"L
or f0(d,2, dL, w)"0 f'(d,2, dL, w)"0 (M(s)3RHL ;L ) (1b) wheredG3[dG, dG]; i"1,2, n; w3[0,R).
Our aim is to determine the maximum uncertainty bounddG at which the systems of equations (1) have no solutions. In the next section, we introduce a method for using interval analysis to solve this problem.
3. INTERVAL ALGORITHM
Consider a function f (x,2, xL), xG3XG (i"1,2, n), we expand f with Taylor's theorem [2].
f (y,2, yL)"f (x,2, xL)# L
GgG(f,2, fG, xG>,2, xL) (yG!xG)
where gG"(*f/*xG) (i"1,2, n). If xG3XG and yG3XG, then this holds for some number fG3XG. In our applications, we sometimes want a linear bound on f (y,2, yL) for all yG3XG (i"1,2, n). Thus, we replacefG with the bounding interval XG (i"1,2, n) and obtain
f (y,2, yL)3f (x,2, xL)# L
GgG (X,2, XG, xG>,2, xL) (yG!xG)
(2) If y is a zero of f, then f (y)"0 and we replace Equation (2) with
f (x)#g(x, X) (y!x)"0. (3)
We de"ne the solution set of Equation (3) to be S"+y: f (x)#g(x, f) (y!x)"0, for all f3X. This set contains any point y3X for which f (y)"0. From Equation (3), we let
>G"xG!f (x,2, xL)#gG(X,2, XG, xG>,2, xL) G\HgH) (yH!xH)# LHG>gH) (yH!xH) (i"1,2, n) (4) and the set Y"+(>,2, >L),. Thus, the set SLY, where the right-hand member of (4) is obtained by simple evaluation using interval arithmetic [3].
For future reference, it is desirable to have a distinctive notation for the solution of Equation (3). In place of >G and Y, we shall use the notation NG(x, X) and N(x, X), which emphasizes the dependence on X and x.
From Equation (3), we de"ne an iterative algorithm of the form
f (xI)#g(xI, XI) [N(xI, XI)!xI]"0 (5a)
XI>"XI5N(xI, XI) (5b)
for k"0, 1, 2,2, where xI is the centre of XI.
The components of N(xI, XI) will be computed sequentially. The intersection in (5b) should be performed as soon as a new component is obtained so that components computed later will be narrower intervals. And proof of convergence for Equation (5) can be found in Moore [1], Krawczky [7] and Alefeld [8].
Figure 3. The steps of interval method.
Now, we consider the robustness stability problem in Section 2. Ifd
G,dG, i"1,2, n, is given in Equation (1), we can use the iterative algorithm given in Equations (5) to solve Equation (1). If there are solutions in boxdG3[dG, dG](i"1,2, n),then wewilldecreasethesizeofbox,otherwise, we will increase the size of the box. We can tune the size of box successively until we get a maximum box from numerical computation such that no solutions exist in the box for Equation (1).
To illustrate the method described above, suppose n"2, due to this being a two-dimensional case, hence there are four variablesd
,d, d and d to be determined. First, let d"d"!r, d"d"r, as shown in Figure 3(a). If there are no solutions in box d3[d, d] and d3[d, d], then we increase the magnitude of r. Otherwise, we decrease the magnitude of r. Until r"dN*, as shown in Figure 3(b). Hold d"d* and let d"d"!r, d"r, then tune the magnitude of r until r"dN*, as shown in Figure 3(c). Then hold d"d* and tune the magnitude of "d
" and "d " again, until "d """d """dN*", as shown in Figure 3(d). Thus, we can decide the magnitude of d
, d, d and d and determine the maximum uncertain bounds for d and d. Simultaneously, to avoid the maximum uncertain bound values tend to in"nite, we will limit the maximum uncertain bound values less thana which is a very large number.
Finally, we have to note that M(s)3RHL;L
in the real uncertainty problem. Let M( jw)"[RG H(w)#IG H(w)j]; i, j"1,2, n; w3="[0,R), where RG H(w) and IG H(w) are the real
and imaginary parts of M( jw), and M( j=)"[RG H(=)#IG H(=)j] (i, j"1,2, n) is an interval matrix.
Although, we cannot use interval arithmetic to compute RG H(=) and IG H(=) (i, j"1,2, n) since ="[0,R] is an unbounded interval. We can restrict = within [0,wN ], where wN is a very large number or use the eigenvalue technique [9] to obtain the interval of RG H(=) and IG H(=) for ="[0,R].
We now describe the steps in the interval algorithm. The subroutine for the iterative algorithm (5) solving Equation (1) is omitted [9, 10].
Step 1: Input r,*r, m and a, let r
N G"rN G"r,*rN G"*rN G"*r; i"1,2, n.
Step 2: Let I"J"0. Step 3: Letd
G"!rN G, dG"rNG; i"1,2, n.
Step 4: Using the iterative algorithm (5a) and (5b), if there are solutions in the boxdG3[dG, dG]
(i"1,2, n), w3[0,R], then decrease the size of box: let I"1,*r
N G"(0.5)(*rN G and
*rN G"(0.5)(*rN G, rNG"rNG!*r
N Gand rN G"rNG!*rN G (i"1,2, n).
Step 5: Otherwise, increase the size of box, let J"1, *r
N G"(0.5)'*rN G and *rN G"(0.5)'*rNG,
r
N G"rN G#*rN Gand rN G"rNG#*rN G (i"1,2, n).
Step 6: If max(*r
N ,2,*rN L,*rN ,2,*rN L))m, then determine which rNH or rNH shouldbeheld; let the
r
N Hor rN H be held at rNH!m or rNH!m, and set the corresponding*rN Hor*rN H to zeros, and let the other*r
N G"*r and*rN G"*r (i"1,2, n); reset I"J"0.
Step 7: If every*r
N Gand*rN G (i"1,2, n) are zeros or max(rN,2, rNL, rN,2, rNL)*a, then stopand print outd
GanddG; i"1,2, n.
Step 8: Otherwise, go to step 3.
4. EXAMPLE
In this section, we give a comparison between using k-analysis and the interval method in determining the stability margins of a real uncertainty problem.
Assume the 5;5 transfer matrix M(s)" 1 s#3s#2 1 s#8s#17 1 s#3s#2 1 s#7s#12 1 s#3s#2 1 s#12s#61 s!7 s#9s#22 1 s#20s#96 1 s#18s#82 1 s#27s#170 s#1 s#6s#10s#8 1 s#8s#15 s!1 s#12s#11 s#5 s#6s#11s#6 s#4 s#11s#43s#65 1 s#14 1 s#16s#15 1 s#2s#10 1 s#7 1 s#9s#8 1 s#4s#29 s#1 s#2s#17 1 s#10s#74 1 s#3 1 s#3s#2 ,
and the corresponding real uncertainty matrices *" d 0 0 0 0 0 d 0 0 0 0 0 d 0 0 0 0 0 d 0 0 0 0 0 d , dG31, i"1,2, 4.
Becausek*(M(s))"0.8471, the closed-loop system consists of M(s) and* being stable in #*#)
(k*(M(s)))\"1.1804. So, the stability margin yielded by k-analysis is "d")1.1804,
"d")1.1804, "d")1.1804 and "d")1.1804. But, if we use the interval algorithm described in the preceding section to compute the stability margin of d, d, d and d. We can obtain the maximum uncertain bounds !3.1049)d)1.1804, !3.1049)d)1.1804, !3.1049)d) 1.1804 and !a)d)1.1804, where a"10000, within which the system is stable. So, the stability margin yielded by k-analysis is more conservative than that obtained by using the interval algorithm.
5. CONCLUSIONS
Recently, the following two important problems have attracted a lot of attention [11}13].
Problem 1 (Stability radii problem). For given matrices A3"L;L, B3"L;K, C3"N;L and a
nontrival partition of the complex plane"""E6"@, where "E is open region, and *3K is an unknown disturbance matrix belonging to a given perturbation set K measure the distance of the stable matrix A to instability, i.e. "nd
BK(A; B, C; "@):"inf+#*#; *3KK;N, p(A#B*C)5"@O, where#*# is any operator norm.
Problem 2. For the given stable matrices A, B and C "nd the largest interval matrix*' with
elements belonging to a given perturbation set K such that the interval matrix A(*')" A#B*'C is stable.
This two above problems can also be solved by interval method. The real stability radii problem (K"1) BK(A; B, C;"@) can be represented as follows:
BK(A; B, C;"@):"inf+#*#; *3KK;N, p(A#B*C)5"@O, "sup+r : p(A#B*C)5"@", ∀*3*P, where*P is de"ned by
*P:"+*:*3KK;N, #*#)r,, r31>
Let*"[dGH] (1)i)m, 1)j)p), and det(jI!A*(*)A(*))"D(*, j)"D(dG H, j)"D(dG H, p), wherep3"@51> and A(*):"A#B*C.
If the following two conditions hold:
(a) the intersection"@51> can be formulated in interval-form;
(b) ""*"")r can also be formulated as dG H3[dG H, dG H], where * are the real perturbation block, *"[dG H] (1)i)m, 1)j)p); then the interval method can solve real stability radii problem.
In addition, because of A, B and C are stable matrices for Problem 2, then the necessary and su$cient condition for A(*') to be stable is det(I!A*(*)A(*))O0, ∀*3*'. Let *"[dG H], dG H3[dG H, dG H], it follows that det(I!A*(*)A(*))"f0(dG H)#f'(dG H) j, where
f0 and f' are the real and imaginary parts of det(I!A*(*)A(*)) depend on dG H (1)i)m,
1)j)p).
Therefore, Problem 2 can be transformed into as a problem of systems of equations as Equation (1). We can use the interval method (illustrated in Figure 3) to determine the maximum uncertainty bounddG H at which det(I!A*(*)A(*))O0, ∀*3*'. Therefore, the largest interval matrix *' can be determined by the interval method within which A(*') is stable. The above results dealing with Problems 1 and 2 are obvious. The reader can easily get the results as the derivation in Section 2. Hereby, we will not derive them in further detail.
APPENDIX: NOMENCLATURE
det determinant
diag+ ) , diagonal matrix
I unit matrix
sup supremum
inf in"mum
Greek letters
* structured uncertainty matrix
k structured singular value with structured uncertainty matrix* ACKNOWLEDGEMENTS
This research was supported by the National Science Council, Republic of China, under Grant Number NSC86-2212-E-009-002.
REFERENCES 1. Moore RE. Interval Analysis, Prentice-Hall: Englewood Cli!s, NJ, 1966.
2. Hansen E. Global Optimization ;sing Interval Analysis, Marcel Dekker: New York, 1992.
3. Neumaier A. Interval Methods for Systems of Equations, Cambridge University Press: New York, 1990.
4. Barmish BR. New tools for robustness analysis. Proceedings of the IEEE Conference on Decision and Control, Austin, 1988; 1399}1404.
5. de Gaston RRE, Safonov MG. Exact calculation of the multiloop stability margin. IEEE ¹ransactions on Automatic
Control 1988; 33: 156}171.
6. Cheng Y, DeMoor B. A multidimensional realization algorithm for parametric uncertainty modelling and multi-parameter margin problem. International Journal of Control 1994; 60:789}807.
8. Alefeld G. On the convergence of some interval arithmetic modi"cations of Newton's method. SIAM Journal on
Numerical Analysis 1984; 21:363}372.
9. Ge ZM, Chu LW. An interval method for computing the structured singular value of a transfer matrix with real uncertainty. JSME International Journal, Series C, 1998; 41:741}750.
10. Moore RE. (ed.). Reliability in Computing, Academic Press: San Diego, 1988; 269}286.
11. Qiu L, Bernhardsson B, Rantzer A, Davison EJ, Young PM, Doyle JC. A formula for computation of the real structured stability radius. Automatica 1995; 31:879}890.
12. Son NK, Hinrichsen D. Robust stability of positive linear systems. Proceedings of the 34th Conference on Decision and
Control, New Orleans, 1994; 2985}2987.
