Parallel Branch and Bound Algorithm for Computing Maximal Structured Singular Value
aa r X i v : . [ m a t h . O C ] M a y Parallel Branch and Bound Algorithm for ComputingMaximal Structured Singular Value ∗ Xinjia Chen and Kemin ZhouDepartment of Electrical and Computer EngineeringLouisiana State UniversityBaton Rouge, LA [email protected] [email protected] 1, 2002, Revised on March 18, 2003
Abstract
In this paper, we have developed a parallel branch and bound algorithm which computes themaximal structured singular value µ without tightly bounding µ for each frequency and thussignificantly reduce the computational complexity. Keywords:
Robust control, structured singular value, branch and bound.
It is well known that the analysis of robust stability and performance with structured uncertaintyboils down to the problem of computing the supremum of the structured singular value over allfrequency [3, 11]. That is, µ max := sup ω ∈ R µ ∆ ( M ( jω )) where M ( s ) is the transfer function of thegeneralized system and ∆ is a set of block structured uncertainties. Related to this problem arethe method proposed by Lawrence, Tits and Dooren [7, 8] and the approach established by Ferreresand Biannic [5, 6]. These interesting techniques can be applied to compute a µ upper bound overa frequency interval without gridding of frequency. For the precise computation of the maximalstructured singular value µ max (i.e., a tight lower bound is also expected in addition to an upperbound), the conventional method is to grid a range of frequency and compute the maximal µ amongall the frequencies [1]. Since the exact computation is in general impossible, µ is obtained for eachfrequency by tightly bounding. Sophisticated upper bounds and lower bounds have been derived for ∗ This research was supported in part by grants from AFOSR (F49620-94-1-0415), ARO (DAAH04-96-1-0193), andLEQSF (DOD/LEQSF(1996-99)-04). µ max lack of efficiency because of the tedious frequency sweeping. In this paper, we investigate a smartfrequency sweeping strategy. More specifically, we apply branch and bound scheme to compute µ for N > ˆ µ − ǫ where ˆ µ is the maximum record of the lower boundsof all branches ever generated and ǫ > µ is returned as the maximalstructured singular value µ max . Since ˆ µ is the maximum record of the lower bounds obtained in allbranches generated (no matter belong to the same frequency or not), it will increase much fasterthan its counterpart in the conventional frequency sweeping algorithms. Note that the raise of ˆ µ results in a significant number of branches to be pruned. Thus ˆ µ convergences quickly to the maximalstructured singular value µ max .The paper is organized as follows. Section 2 discusses existing techniques for computing themaximal structured singular value. Section 3 presents our Parallel Branch and Bound Algorithm.An illustrative example is provided in Section 4 and Section 5 is the conclusion. Consider an M − ∆ set up as follows. ∆ M Figure 1: Uncertain System.Let ∆ be a set of block structured uncertainties. We consider the computation of µ max := sup ω ∈ R µ ∆ ( M ( jω )) . For notation simplicity, let M ( ω ) = M ( jω ). Then µ max = sup ω ∈ R µ ∆ ( M ( ω )).In practice, it is impossible to search µ max over all frequencies. However, we can estimate µ max as follows.Choose a range of frequency [ a, b ] ∈ R and grid it as ω j = a + ( b − a )( j − N K − , j = 1 , · · · , N K (1)2here N ≥ K ≥ µ max can be defined as ˜ µ max := max j =1 , ··· ,NK µ ∆ ( M ( ω j )) . Define the (maximum positive real eigenvalue) function ¯ λ R : C n × n → R as¯ λ R ( M ) := max { λ : λ is a positive real eigenvalue of M } with ¯ λ R ( M ) = 0 if M has no positive real eigenvalues. Let B∆ := { ∆ ∈ ∆ : ¯ σ (∆) ≤ } . Then µ ∆ ( M ) = max ∆ ∈ B∆ ¯ λ R ( M ∆) . Let Q ⊂ B∆ . Define µ on a box [9] µ ( M, Q ) := max ∆ ∈ Q ¯ λ R ( M ∆) . There exists techniques in [9] for computing an upper bound
U B ( M, Q ) and a lower bound LB ( M, Q )for µ ( M, Q ). Thus a branch and bound scheme can be applied to compute µ ∆ ( M ) with parameterspace B∆ .The conventional methods work essentially as follows. For j = 1 , · · · , N K , apply the followingAlgorithm 1 to compute an upper bound U B j and a lower bound LB j for µ ∆ ( M ( ω j )) such that U B j − LB j ≤ ε . Then ˜ µ max satisfiesmax j =1 , ··· ,NK LB j ≤ ˜ µ max ≤ max j =1 , ··· ,NK U B j where max j =1 , ··· ,NK U B j − max j =1 , ··· ,NK LB j ≤ ε. Algorithm — Branch and Bound ([9]) Initialize
Let U j = { Q k } = B∆ . Let
U B j = max k U B ( M ( ω j ) , Q k ) ,LB j = max k LB ( M ( ω j ) , Q k ) . (2) while U B j − LB j > ε • Choose Q to be any element of U j with U B ( M ( ω j ) , Q ) = U B j . • Partition Q into Q a and Q b by bisecting along one of its longest edges. • Add Q a and Q b into U j . Remove Q from U j .3 Remove from U j any Q with U B ( M ( ω j ) , Q ) < LB j . (3) endwhile The most important mechanism of Algorithm 1 is “pruning” [9]. That is, any element of U j forwhich (3) is satisfied will never again be partitioned and need not be considered further. We callinequality (3) as the “pruning condition”.We can see that existing techniques for computing ˜ µ max employ branch and bound techniques foreach frequency independently. In particular, the pruning process for one frequency is independentof another. µ ∆ ( M ) is bounded tightly for each frequency. Note that we usually need to evaluate µ ∆ ( M ) for many frequencies in order to obtain a reasonably good estimate of the maximal structuredsingular value µ max . Thus the overall computation is still a heavy burden, even the computation of µ ∆ ( M ) for each frequency is very efficient.Thus for the sake of efficiency, there is a strong motivation to conceive a smart frequency sweepingstrategy. More specifically, we would raise the following question, Is it possible to obtain the maximal structured singular value µ max without tightly bounding µ ∆ ( M ( ω j )) for each frequency ω j ? The following section is devoted to answering this question.
It is fair to compare the performance of different algorithms on the same set of frequencies. Therefore,we consider again frequencies ω j , j = 1 , · · · , N K defined by (1) and relabel them as ω ij := a + ( b − a )[ K ( i −
1) + ( j − N K − , i = 1 , · · · , N, j = 1 , · · · , K. Now we are in a good position to present our Parallel Branch and Bound Algorithm as follows.
Algorithm — Parallel Branch and Bound Algorithm • Step 1: Initialize. Set j = 1. Set ˆ µ = 0. Set tolerance ǫ >
0. Set maximal iteration number IT . • Step 2: Update ˆ µ and record the number of iterations r ( j ) for frequency ω ij by the followingsteps. – Step 2–1: Let U ij = { Q k } = B∆ , i = 1 , · · · , N . Set r = 1. – Step 2–2: If r = IT + 1 or U ij is empty for any i ∈ { , · · · , N } then record r ( j ) = r andgo to Step 3, else do the following for all i such that U ij is not empty.4 Choose Q to be any element of U ij with U B ( M ( ω ij ) , Q ) = max Q k ∈U ij U B ( M ( ω ij ) , Q k ) . ∗ Partition Q into Q a and Q b by bisecting along one of its longest edges. ∗ Add Q a and Q b into U ij . Remove Q from U ij . ∗ Update ˆ µ = max { ˆ µ, LB ( M ( ω ij ) , Q a ) , LB ( M ( ω ij ) , Q b ) } . (4) ∗ Remove from U ij any Q with U B ( M ( ω ij ) , Q ) < ˆ µ − ǫ . (5) – Step 2–3: Set r = r + 1 and go to Step 2–2. • Step 3: If j = K then STOP, else set j = j + 1 and go to Step 2.In Algorithm 2, N branches of frequency sweeping are performed in parallel with starting fre-quencies ω i , i = 1 , · · · , N and step size b − aNK − . Also, a branch and bound scheme is applied tocompute µ for N frequencies in parallel. Any branch with upper bound smaller than ˆ µ − ǫ will bepruned, where ˆ µ is the maximum record of the lower bounds of all branches ever generated. Thefinal ˆ µ is returned as the maximal structured singular value µ max . Algorithm 2 is visulized in thefollowing Figure 2. Remark 1
Note that Algorithm provides a substantial improvement on efficiency than conven-tional methods in computing the maximal structured singular value. This can be explained by thesignificant relaxation in the “pruning condition” of Algorithm . To see the difference of the two“pruning conditions”, we can compare the right hand sides of inequalities (5) and (3). By (4) and(2), we can see that ˆ µ − ǫ can be much larger than LB j . This is because ˆ µ is the maximum recordof the lower bounds obtained in all branches of all frequencies evaluated and being evaluated, while LB j is only the maximum record of the lower bounds obtained in branches of the frequency beingevaluated. Moreover, ˆ µ is enlarged to ˆ µ − ǫ in the “pruning condition” (5) and hence the “pruning”process is further facilitated. The significant relaxation of the ‘pruning condition” leads to a substan-tial decrease of the number of total subdomains needed to be evaluated. Therefore, our algorithm ismuch more efficient than those previously available to control engineers. Remark 2
It is important to note that Algorithm involves only one CPU processor. It is funda-mentally different from the parallel algorithms which involves more than one CPU processors. B B ∆∆∆
Frequencyr=1r=2r=3 ω ω ω Figure 2: A Picture of Parallel Branch and Bound Algorithm. N = 3 , K = 4. Remark 3
A substantial amount of computation can be saved by the following mechanisms. First,further computation of the lower bound on a domain is not needed once it is determined that thelower bound is smaller than the existing global lower bound. This can be seen from equation (4).Second, the computation of the upper bound should be terminated once condition (5) is satisfied. Theidea of these two mechanisms is to avoid as much as possible the tightly computation of the lowerbound and the upper bound.
In addition to the novel frequency sweeping strategy, another character of Algorithm 2 is thatthere is no tolerance criteria directly forced on the final result, however, the final result falls intotolerance automatically.
Theorem 1
Suppose that the maximal iteration number
IT < ∞ and that Algorithm stops with r ( j ) ≤ IT, j = 1 , · · · , K . Then the final ˆ µ satisfies ≤ ˜ µ max − ˆ µ ˜ µ max < ǫ. Proof . Since ˆ µ is the maximal record of the lower bounds, we have ˜ µ max − ˆ µ ≥
0. We only need toshow that ˜ µ max − ˆ µ ˜ µ max < ǫ . By the assumption that Algorithm 2 stops with r ( j ) ≤ IT, j = 1 , · · · , K , weknow that all subdomains ever generated are finally removed because the “pruning condition” (5) is6atisfied. Note that there exists a subdomains Q ij for frequency ω ij such that µ ( M ( ω ij ) , Q ij ) = ˜ µ max .Let ˆ µ = ¯ µ when Q ij is removed. Then ˜ µ max ≤ U B ( M ( ω ij ) , Q ij ) < ¯ µ − ǫ . Note that ˆ µ is nondecreasingthus the final ˆ µ ≥ ¯ µ . It follows that˜ µ max < ˆ µ − ǫ = ⇒ ˜ µ max − ˆ µ ˜ µ max < ǫ. The proof is thus completed. ✷ Note that one important concern of an algorithm is convergence. It is usually desirable that,given any tolerance ǫ >
0, an algorithm stops and returns the result within tolerance in a finitenumber of iterations. Obviously, the convergence requirement imposes condition of the quality ofbounds.
Definition 1
The upper bound
U B ( M, . ) and lower bound LB ( M, . ) are said to be continuous if lim d ( Q ) → U B ( M, Q ) − LB ( M, Q ) = 0 where d ( Q ) := max q, q ′ ∈ Q || q − q ′ || with Q ⊆ B∆ . Theorem 2
Suppose that all the upper bounds and lower bounds are continuous and that at least onenonzero lower bound appears after a finite number of iterations. Let the maximal iteration number IT = ∞ . Then, for arbitrary tolerance ǫ > , Algorithm stops with a finite number of domainpartitions for each j , i.e., r ( j ) < ∞ , j = 1 , · · · , K . Moreover, the final ˆ µ satisfies ≤ ˜ µ max − ˆ µ ˜ µ max < ǫ. Proof . Suppose that Algorithm 2 does not stop with a finite number of domain partitions for each j . Then ∃ ω ij and an infinite sequence of nested subdomains { Q ijr } associated with frequency ω ij such that Q ij ⊃ Q ij ⊃ · · · ⊃ Q ijr ⊃ · · · . Note that by the assumption ∃ r < ∞ , µ > µ ≥ µ , ∀ r > r . By the continuity, ∃ r such that U B ( M ( ω ij ) , Q ijr ) − LB ( M ( ω ij ) , Q ijr ) < ǫ − ǫ µ , ∀ r > r . Let r = max { r , r } + 1. Then U B ( M ( ω ij ) , Q ijr ) − LB ( M ( ω ij ) , Q ijr ) < ǫ − ǫ µ . Thus
U B ( M ( ω ij ) , Q ijr +1 ) − ˆ µ < ǫ − ǫ ˆ µ = ⇒ U B ( M ( ω ij ) , Q ijr +1 ) < ˆ µ − ǫ Q ijr +1 is removed. This is a contradiction. Therefore Algorithm 2 stops with afinite number of domain partitions for each j and hence by the same argument of Theorem 10 ≤ ˜ µ max − ˆ µ ˜ µ max < ǫ. The proof is thus completed. ✷ Consider an M − ∆ set up as shown in Figure 1 where ∆ = diag( δ , δ ) ∈ R × and M ( s ) = C ( sI − A ) − B with A = − − − − . − . . − − − −
10 0 . − . , B = , C = " − . − . . To compute the supremum of µ , we uniformly grid frequency interval [0 . , .
01] and obtain 1 , ω j = 0 .
01 + (15 . − . j − − , j = 1 , · · · , . In Algorithm 2, we choose the relative error ǫ = 0 .
01 and N = 30 , K = 50. The 1 ,
500 frequenciesare regrouped as ω ij = 0 .
01 + (15 . − . i −
1) + ( j − − , i = 1 , · · · , j = 1 , · · · , . The execution of Algorithm 2 is terminated at r = 1 with ˆ µ = 0 . ω ij =9 . i = 19 , j = 16. It is observed that, for any frequency, no partition is performed forthe original domain B∆ = [ − , × [ − , ˆ µ − ǫ is greater than the upper boundsof singular values for B∆ at other frequencies. The bounds of singular values for B∆ are shown byFigures 3-4. It can be seen from these figures that the upper bounds and lower bounds of singularvalues are far apart for most of the frequencies. To compute the maximal singular value using theconventional branch and bound method, substantial computational effort will be wasted on reducingthe gap between the upper bounds and lower bounds of singular values for most of the frequencies.This example demonstrates that branch and bound should not be applied extensively for any fixedfrequency. Quiet contrary, it should be employed in parallel and in a cooperative manner. This spirithas been reflected in Algorithm 2. 8 B ound s o f S i ngu l a r V a l ue m u Upper boundLower bound
Figure 3: Bounds of Singular Values
Efficient computation of the maximal structured singular value is of fundamental importance inrobustness analysis and robust synthesis with structured uncertainty. Motivated by this, we havedeveloped a parallel branch and bound algorithm for computing the maximal structured singularvalue, which significantly reduce the computational complexity.
References [1] Balas, G., Doyle, J. C., Glover, K., Packard, A. and Smith, R., µ -Analysis and Synthesis Toolbox ,MUSYN Inc. and The MathWorks, Inc. 1995.[2] Balas, G. and Packard A., “The structured singular value ( µ ) framework,” The Control HandBook , pp. 671-687, CRC Press, Inc, 1996.[3] Doyle, J. C., “Analysis of feedback system with structured uncertainty,”
IEE Proc. , pt. D, vol.129, no. 6, pp. 242-250, 1982.[4] Fan, M. K. H. and Tits, A. L., “Characterization and efficient computation of the structuredsingular value,”
IEEE Trans. Autom. Control , vol. 31, pp. 734-743, 1986.9 B ound s o f S i ngu l a r V a l ue m u Upper boundLower bound
Figure 4: Bounds of Singular Values[5] Ferrers and Biannic, “A µ -analysis technique without frequency gridding,” Proceeding of Amer-ican Control Conference , pp. 2294-2298, 1998.[6] Ferrers and Biannic, “Reliable computation of the robustness margin for a flexible aircraft,”
Control Engineering Practice , pp. 1267-1278, December 2001.[7] A. Helmersson, “A finite frequency method for µ -analysis,” Proceeding of European ControlConference , pp. 171-176, 1995.[8] Lawrence, C. T., Tits, A. L. and Dooren, P. V., “A fast algorithm for the computation of anupper bound on the µ -norm,” Automatica , vol. 36, pp. 449-456, 2000.[9] Newlin M. P. and Young P. M., “Mixed µ problems and branch and bound techniques,” Proc.IEEE Conference on Decision and Control , pp. 3175-3180, Tucson, Arizona, 1992.[10] Packard, A. and Doyle, J. C., “The complex structured singular value,”
Automatica , 29(1), pp.71-109, 1993.[11] Zhou, K., Doyle, J. C. and Glover, K.,