A detailed comparative analysis of μ lower bound algorithms

There is something striking with $\mu$ analysis. On the one hand, it appears to be a very mature technique. Several different methods have been developed in the last 30 years in order to tackle the problem of computing the structured singular value $\mu$, and very relevant results are reported in a number of papers. On the other hand, only few well-documented softwares are widely available, and many practical algorithms have only been evaluated by their respective authors. As a result, control engineers and researchers often have to rely on the algorithms implemented in the Robust Control toolbox for Matlab. This observation raizes an important issue: are these tools state-of-the-art or can they be outperformed? 

Another difficulty is that the exact computation of $\mu$ is NP hard in the general case, so both lower and upper bounds have to be computed instead. Most of the computationally tractable methods to compute $\mu$ upper bounds rely on the so-called $D$ and $G$ scaling matrices. By contrast, there is no such consensus regarding $\mu$ lower bounds. Indeed, a wide number of very different approaches have been developed, as confirmed by the substantial bibliography detailed below. But these bounds are of little help if the gap with respect to the exact value of $\mu$ is large or if the computational cost is excessive. A good method should thus offer a satisfactory compromise between conservatism and CPU time. 

In this context, all practical algorithms to compute lower bounds on the structured singular value $\mu$ (i.e. upper bounds on the robust stability margin $k_r$) are evaluated here on a wide set of real-world applications. This is probably the best way to highlight the most promising tools from an engineering perspective, but also to propose appropriate improvements and combinations to get the most out of them. The ultimate goal is of course to bridge the gap between the real value of $\mu$ and its bounds. Another ambition of this work is to provide control engineers with a well-documented software, which can be used without mastering all the underlying theory. In this perspective, new and efficient implementations of the most relevant algorithms and strategies are available in the Skew Mu Analysis based Robustness Tools (SMART) library of the SMAC toolbox, which is currently developed by ONERA (The French Aerospace Lab).

All numerical data related to the benchmarks can be downloaded here. Each benchmark is described either by an uncertain state-space model usys, or by standard matrices sys and blk which can be directly passed to low-level routines such as mussv, mulb_nreal or mulb_mixed. Both descriptions are of course equivalent and can be used with Matlab R2012b or higher.

This page is organized as follows. All approaches to compute lower bounds on the structured singular value are first surveyed. Detailed numerical results are then presented and analyzed. An extension to skew-μ lower bound computation is finally proposed.

Exhaustive survey of $\mu$ lower bound algorithms

All methods that can reasonably be applied to real-world benchmarks are mentioned. Only a few exponential-time algorithms are omitted, as well as techniques with a very limited range of application. In all cases, both $\mu$ lower bounds and associated destabilizing values of the uncertainties are obtained. For most methods, frequency is fixed during optimization. In this case, $M(j\omega)$ and $\Delta(j\omega)$ are constant matrices, which are simply denoted by $M$ and $\Delta$.

1. Exponential-time methods

The first methods to compute $\mu$ lower bounds were proposed in the 1980's and most of them are exponential-time. One of the best known was introduced in [3] for non-repeated ($n_i=1$) real uncertainties and later generalized in [4,5,6]. It uses the mapping theorem of [7]: the image of $k\mathcal{B}_{{\bf\Delta}}$ by the operator $\Delta\rightarrow\textrm{det}(I-\Delta M)$ is included into the convex hull of the images of the $2^N$ vertices of $k\mathcal{B}_{{\bf\Delta}}$. A $\mu$ lower bound is obtained by considering the position of these images with respect to the origin in the complex plane. Like the mapping theorem, the Kharitonov's theorem [8] and the edge theorem [9] can also be interpreted in a frequency domain approach and used to derive $\mu$ lower bound algorithms. These are not detailed here, but a few interesting references can be found in [5,6].

Then, the algebraic method proposed in [10] for non-repeated real uncertainties searches for a destabilizing uncertainty $\Delta \in k \mathcal{B}_\mathbf{\Delta}$ such that all the $\delta_i$ except 2 attain the maximal magnitude $k$, which can be achieved by simple matrix algebra operations. A one-dimensional halving search on $k$ is performed to compute a $\mu$ lower bound. There are $N(N-1)/2$ ways to choose $2$ uncertainties among $N$ and $2^{N-2}$ ways to set each of the $N-2$ others to $-k$ or $k$, so $N(N-1)2^{N-3}$ searches are performed for each value of $k$. The method presented in [11] is quite similar. It is based on the computation of the smallest common real root of a system of linear polynomials using the Sylvester's theorem [12]. The main difference with [10] is that both $k$ and the two $\delta_i$ which do not attain the maximal magnitude are computed in a single step. The optimal value of $k$ is thus determined directly and no halving search is required.

More recently, a method was presented in [13] for (possibly repeated) real uncertainties. It consists of introducing $b=\prod_{i=1}^{N}(n_i+1)-1$ fictitious uncertainties $\overline{\delta}_1,\dots,\overline{\delta}_b$ and rewriting the numerator of $\textrm{det}(I-M(s)\Delta)$ as $f(s)=\hat{p}(s)+\left[\,\overline{\delta}_1\dots\overline{\delta}_b\right] \left[p_1(s)\dots p_b(s)\right]^T$, where $\hat{p}(s),p_1(s),\dots,p_b(s)$ are fixed real polynomials. The stability of $f(s)$ is then evaluated using the stability feeler introduced in [14], which leads to $\mu$ upper and lower bounds.

The algorithms mentioned in this section have different theoretical foundations. But they are all exponential-time, which means that they cannot be applied to most real-world benchmarks with a reasonable computational time. Thus, only polynomial-time algorithms are considered in the sequel.

2. Power algorithm

This technique can be used for all kinds of uncertainties. It is based on the following characterization [2]:

where $\rho_R(QM)$ is the largest magnitude of the real eigenvalues of $QM$, while $\mathcal{Q}$ is the set of all $\Delta=\textrm{diag}(\Delta_1,\dots,\Delta_N)\in{\bf\Delta}$ such that $\delta_i\in [-1,1]$ if $\Delta_i=\delta_iI_{n_i}$ is real and $\Delta_i^*\Delta_i=I_{n_i}$ if $\Delta_i$ is complex. Equation (1) defines a non-concave optimization problem, which is solved in [2] using a fixed point iteration usually referred to as the power algorithm. No theoretical convergence guarantee exists. But if the algorithm converges, a local maximum, i.e. a $\mu$ lower bound, is obtained. In practice, convergence is usually ensured for purely complex and mixed real/complex uncertainties. But in the purely real case, limit cycles often appear and no $\mu$ lower bound can be obtained. Nevertheless, a technique is proposed in [15] to get a nontrivial $\mu$ lower bound when the power algorithm does not converge but returns nonzero values. Further improvements are also reported in [16,17,18].

A different way to solve problem (1) is proposed in [19] in the particular case of purely complex uncertainties ($\rho_R$ is then replaced with the spectral radius $\rho$). The steepest ascent and the conjugate gradient algorithms are used instead of a fixed-point iteration. Unlike the power algorithm, which only provides a $\mu$ lower bound in case of convergence, a nontrivial bound is obtained at each iteration at the price of a higher computational cost.

3. Gain-based algorithm

This method was introduced in [20] to address the convergence issues of the power algorithm when purely real uncertainties are considered. The initial $\mu$ problem is replaced with a series of worst-case $H_\infty$ performance problems, which are efficiently solved using the approach of [21]. In the mixed real/complex case, the real blocks of the destabilizing uncertainties are obtained using the aforementioned technique, while the complex blocks are computed using the power algorithm.

4. Poles migration techniques

These techniques are mainly dedicated to purely real problems, although they can in some cases be extended to handle mixed and complex uncertainties. The idea is to move an eigenvalue of the interconnection between $M(s)$ and $\Delta$ towards the imaginary axis, which requires to work in the time domain. The state matrix of the interconnection is $A_0=A+B\Delta(I-D\Delta)^{-1}C$, where $(A,B,C,D)$ is a state-space representation of $M(s)$. An equivalent definition of the structured singular value can thus be derived:

\begin{equation} \mu_{\bf\Delta}(M(j\omega)) = \Big[\min_{\Delta\in{\bf\Delta}}\left\{\overline{\sigma}(\Delta),\,\lambda_{max}(A_0)=0\right\}\Big]^{-1} \end{equation}

where $\lambda_{max}(A_0)=\displaystyle\max_i\,\Re\left(\lambda_i(A_0)\right)$ and $\lambda_i(A_0)$ denotes the $ith$ eigenvalue of $A_0$.

A first approach is to use a first-order characterization of the variation $d\lambda_i$ of $\lambda_i(A_0)$ caused by a small variation $d\Delta$ of $\Delta$. This can be expressed as $d\lambda_i=(uB+tD)d\Delta (Cv+Dw)$, where $u$/$v$ are the left/right eigenvectors associated to $\lambda_i(A_0)$, $t=uB\Delta(I-D\Delta)^{-1}$ and $w=\Delta(I-D\Delta)^{-1}Cv$. A series of perturbations $d\Delta$ are computed, which progressively move the eigenvalues of $A_0$ towards the imaginary axis. Two algorithms exist. On the one hand, a series of quadratic programming problems are solved in [22] for each $\lambda_i(A_0)$. The uncertainty with minimum Froebenius norm which brings $\lambda_i(A_0)$ on the imaginary axis is first determined, and then the one with minimum $\overline{\sigma}$ norm such that the interconnection remains at the limit of stability. On the other hand, [1] notes that the power algorithm of [2] is quite fast, but often suffers from convergence problems when purely real uncertainties are considered. A three-step procedure is proposed in this case. The initial problem is first regularized by adding a small amount $\epsilon\ll 1$ of complex uncertainty, i.e. $\Delta$ is replaced with $\Delta+\epsilon\Delta_C$, where $\Delta_C$ has the same structure as $\Delta$ but can take on complex values [23]. The power algorithm is then run at each point of a rough frequency grid, usually with good convergence properties. Using the resulting uncertainties as an initialization, a series or linear programming problems are finally solved, so as to find the value $\hat{\Delta}$ of $\Delta$ with minimum $\overline{\sigma}$ norm for which one of the eigenvalues of the initial interconnection becomes unstable. The resulting $\mu$ lower bound is equal to $\overline{\sigma}(\hat{\Delta})^{-1}$.

More recently, a direct approach has been proposed by [24,25]. It considers the optimization problem:

which can be recast as a sequential quadratic programming problem and solved for example using the Matlab Optimization toolbox. In practice, a relaxed condition is considered to avoid convergence issues: the equality constraint is replaced with $|\lambda_{max}(A_0)|\le\epsilon$, where $\epsilon$ is a small user-defined threshold. Problem (2) is non-convex and the result thus strongly depends on the initial value of $\Delta$. Several strategies are combined in [24] to improve the resulting $\mu$ lower bounds.

5. Optimization-based techniques

These algorithms are directly inspired by the definition of $\mu$. A first approach consists of directly solving the following non-convex optimization problem using standard nonlinear optimization tools such as the fmincon function of the Matlab Optimization toolbox:

In practice, a relaxed condition is considered to avoid convergence issues: the equality constraint is replaced with $\underline{\sigma}(I-M\Delta)\le\epsilon$ in [26,27], and with $|\textrm{det}(I-M\Delta)|\le\epsilon$ in [28], where $\epsilon$ is a small user-defined threshold. Any kind of uncertainties can be considered, but this approach is especially relevant for purely real problems, since the number of decision variables significantly increases when complex uncertainties are considered.

A formulation similar to (3) is considered in [29] in the case of (possibly repeated) real uncertainties:

\begin{equation} \min_{\scriptsize\begin{array}{c}\Delta\hspace{-0.8mm}\in\hspace{-0.8mm}{\bf\Delta}\\[-3pt]\lambda_1\hspace{-0.3mm},\hspace{-0.3mm}\lambda_2\hspace{-0.6mm}\in\hspace{-0.3mm}\mathbb{R}\end{array}}\overline{\sigma}(\Delta)+(\lambda_1^q+\lambda_2^q)\lambda\ \ \ \textrm{such that}\ \ \left\{\hspace{-1mm}\begin{array}{ll} \Re (\textrm{det}(I-M\Delta))=\lambda_1^p\\\Im (\textrm{det}(I-M\Delta))=\lambda_2^p \end{array}\right. \end{equation}

where $p$ and $q$ are odd and even positive integers respectively, and $\lambda$ is a large penalty parameter. This optimization problem is solved using the modified subgradient algorithm based on feasible values (F-MSG) introduced in [30].

Another variation is proposed in [31]. The following optimization problem is first solved for different values of $k$ using the fmincon function of the Matlab Optimization toolbox:

\begin{equation} g(k)=\min_{\Delta\in k\mathcal{B}_{{\bf\Delta}}}\Re\left(\textrm{det}(I-M\Delta)\right)\ \ \ \textrm{such that}\ \ \ |\Im\left(\textrm{det}(I-M\Delta)\right)|\le\epsilon \end{equation}

where $\epsilon$ is a user-defined threshold. A function $g$ is thus defined. A modification of the Newton-Raphson method called the secant method is then applied to determine the smallest value $\underline{k}$ of $k$ such that $g(k)=0$. A $\mu$ lower bound is finally obtained as the inverse of $\underline{k}$.

A common feature of the aforementioned optimization-based algorithms is that they all use standard nonlinear optimization tools, although the objective to be minimized is a nonsmooth function. Breakdowns might thus be encountered at points that are not local optima, because the latter are typically nonsmooth points in practice. In this context, [32] proposes a nonsmooth optimization technique to solve problem (3), whose convergence to a local minimum is ensured.

6. Geometrical approach

This method dedicated to (possibly repeated) real uncertainties combines randomization and nonlinear optimization [33]. The signs of the real and the imaginary parts of $\textrm{det}(I-M\Delta)$ are computed for randomly selected points on the surface of a given hyperbox in the uncertainty space $\mathbb{R}^N$. This hyperbox is enlarged until the four possible sign combinations are found, which means that it might contain values of $\delta_1,\dots,\delta_N$ such that $\textrm{det}(I-M\Delta)=0$. A series of contractions and expansions are then performed to approach the singular region until the size of the box becomes smaller than some tolerance value.

Remark 1. Except for the poles migration techniques, a classical strategy is to compute $\mu$ lower bounds on a predefined frequency grid. An alternative is to compute bounds on a set of frequency intervals whose union covers the whole frequency range. This can be achieved by considering frequency as an additional uncertainty $\delta_\omega$. An augmented interconnection $\tilde{M}-\tilde{\Delta}$ is obtained, where $\tilde{\Delta}=\textrm{diag}(\Delta,\delta_\omega I_m)$, $m$ is the state dimension of $M(s)$, and the frequency interval of interest is covered when $\delta_\omega\in\mathbb{R}$ varies between $-1$ and $1$ (see e.g. [28,32]). A skew-$\mu$ problem is then to be solved, since $\delta_\omega$ is bounded.

$\mu$ upper bound algorithms

Most of the algorithms which can reasonably be applied to real-world benchmarks are based on the following result [15]:

Theorem 1. Let $\beta>0$. If there exist matrices $D \in {\cal D}$ and $G \in {\cal G}$ which satisfy one of the following relations:

\begin{equation} M^*DM+j(GM-M^*G)\leq \beta^2D \end{equation} \begin{equation} \overline{\sigma}\left((I+G^2)^{-\frac{1}{4}}\left(\frac{DMD^{-1}}{\beta}-jG \right)(I+G^2)^{-\frac{1}{4}}\right)\leq 1 \ \end{equation}

where $\mathcal{D} = \{ D=D^*>0:\forall \Delta\in{\bf\Delta},\,D\Delta = \Delta D\}$ and $\mathcal{G} = \{ G=G^*:\forall \Delta\in{\bf\Delta},\,G\Delta = \Delta^* G \}$, then $\mu_{\bf\Delta}(M) \leq \beta$. The problem of minimizing $\beta$ can be solved either optimally using an LMI solver or faster but suboptimally using a gradient descent algorithm.

Computing an upper bound on $\mu_{\bf\Delta}(M(j\omega))$ on the whole frequency range of interest $\Omega$ is thus a challenging problem with an infinite number of frequency-domain constraints and optimization variables, since Theorem 1 must be applied for each $\omega\in\Omega$. In practice, it is usually solved on a finite frequency grid only (see e.g. the function mussv of the Matlab Robust Control toolbox). However, a crucial problem appears in this procedure: the grid must contain the most critical frequency point for which the maximal value of $\mu_{\bf\Delta}(M(j\omega))$ is reached. If not, the resulting lower bound on $k_r$ can be over-evaluated, i.e. be larger than the real value of $k_r$. Unfortunately, the aforementioned critical frequency is usually unknown! To overcome this issue, several approaches have been proposed to compute $\mu$ upper bounds which are guaranteed on the whole frequency range (see e.g. [34,35,36,37,38]). Most of them first use the characterizations of Theorem 1 to compute bounds at given frequencies. A hamiltonian-based technique is then implemented to determine the frequency intervals on which these bounds remain valid. A guaranteed $\mu$ upper bound is finally obtained as soon as the union of all intervals covers the whole frequency range of interest.

Several other techniques exist to compute $\mu$ upper bounds (see e.g. [3,39,40,41]) Nevertheless, they usually require a very high computational time when realistic applications are considered, so they are not detailed here. Note also that the above enumeration is not exhaustive, since this page mainly focuses on $\mu$ lower bound algorithms.

Testing framework

Considered algorithms 

The most relevant  lower bound algorithms are summarized in Table 1 and compared below. The first two are implemented in the function mussv of the standard Matlab Robust Control toolbox, which is called with the optional argument 'f' and 'fg' respectively. The fourth one is implemented in the function mulb_nreal of the SMART library and is used with default values of the optional parameters. It is based on a function already available in the Skew Mu toolbox for Matlab [42], which has been improved in [43]. For the other five algorithms, Matlab functions provided by the respective authors are used, with default values of the optional parameters in most cases (some slight modifications have been performed to handle the trade-off between accuracy and computational time). It is worth being emphasized that each of the eight considered algorithms is called in a similar fashion for all benchmarks.

Grid-based methods (1-2-6-7-8) are applied at each point of a 100-point frequency grid generated with the function gen_grid of the SMART library of the SMAC toolbox [43]. This grid is composed of 50 logarithmically-spaced points within the system bandwidth, and 50 additional points used to refine the grid in some frequency regions corresponding to weakly damped modes. Interval-based implementations are available for methods 6 and 7 (see Remark 1). For these algorithms, the system bandwidth is divided into 6 and 4 frequency intervals of equal size on a logarithmic scale respectively. A coarse 10-point frequency grid generated using the function gen_grid is used for method 4. Finally, methods 3 and 5 do not require any frequency grid. 

Some $\mu$ lower bound algorithms mentioned in the above survey are not considered in the sequel. First, the exponential-time methods cannot be applied to all benchmarks with a reasonable computational time, so only polynomial-time algorithms are compared. Then, the optimization-based technique of [19] has a very limited range of application and is omitted, since it can only be applied to complex uncertainties. Moreover, [18] shows that the power algorithm is usually more efficient to solve (1) in terms of accuracy and computational time. Finally, the approaches of [29,31] are very similar to the one of [26,27,28] which is already included in the comparison.

Note also that some techniques could be used in addition to any of the lower bound algorithms : the $\mu$-sensitivities introduced in [44] allow to reduce the number of uncertainties and thus the computational time at the price of a loss of accuracy, while branch-and-bound schemes [45] allow to tighten the gap between $\mu$ upper and lower bounds at the price of an increase in the computational time. These techniques are not considered here.

For each benchmark, a $\mu$ upper bound is computed using the algorithm proposed in [34], improved in [35,36] and implemented in the SMART library [43]. It provides an upper bound on the gap between the true value of $k_r$ and the best $k_r$ upper bound obtained with the $\mu$ lower bound algorithms listed in Table 1.

List of benchmarks

As shown in Table 2, a total of 36 challenging benchmarks are considered in this work, corresponding to various fields of application, system dimensions and structures of the uncertainties. Some of them contain poorly damped modes, which usually produce extremely sharp peaks on the $\mu$ plot, while others are characterized by large state vectors as well as numerous and/or highly repeated uncertainties. Purely real uncertainties are majority (benchmarks 1-29), since they are more common in engineering problems than complex or mixed ones (benchmarks 30-36). Moreover, the presence of complex uncertainties usually simplifies the computation of $\mu$ upper and lower bounds, making them of less interest in the present study.

Detailed numerical results

Computations are performed using Matlab R2010b on a Windows 7 Workstation with a CPU Intel Xenon W3530 running at 2.8 GHz and 6 GB of RAM. All results are shown in Table 3. Two lines are used for each benchmark to display first the lower bounds $\tilde{\mu}$ and then the corresponding computational times in seconds. A lower bound $\tilde{\mu}$ is considered to be valid if a frequency $\tilde{\omega}$ and a perturbation $\tilde{\Delta}$ of norm $\tilde{\mu}^{-1}$ are found, which satisfy $|\textrm{det}(I-M(j\tilde{\omega})\tilde{\Delta})|<10^{-8}$. The highest lower bound $\check{\mu}$ out of all algorithms is displayed in red. Bounds which are less than 1% and 5% lower than $\check{\mu}$ are highlighted in green and yellow respectively. NA means that the considered algorithm is not applicable. The symbols (g) and (i) denote the grid-based and the interval-based versions of the algorithms respectively. Finally, the gap between $\check{\mu}$ and the upper bound $\hat{\mu}$ is computed as $(\hat{\mu}-\check{\mu})/\check{\mu}$.

Image

Analysis of the results

Purely real problems 

All algorithms are applied to benchmarks 1-29. In Table 4, the mean computational time is computed for each algorithm over all benchmarks except the biochemical network (benchmark 29), which is significantly more complicated than the others and requires a high computational effort. Figure 1 must be interpreted as follows. Let $(x,y)$ be the coordinates of any star sign. Then $y$ corresponds to the number of benchmarks for which the gap $(\check{\mu}-\tilde{\mu})/\check{\mu}$ (upper plot) or the computational time (lower plot) is lower than $x$.

Image

Image

Image

The most relevant algorithm is the poles migration technique of [1], with by far the highest accuracy (less than 1% on average for all 29 benchmarks) and also the lowest computational time (less than 1 second on average for benchmarks 1-28 and 47 seconds for the biochemical network). The poles migration technique of [30] gives very good results too for all except 3 benchmarks, and the computational time remains reasonable (but still 20 times higher). This shows that poles migration is probably the most efficient way to compute accurate $\mu$ lower bounds for purely real problems. Moreover, an attractive feature of this approach is that frequency is not fixed, which allows to detect critical frequencies corresponding to peak values of $\mu$ more easily. The gain-based algorithm of [20], the direct nonlinear optimization technique (grid-based version) of [28] and the direct nonsmooth optimization technique (interval-based version) of [32] also give quite satisfactory results in most cases, since the accuracy is about 10% on average. Nevertheless, these algorithms are much slower (especially the last two). Note also that the nonsmooth technique gives the best results in most cases but sometimes proves very conservative, whereas the nonlinear one usually produces bounds which are not the best ones but yet good ones. Moreover, interval-based versions allow to drastically reduce the computational time, but do not bring any significant improvement in terms of accuracy. Then, the poles migration technique of [22] exhibits both a reasonable accuracy (less than 20% on average for all 29 benchmarks) and a very low computational time (about 1.0 seconds on average for benchmarks 1-28 and only 63 seconds for the biochemical network). Finally, the power algorithm of [2] often suffers from convergence problems, especially for frequencies corresponding to peak values of $\mu_{\bf\Delta}(M(j\omega))$. This is in stark contrast to the complex and mixed cases for which the convergence properties are excellent.

Remark 2. The gain-based algorithm of [20], the direct optimization-based techniques of [28,32] and the geometrical approach of [33] use randomization. Therefore, several runs of these algorithms might provide different $\mu$ lower bounds.

Purely complex and mixed real/complex problems

The poles migration technique of [1] and the geometrical approach of [33] cannot be applied in the presence of complex uncertainties. Moreover, the gain-based algorithm of [20] is not relevant when purely complex uncertainties are considered, since it reduces to the power algorithm of [2], and no improvement over the power algorithm is observed for mixed real/complex problems. Therefore, it is not considered in Table 5 and Figure 2.

Image

Image

Image

The most relevant algorithm is the power algorithm of [2], with the highest accuracy and almost the lowest computational time. This confirms its good convergence properties when mixed or complex uncertaintes are considered. The direct nonlinear optimization technique (grid-based version) of [28] also gives very accurate results, but the mean computational time is 2400 times higher. More generally, optimization-based techniques (algorithms 5-6-7) are not very attractive in terms of computational time, since full complex blocks generate a large number of optimization variables ($2n_i^2$ variables for a $n_i\times n_i$ block).

Approaching the exact value of $\mu$

Purely real problems

The results presented above show that the poles migration technique of [1] is the most relevant method for purely real problems. The best $\mu$ lower bounds are obtained for 26 benchmarks out of 29. Moreover, the bounds computed for the other 3 benchmarks are not very far from the highest values over all algorithms. Improving the poles migration technique to obtain the best results in all cases does not seem to be a trivial issue. Here, the most promising idea is to combine different algorithms to get the most out of them. Apart from the poles migration techniques, the most efficient methods are the gain-based algorithm of [20] and the global optimization tools of [28] and [32]. The following three-step strategy is thus proposed:

• The poles migration technique of [1] (algorithm 4) is executed first.
• The gain-based algorithm of [20] (algorithm 2) is then executed for a few selected frequencies only using the previous results as an initialization. These frequencies can be those for which the uncertain system becomes unstable after algorithm 4 is applied or if a single uncertainty is considered. The frequency for which the lower bound on $k_r$ (i.e. the highest $\mu$ upper bound) is obtained is also a good candidate.
• A global optimization tool such as particle swarm optimization is finally applied using the previous results as an initialization. This kind of techniques offers a bunch of tuning parameters (number of swarms, of particles in each swarm, of topologies...), which allows to easily handle the trade-off between accuracy and computational time.

The improved $\mu$ lower bounds so obtained for benchmarks 20/26/29 are shown in Table 6. They are compared with the bounds computed in Table 3 by independently applying all algorithms one at a time.

Image

In all cases, the best $\mu$ lower bound is obtained with the proposed strategy. The computational time is of course larger than if algorithm 4 is applied alone. Nevertheless, it remains much smaller than with algorithms 2 and 7, which produced the best results so far. The bound is even improved for benchmarks 26 and 29, for which a single application of any algorithm cannot bring the highest value.

Taking the above results into account, the gap between the best $\mu$ lower bound and the $\mu$ upper bound is 12.71% on average for the 29 purely real problems under consideration. An important issue is now to determine which of the bounds is responsible for this. For 19 benchmarks out of 29, the gap is less than 0.00%. For the other 10, the branch-and-bound algorithm implemented in the function mubb_mixed of the SMART library [62] is applied to improve the $\mu$ upper bound (sometimes with a very high computational cost). Results are presented in Table 7.

Image

The mean value of the gap over all 29 benchmarks is now only 0.95%. It can thus be concluded that the $\mu$ lower bound computed with the poles migration technique of 4 (eventually combined with the above three-step strategy) almost always equals the exact value of $\mu$. Moreover, the computational time remains very reasonable. This strategy is implemented in the function mulb_nreal of the SMART library.

Purely complex and mixed real/complex problems

The results presented above show that the power algorithm of [2] is the most relevant technique for purely complex and mixed real/complex problems. Nevertheless, the resulting $\mu$ lower bounds are conservative in some cases: the power algorithm is slightly outperformed by other techniques for benchmarks 33/34/35. A careful investigation of several benchmarks including those presented in this page shows that the main reasons for this are threefold:

• The power algorithm is applied at each point of a fixed frequency grid. Unless by some stroke of luck, the critical frequency corresponding to the exact value of $\mu$ is usually not part of this grid, even if the latter is very tight.
• The power algorithm almost always converges in the presence of complex uncertainties, but it sometimes requires a quite large number of iterations.
• Optimization problem (1) is non-concave. Thus, the accuracy of the resulting $\mu$ lower bound strongly depends on the way the power algorithm is initialized.

In this context, the following strategy can be implemented to improve the results displayed in Table 3. First, the power algorithm of [2] is classically applied at each point of a fixed frequency grid. In most cases, the $\mu$ plot does not exhibit very sharp peaks in the presence of complex uncertainties. Thus, a reasonably rough grid (e.g. 20 frequency points) is usually enough to capture approximately the main peaks. Then, the grid is gradually tightened around the frequencies of these peaks until improvement in the $\mu$ lower bound becomes marginal. This strategy allows to address the first two issues highlighted before. Indeed, the grid becomes very tight in the most critical frequency regions and the frequency corresponding to the exact value of $\mu$ can be detected with a very good accuracy. Moreover, repeteadly applying the power algorithm with a limited number of iterations (e.g. 100 iterations) at very close frequencies is roughly equivalent to applying it once with a large number of iterations, provided that it is initialized every time with the result obtained at the previous frequency. Finally, in order to address the third issue, the power algorithm is not only applied with the latter initialization, but also with one or more random initializations at each considered frequency.

The improved $\mu$ lower bounds so obtained for benchmarks 30-36 are shown in Table 8. They are compared with the ones computed by applying the standard power algorithm of the function mussv of the Robust Control toolbox on a fixed frequency grid. The $\mu$ upper bounds are the same as in Table 3 except for benchmark 33, where the branch-and-bound algorithm implemented in the function mubb_mixed of the SMART library has been used to reduce conservatism.

Image

The exact value of $\mu$ is obtained in all cases. The mean computational time is 1.1s. It is exactly the same as with mussv (see Table 5). Indeed, the initial frequency grid is quite rough and it is only tightened in a few selected frequency intervals. Moreover, the number of power iterations performed at each frequency is quite low. The total computational time thus remains very reasonable. The improved version of the power algorithm is implemented in the function mulb_mixed of the SMART library.

Remark 3. Several techniques have been proposed to improve the convergence properties of the power algorithm in case of purely real problems [16,17,18]. But convergence is usually not an issue in the presence of complex uncertainties. These techniques are thus of little help in the present context.

Extension to skew-$\mu$ lower bound computation

Let us now focus on the computation of $k_s$ upper bounds, i.e. of skew-$\mu$ lower bounds (see this page for more details). Most of the $\mu$ lower bound algorithms can be extended more or less easily to handle skew uncertainties. The objective here is not to study these generalizations in details, but to show that the most relevant methods highlighted above still perform very well in the skew case.

Purely complex and mixed real/complex problems

The most efficient $\mu$ lower bound algorithm is by far the power algorithm of [2]. It can be easily extended as explained in [63,64,65]. Combined with the strategy detailed in the previous section, it usually provides very good $\nu$ lower bounds. Indeed, a significant number of tests have revealed that the exact value of $\nu$ can be obtained in many cases. This improved version of the power algorithm is implemented in the function mulb_mixed of the SMART library.

Purely real problems

The most efficient $\mu$ lower bound algorithm is the poles migration technique of [1]. It can be easily extended to skew-$\mu$ problems and it is implemented in the function mulb_nreal of the SMART Library. Concerning the nonlinear and the nonsmooth optimization techniques of [28] and [32], Matlab functions compatible with skew uncertainties have been developed by the respective authors, which allows to perform a comparison in the next subsection. It would have been interesting to consider also some extended versions of the gain-based algorithm of [20] and the poles migration technique of [24], since satisfactory results are obtained in Table 3 in the classical $\mu$ case. Unfortunately, no Matlab implementations dedicated to skew-$\mu$ problems are available.

Numerical results

This section focuses on purely real problems. Methods 1-4-6-7 of Table 1 are compared on two realistic applications proposed in [32]:

• The first one describes the longitudinal behavior of a bank-to-turn missile. The nominal closed-loop system $M(s)$ has 6 states, while $\Delta$ contains 4 non-repeated real parametric uncertainties, half of which being fixed range. More details as well as numerical data can be found in [32,50].
• The second one describes the lateral behavior of a passenger aircraft. $M(s)$ has 46 states, while $\Delta$ contains 20 non-repeated real parametric uncertainties, 14 of which being fixed range. More details can be found in [54] and numerical data is available in the Skew Mu toolbox for Matlab [42].

Results are shown in Table 9. Two lines are used for each benchmark to display first the lower bounds $\tilde{\nu}$ and then the corresponding computational times in seconds. A lower bound $\tilde{\nu}$ is considered to be valid if a frequency $\tilde{\omega}$ and a perturbation $\tilde{\Delta}=\mathrm{diag}(\tilde{\Delta}_1,\tilde{\Delta}_2)$ such that $\overline{\sigma}(\tilde{\Delta}_1)\le 1$ and $\overline{\sigma}(\tilde{\Delta}_2)=\tilde{\nu}^{-1}$ are found, which satisfy $|\textrm{det}(I-M(j\tilde{\omega})\tilde{\Delta})|<10^{-8}$. The symbols (g) and (i) denote the grid-based and the interval-based versions of the algorithms respectively. Finally, the gap between the highest lower bound $\check{\nu}$ and the upper bound $\hat{\nu}$ is computed as $(\hat{\nu}-\check{\nu})/\check{\nu}$.

Image

Once again, the most relevant algorithm is the poles migration technique of [1], with both the highest accuracy and the lowest computational time. Indeed, the exact value of $\nu$ is obtained in both cases in only 3.8s on average. The nonsmooth optimization technique of [32] produces mixed results, since the exact value of $\nu$ is obtained for the first benchmark only. Note that much better results are reported in [32] for the second benchmark. Nevertheless, this requires to consider a much larger number of frequency points or intervals, and the computational burden significantly increases. As expected, the power algorithm of [63] fails to converge for a large number of frequencies and returns poor results.

Conclusion

This page provides a thorough comparative analysis of all practical methods to compute lower bounds on the structured singular value $\mu$. The most relevant algorithms are the poles migration technique of [1] for purely real problems and the power algorithm of [2] for purely complex and mixed real/complex problems. Based on these conclusions, simple improvements and combinations are proposed to approach the exact value of $\mu$ with a reasonable computational effort, resulting in state-of-the-art tools. Their evaluation on a wide set of 36 real-world benchmarks available in the litterature shows that non-conservative values can be obtained in almost all cases. Moreover, their extension to handle skew-$\mu$ problems also proves conclusive. All these tools are implemented in the Skew Mu Analysis based Robustness Tools (SMART) library of the SMAC toolbox.

References