Fast Universal Algorithms for Robustness Analysis
aa r X i v : . [ m a t h . O C ] M a y Fast Universal Algorithms for Robustness Analysis ∗ Xinjia Chen, Kemin Zhou, and Jorge L. AravenaDepartment of Electrical and Computer EngineeringLouisiana State UniversityBaton Rouge, LA 70803 { chan,kemin,aravena } @ece.lsu.eduTel: (225)578- { } Fax: (225) 578-5200March 2003
Abstract
In this paper, we develop efficient randomized algorithms for estimating probabilistic robust-ness margin and constructing robustness degradation curve for uncertain dynamic systems. Oneremarkable feature of these algorithms is their universal applicability to robustness analysis prob-lems with arbitrary robustness requirements and uncertainty bounding set. In contrast to existingprobabilistic methods, our approach does not depend on the feasibility of computing deterministicrobustness margin. We have developed efficient methods such as probabilistic comparison, prob-abilistic bisection and backward iteration to facilitate the computation. In particular, confidenceinterval for binomial random variables has been frequently used in the estimation of probabilisticrobustness margin and in the accuracy evaluation of estimating robustness degradation function.Motivated by the importance of fast computing of binomial confidence interval in the contextof probabilistic robustness analysis, we have derived an explicit formula for constructing theconfidence interval of binomial parameter with guaranteed coverage probability. The formulaovercomes the limitation of normal approximation which is asymptotic in nature and thus in-evitably introduce unknown errors in applications. Moreover, the formula is extremely simpleand very tight in comparison with classic Clopper-Pearson’s approach.
In recent years, there have been growing interest on the development of probabilistic methods forrobustness analysis and design problems aimed at overcoming the computational complexity and ∗ This research was supported in part by grants from NASA (NCC5-573) and LEQSF (NASA /LEQSF(2001-04)-01). deterministic robustness margin which is defined as the maximal radius of uncertainty set suchthat the robustness requirement is guaranteed everywhere. However, it should be borne in mindthat the uncertainty set may include worst cases which never happen in reality. Instead of seekingthe worst case guarantee, it is sometimes “acceptable” that the robustness requirement is satisfiedfor most of the cases. It has been demonstrated that the proportion of systems guaranteeing therobustness requirement can be close to 1 even if the radii of uncertainty set are much larger than thedeterministic robustness margin. The idea of sacrificing extreme cases of uncertainty has become anew paradigm of robust control. In particular, the concept of probabilistic robustness margin has beenrecently established [2, 5, 24, 17, 4, 6, 23]. Moreover, to provide more insight for the robustness of theuncertain system, the concept of robustness degradation function has been developed by a numberof researchers [2, 6]. For example, Barmish and Lagoa [5] have constructed a curve of robustnessmargin amplification versus risk in a probabilistic setting. In a similar spirit, Calafiore, Dabbeneand Tempo [6, 7] have constructed a probability degradation function in the context of real andcomplex parametric uncertainty. Such a function describes quantitatively the relationship betweenthe proportion of systems guaranteeing the robustness requirement and the radius of uncertaintyset. Clearly, it can serve as a guide for control engineers in evaluating the robustness of a controlsystem once a controller design is completed. It is important to note that the robustness degradationfunction and the probabilistic robustness margin can be estimated in a distribution free manner. Thiscan be justified by the Truncation Theory established by Barmish, Lagoa and Tempo [2] and canalso be illustrated by relaxing the deterministic worst-case paradigm.Existing techniques for constructing the robustness degradation curve relies on the feasibility ofcomputing the deterministic robustness margin. Obtaining the probabilistic robustness margin ispossible only when the robustness degradation curve is constructed. However, the computation ofdeterministic robustness margin is possible only when the robustness requirement P is simple and theuncertainty bounding set is of some special structure. For example, P is stability requirement anduncertainty bounding sets are spectral normed balls. In that case, deterministic analysis methodssuch as µ -theory can be applied to compute the deterministic robustness margin. The determinis-tic robustness margin is then taken as a starting point of uncertainty radius interval for which therobustness degradation curve is to be constructed. In general, the problem of computing the deter-ministic robustness margin is not tractable. Hence, to construct a robustness degradation curve ofpractical interest, the selection of uncertainty radius interval is itself a question. Clearly, the range ofuncertainty radius for which robustness degradation curve is significantly below 1 is not of practicalinterest since only a small risk can be tolerated in reality. From application point of view, it is onlyneeded to construct robustness degradation curve for the range of uncertainty radius such that the2urve is above an a priori specified level 1 − ǫ where risk parameter ǫ ∈ (0 ,
1) is acceptably small.In this work, we consider robustness analysis problems with arbitrary robustness requirement anduncertainty bounding set. We develop efficient randomized algorithms for estimating probabilisticrobustness margin ρ ǫ which is defined as the maximal uncertainty radius such that the probabilityof guaranteeing the robust requirements is at least 1 − ǫ . We have also developed fast algorithms forconstructing robustness degradation curve which is above a priori specified level 1 − ǫ . In particular,we have developed efficient mechanisms such as probabilistic comparison, probabilistic bisection andbackward iteration to reduce the computational complexity. Complexity of probabilistic comparisontechniques are rigorously quantified.In our algorithms, confidence interval for binomial random variables has been frequently used toimprove the efficiency of estimating probabilistic robustness margin and in the accuracy evaluationof robustness degradation function. Obviously, fast construction of binomial confidence interval isimportant to the efficiency of the randomized algorithm. Therefore, we have derived an explicitformula for constructing the confidence interval of binomial parameter with guaranteed coverageprobability [11]. The formula overcomes the limitation of normal approximation which is asymptoticin nature and thus inevitably introduce unknown errors in applications. Moreover, the formula isextremely simple and very tight in comparison with classic Clopper-Pearson’s approach.The paper is organized as follows. Section 2 is the problem formulation. Section 2 discussesbinomial confidence interval. Section 3 is devoted to probabilistic robustness margin. Section 4presents algorithms for constructing robustness degradation curve. Illustrative examples are givenin Section 5. Section 6 is the conclusion. Proofs are included in Appendix. We adopt the assumption, from the classical robust control framework, that the uncertainty isdeterministic and bounded. We formulate a general robustness analysis problem in a similar way as[10] as follows.Let P denote a robustness requirement. The definition of P can be a fairly complicated combi-nation of the following: • Stability or D -stability; • H ∞ norm of the closed loop transfer function; • Time specifications such as overshoot, rise time, settling time and steady state error.Let B ( r ) denote the set of uncertainties with size smaller than r . In applications, we are usuallydealing with uncertainty sets such as the following:3 l p ball B p ( r ) := { ∆ ∈ R n : || ∆ || p ≤ r } where || . | p denotes the l p norm and p = 1 , , · · · , ∞ . Inparticular, B ∞ ( r ) denotes a box. • Spectral norm ball B σ ( r ) := { ∆ ∈ ∆ : ¯ σ (∆) ≤ r } where ¯ σ (∆) denotes the largest singularvalue of ∆. The class of allowable perturbations is ∆ := { blockdiag[ q I r , · · · , q s I r s , ∆ , · · · , ∆ c ] } (1)where q i ∈ F , i = 1 , · · · , s are scalar parameters with multiplicity r , · · · , r s and ∆ i ∈ F n i × m i , i = 1 , · · · , c are possibly repeated full blocks. Here F is either the complex field C or the real field R . • Homogeneous star-shaped bounding set B H ( r ) := { r (∆ − ∆ ) + ∆ : ∆ ∈ Q } where Q ⊂ R n and ∆ ∈ Q (see [2] for a detailed illustration).Throughout this paper, B ( r ) refers to any type of uncertainty set described above. Define a function ℓ ( . ) such that, for any X , ℓ ( X ) := min { r : X ∈ B ( r ) } , i.e., B ( ℓ ( X )) includes X exactly in the boundary. By such definition, ℓ ( X ) = min (cid:26) r : X − ∆ r + ∆ ∈ Q (cid:27) ,ℓ ( X ) = ¯ σ ( X ) , and ℓ ( X ) = || X || p in the context of homogeneous star-shaped bounding set, spectral norm ball and l p ball respectively.To allow the robustness analysis be performed in a distribution-free manner, we introduce thenotion of proportion as follows. For any ∆ ∈ B ( r ) there is an associated system G (∆). Define proportion P ( r ) := vol( { ∆ ∈ B ( r ) : The associated system G (∆) guarantees P } )vol( B ( r ))with vol( S ) := Z q ∈ S dq, where the notion of dq is illustrated as follows: • (I): If q = [ x rs ] n × m is a real matrix in R n × m , then dq = Q nr =1 Q ms =1 dx rs . • (II): If q = [ x rs + jy rs ] n × m is a complex matrix in C n × m , then dq = Q nr =1 Q ms =1 ( dx rs dy rs ).4 (III): If q ∈ ∆ , i.e., q possess a block structure defined by (1), then dq = ( Q si =1 dq i )( Q ci =1 d ∆ i )where the notion of dq i and d ∆ i is defined by (I) and (II).It follows that P ( r ) is a reasonable measure of the robustness of the system [6, 25]. In the worstcase deterministic framework, we are only interested in knowing if P is guaranteed for every ∆.However, one should bear in mind that the uncertainty set in our model may include worst caseswhich never happen in reality. Thus, it would be “acceptable” in many applications if the robustnessrequirement P is satisfied for most of the cases. Hence, due to the inaccuracy of the model, we shouldalso obtain the value of P ( r ) for uncertainty radius r which exceeds the deterministic robustnessmargin.Clearly, P ( r ) is deterministic in nature. However, we can resort to a probabilistic approach toevaluate P ( r ). To see this, one needs to observe that a random variable with uniform distribution over B ( r ), denoted by ∆ u , guarantees thatPr { ∆ u ∈ S } = vol( S T B ( r ))vol( B ( r ))for any S , and thus P ( r ) = Pr { The associated system G ( ∆ u ) guarantees P } . Define a Bernoulli random variable X such that X takes value 1 if the associated system G ( ∆ u )guarantees P and takes value 0 otherwise. Then estimating P ( r ) is equivalent to estimating binomialparameter P X := Pr { X = 1 } = P ( r ). It follows that a Monte Carlo method can be employed toestimate P ( r ) based on i.i.d. observations of X .Obviously, the robustness problem will be completely solved if we can efficiently estimate P ( r )for all r ∈ (0 , ∞ ). However, this is infeasible from computational perspective. In practice, only asmall risk ǫ can be tolerated by a system. Therefore, what is really important to know is the value of P ( r ) over the range of uncertain radius r for which P ( r ) is at least − ǫ where ǫ ∈ (0 ,
1) is referredas the risk parameter in this paper. When the deterministic robustness margin ρ := sup { r : G (∆) guarantees P for any ∆ ∈ B ( r ) } can be computed, we can construct the robustness degradation curve as follows: (I) Estimate pro-portion for a sequence of discrete values of uncertainty radius which is started from r = ρ andgradually increased; (II) The construction process is continued until the proportion is below 1 − ǫ .This is the conventional approach and it is feasible only when computing ρ is tractable. Unfor-tunately, in general there exists no effective technique for computing the deterministic robustnessmargin ρ . In light of this situation, we establish a new approach which does not depend on the5easibility of computing the deterministic robustness margin. Our strategy is to firstly estimate the probabilistic robustness margin ρ ( ǫ ) := sup { r : P ( r ) ≥ − ǫ } and consequently construct the robust degradation curve in a backward direction (in which r isdecreased) by choosing the estimate of ρ ( ǫ ) as the starting uncertainty radius.To reduce computational burden, the estimation of probabilistic robustness margin relies on thefrequent use of binomial confidence interval. The confidence interval is also served as a validationmethod for the accuracy of estimating robustness degradation function. Hence, it is desirable toquickly construct binomial confidence interval with guaranteed coverage probability. Clopper and Pearson [12] provided a rigorous approach for constructing binomial confidence interval.However, the computational complexity involved with this approach is very high. The standardtechnique is to use normal approximation which is not accurate for rare events. The coverageprobability of the confidence interval derived from normal approximation can be significantly belowthe specified confidence level even for very large sample size. In the context of robustness analysis,we are dealing with rare events because the probability that the robustness requirement is violatedis usually very small. We shall illustrate these standard methods as follows.
Let the sample size N and confidence parameter δ ∈ (0 ,
1) be fixed. We refer an observation of X with value 1 as a successful trial . Let K denote the number of successful trials during the N i.i.d.sampling experiments. Let k be a realization of K . The classic Clopper-Pearson lower confidencelimit L N,k,δ and upper confidence limit U N,k,δ are given respectively by L N,k,δ := ( k = 0 p if k > U N,k,δ := ( k = Np if k < N where p ∈ (0 ,
1) is the solution of the following equation k − X j =0 (cid:18) Nj (cid:19) p j (1 − p ) N − j = 1 − δ p ∈ (0 ,
1) is the solution of the following equation k X j =0 (cid:18) Nj (cid:19) p j (1 − p ) N − j = δ . (3)6 .2 Normal Approximation It is easy to see that the equations (2) and (3) are hard to solve and thus the confidence limits aredifficult to determine using Clopper-Pearson’s approach. For large sample size, it is computationallyintensive. To get around the difficulty, normal approximation has been widely used to develop simpleapproximate formulas (see, for example, [15] and the references therein). Let Φ( . ) denote the normaldistribution function and Z δ denote the critical value such that Φ( Z δ ) = 1 − δ . It follows from theCentral Limit Theorem that, for sufficiently large sample size N , the lower and upper confidencelimits can be estimated respectively as e L ≈ kN − Z δ q kN (1 − kN ) N and e U ≈ kN + Z δ q kN (1 − kN ) N .The critical problem with the normal approximation is that it is of asymptotic nature. It is notclear how large the sample size is sufficient for the approximation error to be negligible. Such anasymptotic approach is not good enough for studying the robustness of control systems. It is desirable to have a simple formula which is rigorous and very tight for the confidence intervalconstruction. Recently, we have derived the following simple formula for constructing the confidencelimits [11] (The proof is provided in Appendix for purpose of completeness).
Theorem 1
Let L ( N, k, δ ) = kN +
34 1 − kN − q θ k (1 − kN )1+ θN and U ( N, k, δ ) = kN +
34 1 − kN + q θ k (1 − kN )1+ θN with θ =
98 ln δ . Then Pr {L ( N, K, δ ) < P X < U ( N, K, δ ) } > − δ . Figures 1 and 2 show the confidence limits derived by different methods (curve A and B representrespectively the upper and lower confidence limits computed by Theorem 1; curve C and D representrespectively the upper and lower confidence limits calculated by Clopper-Pearson’s method). It canbe seen from these figures that our formula is very tight in comparison with the Clopper-Pearson’sapproach. Obviously, there is no comparison on the computational complexity. Our formula issimple enough for hand calculation. Simplicity of the confidence interval is especially important inthe context of our robustness analysis problem since the confidence limits are repeatedly used for alarge number of simulations.
In this section, we shall develop efficient randomized algorithms for constructing an estimate for ρ ( ǫ ).7
100 200 300 400 500 600 700 800 900 100000.10.20.30.40.50.60.70.80.91 Number of Successful Trials C on f i den c e L i m i t s A B C D
Figure 1: Confidence Interval ( N = 1000 , δ = 10 − ) C on f i den c e L i m i t s A B C D
Figure 2: Confidence Interval ( N = 1000 , δ = 0 . .1 Separable Assumption We assume that the robustness degradation curve of the system can be separated into two parts by ahorizontal line with height − ǫ , i.e., P ( r ) < − ǫ for all r ≥ ρ ( ǫ ) . We refer such an assumption as the
Separable Assumption . From an application point of view,this assumption is rather benign. Our extensive simulation experience indicated that, for small riskparameter ǫ , most control systems guarantee the separable assumption. It should be noted that it iseven much weaker than assuming that P ( r ) is non-decreasing (See illustrative Figure 3). Moreover,the non-increasing assumption is rather mild. This can be explained by a heuristic argument asfollows. Let B P ( r ) := { ∆ ∈ B ( r ) : The associated system G (∆) guarantees P } . Then P ( r ) = vol( B P ( r ))vol( B ( r ))and d P ( r ) dr = 1vol( B ( r )) (cid:20) d vol {B P ( r ) } dr − P ( r ) d vol {B ( r ) } dr (cid:21) . (4)Moreover, due to the constraint of robust requirement P , it is true that B P ( r + dr ) \ B P ( r ) is a subsetof B ( r + dr ) \ B ( r ) and it follows that vol {B ( r ) } increases (as r increases) faster than vol {B P ( r ) } ,i.e., d vol {B P ( r ) } dr ≤ d vol {B ( r ) } dr . Hence inequality d vol {B P ( r ) } dr ≤ P ( r ) d vol {B ( r ) } dr is not hard to guaranteein the range of r such that P ( r ) is close to 1. It follows from equation (4) that d P ( r ) dr ≤ ǫ .It is interesting to note that our randomized algorithms for estimating ρ ( ǫ ) presented in the sequelrely only on the following assumption: P ( r ) < − ǫ ∀ r ∈ ( ρ ( ǫ ) , κ ] [ { i : κ < i ≤ } (5)where κ = ⌈ log ρ ( ǫ ) ⌉ . It can be seen that condition (5) is even weaker than the separable assumption.When condition (5) is guaranteed, an interval which includes ρ ( ǫ ) can be readily found by startingfrom uncertainty radius r = 1 and then successively doubling r or cutting r in half based on thecomparison of P ( r ) with 1 − ǫ . Moreover, bisection method can be employed to refine the estimatefor ρ ( ǫ ). Of course, the success of such methods depends on the reliable and efficient comparisonof P ( r ) with 1 − ǫ based on Monte Carlo method. In the following subsection, we illustrate a fastmethod of comparison. 9 obustness(cid:13)Degradation(cid:13)Curve(cid:13) Proportion(cid:13)
Uncertainty Radius(cid:13)
Figure 3: Illustration of Separable Assumption. The robustness degradation curve can be separatedas the upper segment and lower segment by the dashed horizontal line with height 1 − ǫ . In thisexample, the separable assumption is satisfied, while the non-increasing assumption is violated. In general, P X can only be estimated by a Monte Carlo method. The conventional method is todirectly compare KN with 1 − ǫ where K is the number of successful trials during N i.i.d. samplingexperiments. There are three problems with the conventional method. First, the comparison of KN with 1 − ǫ can be very misleading. Second, the sample size N is required to be very large to obtaina reliable comparison. Third, we don’t know how reliable the comparison is. In this subsection,we present a new approach which allows for a reliable comparison with many fewer samples. Thekey idea is to compare binomial confidence limits with the fixed probability 1 − ǫ and hence reliablejudgement can be made in advance.Function name: Probabilistic-Comparison .Input: Risk parameter ǫ and confidence parameter δ .Output: d = Probabilistic-Comparison ( δ, ǫ ).Step 1. Let d ← d = 0 do the following: • Sample X . 10 −3 −2 −1 X S a m p l e S i z e Figure 4: Complexity of Probabilistic Comparison. The horizontal axis represents 1 − P X . Thevertical axis represents sample size. The solid line and the dash-dot line respectively show theaverage sample size and the 95%-quantile of the sample size. • Update N and K . • Compute lower confidence limit L and and upper confidence limit U by Theorem 1. • If U < − ǫ then let d ← −
1. If
L > − ǫ then let d ← δ is used to control the reliability of the comparison. A typical valueis δ = 0 .
01. The implication of output is as follows: d = 1 indicates that P X > − ǫ is true withhigh confidence; d = − P X < − ǫ is true with high confidence.Obviously, the sample size is random in nature. For ǫ = δ = 0 .
01, we simulated the ProbabilisticComparison Algorithm identically and independently 100 times for different values of P X . We observethat, for each value of P X , the Probabilistic Comparison Algorithm makes correct judgement amongall 100 simulations. Figure 4 shows the average sample size and the 95%-quantile of the sample sizeestimated from our simulation. It can be seen from the figure that, as long as P X is not very close to1 − ǫ , the Probabilistic Comparison Algorithm can make a reliable comparison with a small samplesize. Under the separable assumption, an interval [ a, b ] which includes ρ ( ǫ ) can be quickly determined bythe following algorithm. 11unction name: Initial .Input: Risk parameter ǫ and confidence parameter δ .Output: [ a, b ] = Initial ( δ, ǫ ).Step 1. Let r ←
1. Apply Probabilistic-Comparison algorithm to compare P (1) with 1 − ǫ . Let theoutcome be d .Step 2. If d = 1 then let d ← d and do the following: • While d = 1 do the following: – Let r ← r . Apply Probabilistic-Comparison algorithm to compare P ( r ) with 1 − ǫ .Let the outcome be d . • Let a ← r and b ← r .Step 3. If d = − d ← d and do the following: • While d = − – Let r ← r . Apply Probabilistic-Comparison algorithm to compare P ( r ) with 1 − ǫ .Let the outcome be d . • Let a ← r and b ← r . Once an initial interval [ a, b ] is obtained, an estimate b R for the probabilistic robustness margin ρ ( ǫ )can be efficiently computed as follows.Function name: Bisection .Input: Risk parameter ǫ , confidence parameter δ , initial interval [ a, b ], relative tolerance γ .Output: [ b R ] = Bisection ( γ, a, b, δ, ǫ ).Step 1. Input γ, a, b, δ, ǫ .Step 2. While b − a > γa do the following: • Let r ← a + b . Apply Probabilistic-Comparison algorithm to compare P ( r ) with 1 − ǫ . Letthe outcome be d . • If d = − b ← r , else let a ← r .12tep 3. Return b R = b (Note: this is actually a soft upper bound).It should be noted that when applying bisection algorithm to refine the initial interval [ a, b ], theexecution of the algorithm may take very long time if P ( r ) ≈ ǫ . However, such chance is almost 0.This problem can be fixed by the following methods. • We can limit the maximum number of simulations to a number M . When conducting simulationat radius r , simulation results can be saved for uncertainty radius a + r where a is the lowerbound of the current interval with middle point r (See the next section for the idea of samplereuse). After the number of simulations for uncertainty radius r exceeds M , the simulation isswitched for uncertainty radius a + r . • In application, one might want to construct the robustness degradation curve in the backwarddirection by starting from an upper bound of ρ ( ǫ ). In that situation, we don’t need to computea very tight interval for ρ ( ǫ ) and hence the chance of P ( r ) ≈ ǫ is even smaller. We shall develop efficient randomized algorithms for constructing robustness degradation curve,which provide more insight for the robustness of the uncertain system than probabilistic robustnessmargin. First we recall the Sample Reuse algorithm [10] for constructing robustness degradationcurve for a given range of uncertainty radius.
Sample Reuse Algorithm
Input: Sample size N , confidence parameter δ ∈ (0 , a, b ], number ofuncertainty radii l .Output: Proportion estimate b P i and the related confidence interval for r i = b − ( b − a )( i − l − , i =1 , , · · · , l . In the following, m i denotes the number of sampling experiments conducted at r i and m i denotes the number of observations guaranteeing P during the m i sampling experi-ments.Step 1 (Initialization). Let M = [ m ij ] l × be a zero matrix.Step 2 (Backward Iteration). For i = 1 to i = l do the following: • Let r ← r i . • While m i < N do the following: – Generate uniform sample q from B ( r ). Evaluate the robustness requirement P for q .13 Let m s ← m s + 1 for any s such that r s ≥ ℓ ( q ). – If robustness requirement P is satisfied for q then let m s ← m s + 1 for any s suchthat r s ≥ ℓ ( q ). • Let b P i ← m i N and construct the confidence interval of confidence level 100(1 − δ )% byTheorem 1.It should be noted that the idea of the Sample Reuse Algorithm is not simply a save of samplegeneration. It is actually a backward iterative mechanism. In the algorithm, the most importantsave of computation is usually the evaluation of the complex robustness requirements P (See [10] fordetails).Now we introduce the global strategy for constructing robustness degradation curve. The idea isto successively apply the Sample Reuse Algorithm for a sequence of intervals of uncertainty radius.Each time the size of interval is reduced by half. The lower bound of the current interval is defined tobe the upper bound of the next consecutive interval. The algorithm is terminated once the robustnessrequirement P is guaranteed for all N samples of an uncertainty set of which the radius is taken asthe lower bound of an interval of uncertainty radius. More precisely, the procedure is presented asfollows. Global Strategy
Input: Sample size N , risk parameter ǫ and confidence parameter δ ∈ (0 , b P i and the related confidence interval.Step 1. Compute an estimate b R for probabilistic robustness margin ρ ( ǫ ).Step 2. Let ST OP ←
0. Let a ← b R and b ← b R .Step 3 (Backward Iteration). While ST OP = 0 do the following: • Apply Sample Reuse Algorithm to construct robustness degradation curve for uncertaintyradius interval [ a, b ]. • If the robustness property P is guaranteed for all N samples of uncertainty set B ( r ) withradius r = a then let ST OP ←
1, otherwise let b ← a and a ← b .For given risk parameter ǫ and confidence parameter δ , the sample size is chosen as N > − ǫ + αǫ )(1 − α ) ln δ α ǫ (6)14ith α ∈ (0 , (cid:26)(cid:12)(cid:12)(cid:12)(cid:12) P X − KN (cid:12)(cid:12)(cid:12)(cid:12) < αǫ (cid:27) > − δ with P X = 1 − ǫ (See also Lemma 1 in Appendix). It should be noted that Massart’s inequality isless conservative than the Chernoff bounds in both multiplicative and additive forms.We would like to remark that the algorithms we propose for estimating the probabilistic robust-ness margin and constructing robustness degradation curve are susceptible of further improvement.For example, the idea of sample reuse is not employed in computing the initial interval and in theprobabilistic bisection algorithm. Moreover, in constructing the robustness degradation curve, theSample Reuse Algorithm is independently applied for each interval of uncertainty radius. Actually,the simulation results can be saved for the successive intervals.We would also like to note that, for a practitioner, computing the probabilistic robustness marginmight be sufficient for understanding the system robustness. However, when more insight about thesystem robustness is expected, the techniques introduced in Section 5 can be employed. Of course,the price is more computational effort. In this section we demonstrate through examples the power of our randomized algorithms in solvinga wide variety of complicated robustness analysis problems which are not tractable in the classicaldeterministic framework.We consider a system which has been studied in [14] by a deterministic approach. The system isas shown in Figure 5.
P(s)(cid:13)C(s)(cid:13) r(cid:13) +(cid:13)_(cid:13) e(cid:13) c(cid:13)
Figure 5: Uncertain SystemThe compensator is C ( s ) = s +2 s +10 and the plant is P ( s ) = . δ ) s ( s +4+0 . δ )( s +6+0 . δ ) with parametricuncertainty ∆ = [ δ , δ , δ ] T . The nominal system is stable. The closed-loop roots of the nominalsystem are: z = − . , z = − . , z = − . . i, z = − . − . i. P peak = 1 . , t r = 0 . , t s = 3 . D -stability of the system. The robustness requirement P is definedas D -stability with the domain of poles defined as: Real part < − .
5, or fall within one of the twodisks centered at z and z with radius 0 .
3. The uncertainty set is defined as the polytope B H ( r ) := ( r ∆ + (1 − r ) P i =1 ∆ i ∈ conv { ∆ , ∆ , ∆ , ∆ } ) where ‘conv’ denotes the convex hull of ∆ i = [ sin( i − π ) , cos( i − π ) , − √ ] T for i = 1 , , = [0 , , T .Obviously, there exists no effective method for computing the deterministic robustness margin inthe literature. However, our randomized algorithms can efficiently construct the robustness degra-dation curve. See Figure 6.In this example, the risk parameter is a priori specified as ǫ = 0 . N denote the sample size which israndom in nature. Let K denote the number of successful trials among N i.i.d. sampling experimentsas defined in Section 2 and Subsection 3 . δ = 0 .
01 and choose tolerance γ = 0 .
05. Staring from r = 1, after N = 7060 simulations we obtain K = 7060, the probabilisticcomparison algorithm determined that P (1) > − ǫ since the lower confidence limit L > − ǫ . Thesimulation is thus switched to uncertainty radius r = 2. After N = 65 times of simulation, it isfound that K = 61. The probabilistic comparison algorithm detected that P (2) < − ǫ because theupper confidence limit U < − ǫ . So, initial interval [1 ,
2] is readily obtained. Now the probabilisticbisection algorithm is invoked. Staring with the middle point of the initial interval (i.e., r = = ),after N = 613 times of simulations, it is found that K = 607, the probabilistic comparison algorithmconcluded that P ( ) < − ǫ since the upper confidence limit U < − ǫ . Thus simulation is movedto r = = . It is found that K = 9330 among N = 9331 times of simulations. Hence, theprobabilistic comparison algorithm determined that P ( ) > − ǫ since the lower confidence limit L > − ǫ . Now the simulation is performed at r = + = . After N = 6653 simulations, it isdiscovered that K = 6636. The probabilistic comparison algorithm judged that P ( ) < − ǫ basedon calculation that the upper confidence limit U < − ǫ . At this point the interval is [ , ] and thebisection is terminated since tolerance condition b − a ≤ γa is satisfied. The evolution of intervalsproduced by the probabilistic bisection algorithm is as follows:[1 , −→ (cid:20) , (cid:21) −→ (cid:20) , (cid:21) −→ (cid:20) , (cid:21) . Now we have obtained an interval [ , ] which includes ρ (0 . l = 100 and the confidence parameter is chosen as δ = 0 . α = 0 . N = 50 , . The interval from which we start constructing robustness degradation curve is [ , ]. It is deter-mined that K = N = 50 ,
632 at uncertainty radius r = . Therefore, the Sample Reuse Algorithmis invoked only once and the overall algorithm is terminated (If K = N for r = , then the nextinterval would be [ , ]). Although P ( r ) is evaluated for l = 100 uncertainty radii with the samesample size N , the total number of simulations is only 153 , << N l = 100 N . To provide anevaluation of accuracy for all estimates of P ( r ), confidence limits are computed by Theorem 1. P r opo r t i on Upper Confidence Limit Lower Confidence Limit Minimum Variance Unbias Estimate
Figure 6: Robustness Degradation CurveWe now apply our algorithms to a robustness problem with time specifications. Specifically,the robustness requirement P is : Stability, and rise time t r < t r = 0 .
25, settling time t s < t s = 3 .
5, overshoot P peak < P peak = 1 .
7. The uncertainty set is B ∞ ( r ) := { ∆ : || ∆ || ∞ ≤ r } .In this case, the risk parameter is a priori specified as ǫ = 0 .
01. It is well known that, for this typeof problem, there exists no effective method for computing the deterministic robustness margin in theliterature. However, our randomized algorithms can efficiently construct the robustness degradationcurve. See Figure 7.We choose γ = 0 .
25 and δ = 0 .
01 for estimating ρ (0 . r = 1,17he initial interval is easily found as [ , ] through the following interval evolution: (cid:20) , (cid:21) −→ (cid:20) , (cid:21) −→ (cid:20) , (cid:21) . The sequence of intervals generated by the probabilistic bisection algorithm is as follows: (cid:20) , (cid:21) −→ (cid:20) , (cid:21) −→ (cid:20) , (cid:21) . So, we obtained an estimate for the probabilistic robustness margin ρ (0 .
01) as . To constructrobustness degradation curve, the number of uncertainty radii is chosen as l = 100 and the confidenceparameter is chosen as δ = 0 .
01. A constant sample size is computed by formula (6) with α = 0 . N = 24 , . The interval from which we start constructing robustness degradation curve is [ , ]. We foundthat this is also the last interval of uncertainty radius because it is determined that K = N atuncertainty radius r = . P r opo r t i on Upper Confidence Limit Lower Confidence Limit Minimum Variance Unbias Estimate
Figure 7: Robustness Degradation Curve
In this paper, we have established efficient techniques which applies to robustness analysis problemswith arbitrary robustness requirements and uncertainty bounding set. The key mechanisms areprobabilistic comparison, probabilistic bisection and backward iteration. Motivated by the crucial18ole of binomial confidence interval in reducing the computational complexity, we have derived anexplicit formula for computing binomial confidence limits. This formula overcomes the computationalissue and inaccuracy of standard methods.
A Proof of Theorem 1
To show Theorem 1, we need some preliminary results. The following Lemma 1 is due to Massart[18].
Lemma 1 Pr (cid:8) KN ≥ P X + ǫ (cid:9) ≤ exp (cid:16) − Nǫ P X + ǫ ) (1 − P X − ǫ ) (cid:17) for all ǫ ∈ (0 , − P X ) . Of course, the above upper bound holds trivially for ǫ ≥ − P X . Thus, Lemma 1 is actually truefor any ǫ > Lemma 2 Pr (cid:8) KN ≤ P X − ǫ (cid:9) ≤ exp (cid:16) − Nǫ P X − ǫ ) (1 − P X + ǫ ) (cid:17) for all ǫ > . Proof . Define Y = 1 − X . Then P Y = 1 − P X . At the same time when we are conducting N i.i.d. experiments for X , we are also conducting N i.i.d. experiments for Y . Let the numberof successful trials of the experiments for Y be denoted as K Y . Obviously, K Y = N − K . Ap-plying Lemma 1 to Y , we have Pr n K Y N ≥ P Y + ǫ o ≤ exp (cid:16) − Nǫ P Y + ǫ ) (1 − P Y − ǫ ) (cid:17) . It follows thatPr (cid:8) N − KN ≥ − P X + ǫ (cid:9) ≤ exp (cid:16) − Nǫ − P X + ǫ ) [1 − (1 − P X ) − ǫ ] (cid:17) . The proof is thus completed by observ-ing that Pr (cid:8) N − KN ≥ − P X + ǫ (cid:9) = Pr (cid:8) KN ≤ P X − ǫ (cid:9) . ✷ The following lemma can be found in [13].
Lemma 3 P kj =0 (cid:0) Nj (cid:1) x j (1 − x ) N − j decreases monotonically with respect to x ∈ (0 , for k = 0 , , · · · , N . Lemma 4 P kj =0 (cid:0) Nj (cid:1) x j (1 − x ) N − j ≤ exp (cid:18) − N ( x − kN ) x + k N ) (1 − x − k N ) (cid:19) ∀ x ∈ ( kN , for k = 0 , , · · · , N . Proof.
Consider binomial random variable X with parameter P X > kN . Let K be the number ofsuccessful trials during N i.i.d. sampling experiments. Then P kj =0 (cid:0) Nj (cid:1) P jX (1 − P X ) N − j = Pr { K ≤ k } .Note that Pr { K ≤ k } = Pr (cid:8) KN ≤ P X − (cid:0) P X − kN (cid:1)(cid:9) . Applying Lemma 2 with ǫ = P X − kN >
0, wehave k X j =0 (cid:18) Nj (cid:19) P jX (1 − P X ) N − j ≤ exp − N ( P X − kN ) P X − P X − kN ) (1 − P X + P X − kN ) = exp − N ( P X − kN ) P X + k N ) (1 − P X − k N ) ! . X with P X > kN , thus the proof ofthe lemma is completed. (cid:3) Lemma 5 P k − j =0 (cid:0) Nj (cid:1) x j (1 − x ) N − j ≥ − exp (cid:18) − N ( x − kN ) x + k N ) (1 − x − k N ) (cid:19) ∀ x ∈ (0 , kN ) for k = 1 , · · · , N . Proof.
Consider binomial random variable X with parameter P X < kN . Let K be the number ofsuccessful trials during N i.i.d. sampling experiments. Then k − X j =0 (cid:18) Nj (cid:19) P jX (1 − P X ) N − j = Pr { K < k } = Pr (cid:26) KN < P X + ( kN − P X ) (cid:27) . Applying Lemma 1 with ǫ = kN − P X >
0, we have that k − X j =0 (cid:18) Nj (cid:19) P jX (1 − P X ) N − j ≥ − exp − N ( kN − P X ) P X + kN − P X ) (1 − P X − kN − P X ) = 1 − exp − N ( P X − kN ) P X + k N ) (1 − P X − k N ) ! . Since the argument holds for arbitrary binomial random variable X with P X < kN , thus the proof ofthe lemma is completed. (cid:3) Lemma 6
Let ≤ k ≤ N . Then L N,k,δ < U
N,k,δ . Proof.
Obviously, the lemma is true for k = 0 , N . We consider the case that 1 ≤ k ≤ N − S ( N, k, x ) := P kj =0 (cid:0) Nj (cid:1) x j (1 − x ) N − j for x ∈ (0 , S ( N, k, p ) = S ( N, k − , p ) + (cid:0) Nk (cid:1) p k (1 − p ) N − k = δ . Thus S ( N, k − , p ) − S ( N, k − , p ) = 1 − δ − (cid:20) δ − (cid:18) Nk (cid:19) p k (1 − p ) N − k (cid:21) . Notice that δ ∈ (0 ,
1) and that p ∈ (0 , S ( N, k − , p ) − S ( N, k − , p ) = 1 − δ + (cid:18) Nk (cid:19) p k (1 − p ) N − k > . By Lemma 3, S ( N, k − , x ) decreases monotonically with respect to x , we have that p < p . (cid:3) We are now in the position to prove Theorem 1. For national simplicity, let p = U N,k,δ , q = U ( N, k, δ ) . It can be easily verified that p ≤ q for k = 0 , N . We need to show that p ≤ q for 0 < k < N .Straightforward computation shows that q is the only root of equationexp − N ( x − kN ) x + k N ) (1 − x − k N ) ! = δ x ∈ ( kN , ∞ ). There are two cases: q ≥ q <
1. If q ≥ p ≤ q is triviallytrue. We only need to consider the case that kN < q <
1. In this case, it follows from Lemma 4 that k X j =0 (cid:18) Nj (cid:19) q j (1 − q ) N − j ≤ exp − N ( q − kN ) q + k N ) (1 − q − k N ) ! = δ . Recall that k X j =0 (cid:18) Nj (cid:19) p j (1 − p ) N − j = δ , we have k X j =0 (cid:18) Nj (cid:19) p j (1 − p ) N − j ≥ k X j =0 (cid:18) Nj (cid:19) q j (1 − q ) N − j . Therefore, by Lemma 3, we have that p ≤ q for 0 < k < N . Thus, we have shown that U N,k,δ ≤ U
N,k,δ for all k .Similarly, by Lemma 5 and Lemma 3 we can show that L N,k,δ ≥ L ( N, k, δ ). By Lemma 6, wehave L ( N, k, δ ) < L N,k,δ < U
N,k,δ < U ( N, k, δ ). Finally, the proof of Theorem 1 is completed byinvoking the probabilistic implication of the Clopper-Pearson confidence interval.
References [1] E. W. BAI, R. TEMPO, AND M. FU, “Worst-case properties of the uniform distribution andrandomized algorithms for robustness analysis,”
Mathematics of Control, Signals and Systems ,vol. 11(1998), pp.183-196.[2] B. R. BARMISH, C. M. LAGOA, AND R. TEMPO, “Radially truncated uniform distribu-tions for probabilistic robustness of control systems,”
Proc. of American Control Conference ,Albuquerque, New Mexico, June, 1997, pp. 853-857.[3] B. R. BARMISH, “A probabilistic result for a multilinearly parameterized H ∞ norm,” Proc. ofAmerican Control Conference , Chicago, Illinois, June, 2000, pp. 3309-3310.[4] B. R. BARMISH AND B. T. POLYAK, “A new approach to open robustness problems basedon probabilistic predication formulae,” IFAC’1996, San Francisco, Vol. H, pp. 1-6.[5] B. R. BARMISH AND C. M. LAGOA, “The uniform distribution: a rigorous justification forits use in robustness analysis,”
Mathematics of Control, Signals and Systems , vol. 10 (1997),pp.203-222.[6] G. CALAFILORE, F. DABBENE, AND R. TEMPO, “Randomized algorithms for probabilisticrobustness with real and complex structured uncertainty,”
IEEE Transaction on AutomaticControl , Vol. 45, No.12, December, 2000, pp. 2218-2235.217] G. CALAFILORE, F. DABBENE, AND R. TEMPO, “Radial and uniform distributions invector and matrix spaces for probabilistic robustness,”
Topics in Control and Its Applications ,Miller D. and Qiu L., Eds, Spring-Verlag, 1999, pp. 17-31.[8] X. CHEN AND K. ZHOU, “Order statistics and probabilistic robust control,”
Systems andControl Letters , vol. 35(1998), pp. 175-182.[9] X. CHEN AND K. ZHOU, “Constrained robustness analysis and synthesis by randomized al-gorithms,”
IEEE Transaction on Automatic Control , Vol. 45, No.6, June, 2000, pp. 1180-1186.[10] X. CHEN, K. ZHOU AND J. L. ARAVENA, “Fast construction of robustness degradationfunction”,
Proceedings of Conference on Decision and Control , Las Vagas, Nevada, December2002.[11] X. CHEN, K. ZHOU AND J. L. ARAVENA, “Explicit formula for constructing binomial confi-dence interval with guaranteed coverage probability,”
Communications in Statistics — Theoryand Methods , vol. 37, pp. 1173–1180, 2008.[12] C. J. CLOPPER and E. S. PEARSON, “The use of confidence or fiducial limits illustrated inthe case of the binomial,”
Biometrika , vol. 26, pp.404-413, 1934.[13] C. W. CLUNIES-ROSS, “Interval estimation for the parameter of a binomial distribution”,
Biometrika , vol. 45, pp. 275-279, 1958.[14] R. R DE GASTON AND M. G. SAFONOV, “Exact calculation of the multiloop stability mar-gin,”
IEEE Trans. Autom. Control , vol. 33(1988), pp. 156-171.[15] A. HALD,
Statistical Theory with Engineering Applications , pp. 697-700, John Wiley and Sons,1952.[16] P. P. KHARGONEKAR AND A. TIKKU, “Randomized algorithms for robust control analysisand synthesis have polynomial complexity,”
Proceedings of the 35th Conference on Decision andControl , Kobe, Japan, December 1996, pp. 3470-3475.[17] C. M. LAGOA, “Probabilistic enhancement of classic robustness margins: A class of nonesymmetric distributions ,”
Proc. of American Control Conference , Chicago, Illinois, June, 2000,pp. 3802-3806.[18] P. MASSART, “The tight constant in the Dvoretzky-Kiefer-Wolfowitz inequality”,
The Annalsof Probability , pp. 1269-1283, Vol. 18, No. 3, 1990.2219] C. MARRISON AND R. F. STENGEL, “Robust control system design using random searchand genetic algorithms,”
IEEE Transaction on Automatic Control , Vol. 42, No.6, June, 1997,pp. 835-839.[20] S. R. ROSS AND B. R. BARMISH, “Distributionally Robust Gain Analysis for systems contain-ing complexity,”
Proceedings of the 40th Conference on Decision and Control , Orlando, Florida,December 2001, pp. 5020-5025.[21] L. R. RAY AND R. F. STENGEL, “A monte carlo approach to the analysis of control systemsrobustness,”
Automatica , vol. 3(1993), pp. 229-236.[22] R. F. STENGEL AND L. R. RAY, “Stochastic robustness of linear time-invariant systems,”
IEEE Transaction on Automatic Control , AC-36(1991), pp. 82-87.[23] B. T. POLYAK AND P. S. SHCHERABAKOV, “Random spherical uncertainty in estimationand robustness,”
IEEE Trans. Autom. Control , vol. 45(2000), pp. 2145-2150.[24] R. TEMPO, E. W. BAI, AND F. DABBENE, “Probabilistic robustness analysis: explicit boundsfor the minimum number of samples,”
Systems and Control Letters , vol. 30 (1997), pp. 237-242.[25] R. TEMPO AND F. DABBENE, “Probabilistic robustness analysis and design of uncertainsystems,”
Progress in Systems and Control Theory , (ed. G. Picci) Birkhauser, March 1999, pp.263-282.[26] M. VIDYASAGAR AND V. D. BLONDEL, “Probabilistic solutions to NP-hard matrix prob-lems”,
Automatica , vol. 37(2001), pp. 1597-1405.[27] M. VIDYASAGAR, “Randomized algorithms for robust controller synthesis using statisticallearning theory”,