Getting Reliable Estimates of Robustness Margins
Computing the smallest perturbation that causes instability at a given frequency is the cornerstone of most robustness analysis algorithms in Robust Control Toolbox™. To estimate, for example, the robust stability margin over the entire frequency range, the function robuststab performs this basic computation for a finite set of frequencies and looks at the worst case over this set.
Under most conditions, the robust stability margin is continuous with respect to frequency, so this approach gives good estimates provided you use a sufficiently dense frequency grid. However in some problems with only parameter uncertainty (UREAL), the migration of poles from stable to unstable can occur at isolated frequencies (generally unknown to you). As a result, any frequency grid that excludes these particular frequencies fail to detect worst-case perturbations and give over-optimistic stability margins.
In this demo we look at how to modify the uncertainty description to avoid discontinuities and get reliable estimates of the margin of uncertainty.
Contents
How Discontinuities Can Hide Robustness Issues
In this example, we consider a spring-mass-damper system with 100% parameter uncertainty in the damping coefficient and 0% uncertainty in the spring coefficient. Note that all uncertainty here is of ureal type:
m = 1; k = 1; c = ureal('c',1,'plusminus',1); sys = tf(1,[m c k]);
As the uncertain element c varies, the only place where the poles can migrate from stable to unstable is at s = j*1 (1 rad/sec). No amount of variation in c can cause them to migrate across the jw-axis at any other frequency.
When computing the robust stability margins with the robuststab function, almost any frequency grid will exclude f=1, leading to the incorrect conclusion that the margin of uncertainty on c is infinite.
omega = logspace(-1,1,40); % one possible grid
sysg = frd(sys,omega);
[stabmarg,du,rep,info] = robuststab(sysg);
stabmarg
stabmarg = UpperBound: Inf LowerBound: 2.2675e+015 DestabilizingFrequency: 0.1000
When we look at the mussv bounds computed by the robuststab function, the structured singular value mu is zero at all test frequencies in omega:
semilogx(info.MussvBnds,'*') ylim([0 1]) title('Original \mu bounds, 40 frequency points') legend('Upper bound','Lower bound')

Figure 1: Original mu bounds, 40 frequency points.
Note that making the grid denser would not help. Only by adding f=1 to the grid will we find the true margin.
f = 1; stabmarg = robuststab(frd(sys,f))
stabmarg = UpperBound: 1.0000 LowerBound: 1.0000 DestabilizingFrequency: 1
Modifying the Uncertainty Model to Eliminate Discontinuities
The example above shows that the robust stability margin can be a discontinuous function of frequency. In other words, it can have jumps. We can eliminate such jumps by adding a small amount of uncertain dynamics to every uncertain real parameter. This amounts to adding some dynamics to pure gains. Importantly, as the size of the added dynamics goes to zero, the estimated margin for the modified problem converges to the true margin for the original problem.
In the spring-mass-damper example, we model c as a ureal with the range [0.05,1.95] rather than [0,2], and add a ultidyn perturbation with gain bounded by 0.05. This combination covers the original uncertainty in c and introduces only 5% conservatism.
cc = ureal('cReal',1,'plusminus',0.95) + ultidyn('cUlti',[1 1],'Bound',0.05); sysreg = usubs(sys,'c',cc);
Now let's recompute the robust stability margin:
[stabmarg,du,report,info] = robuststab(frd(sysreg,omega)); stabmarg report
stabmarg = UpperBound: 2.3630 LowerBound: 2.0016 DestabilizingFrequency: 0.9427 report = Assuming nominal UFRD system is stable ... Uncertain System is robustly stable to modeled uncertainty. -- It can tolerate up to 200% of the modeled uncertainty. -- A destabilizing combination of 236% of the modeled uncertainty exists, causing an instability at 0.943 rad/s. -- Sensitivity with respect to uncertain element ... 'cReal' is 35%. Increasing 'cReal' by 25% leads to a 9% decrease in the margin. 'cUlti' is 63%. Increasing 'cUlti' by 25% leads to a 16% decrease in the margin.
Now the calculation determines that the margin is not infinite. The value 2.36 is still greater than 1 (the true margin) because the density of frequency points is not high enough.
If we plot the mu bounds (reciprocal of robust stability margin) as a function of frequency, a peak clearly appears around 1 rad/s:
semilogx(info.MussvBnds) ylim([0 1]) title('Regularized \mu, 40 frequency points, 5% added dynamics'); legend('Upper bound','Lower bound')

Figure 2: Regularized mu, 40 frequency points, 5% added dynamics
Increasing the Frequency Density for Sharper Estimates.
The sharpness of the peak of the mussv plot as well as its shape suggests that a denser frequency grid may be necessary. We'll increase the frequency-grid density by a factor of 5, and recompute:
OmegaDense = logspace(-1,1,200); [stabmarg,du,report,info] = robuststab(frd(sysreg,OmegaDense)); stabmarg report
stabmarg = UpperBound: 1.0056 LowerBound: 1.0026 DestabilizingFrequency: 1.0116 report = Assuming nominal UFRD system is stable ... Uncertain System is robustly stable to modeled uncertainty. -- It can tolerate up to 100% of the modeled uncertainty. -- A destabilizing combination of 101% of the modeled uncertainty exists, causing an instability at 1.01 rad/s. -- Sensitivity with respect to uncertain element ... 'cReal' is 95%. Increasing 'cReal' by 25% leads to a 24% decrease in the margin. 'cUlti' is 6%. Increasing 'cUlti' by 25% leads to a 2% decrease in the margin.
The computed robust stability margin has decreased significantly. It helped to use a denser grid, and the margin estimate is now close to the true value 1. The mu plot shows a pronounced peak at 1 rad/sec:
semilogx(info.MussvBnds) ylim([0 1]) title('Regularized \mu, 200 frequency points, 5% added dynamics'); legend('Upper bound','Lower bound')

Figure 3: Regularized mu, 200 frequency points, 5% added dynamics
For the modified uncertainty model, the robust stability margin is 1, which is equal to that of the original problem. In general, the stability margin of the modified problem is less than or equal to that of the original problem. If it is significantly less, then the answer to the question "What is the stability margin?" is very sensitive to the uncertainty model. In this case, we put more faith in the value that allows for a few percents of unmodeled dynamics. Either way, the stability margin for the modified problem is most trustworthy.
Automating Subsitution of Uncertainty Models
The command complexify automates the procedure of replacing a ureal with the sum of a ureal and ultidyn. The analysis above can be repeated using complexify obtaining the same results.
sysreg = complexify(sys,0.05,'ultidyn'); [stabmargc,duc,reportc,infoc] = robuststab(frd(sysreg,OmegaDense)); stabmargc reportc semilogx(info.MussvBnds,'b',infoc.MussvBnds,'r') ylim([0 1]) title('Regularized \mu, 200 frequency points, 5% added dynamics'); legend('Upper bound','Lower bound','(complexify) Upper bound',... '(complexify) Lower bound','Location','NorthEast');
stabmargc = UpperBound: 1.0056 LowerBound: 1.0026 DestabilizingFrequency: 0.9885 reportc = Assuming nominal UFRD system is stable ... Uncertain System is robustly stable to modeled uncertainty. -- It can tolerate up to 100% of the modeled uncertainty. -- A destabilizing combination of 101% of the modeled uncertainty exists, causing an instability at 0.988 rad/s. -- Sensitivity with respect to uncertain element ... 'c' is 95%. Increasing 'c' by 25% leads to a 24% decrease in the margin. 'c_cmpxfy' is 6%. Increasing 'c_cmpxfy' by 25% leads to a 2% decrease in the margin.

Figure 4: Regularized mu, 200 frequency points, 5% added dynamics
References
In 1989, a paper appeared in Systems and Control Letters, entitled "Robustness Margin need not be a continuous function of the problem data" (vol 15, page 91-98, Barmish, Khargonekar, Shi and Tempo) raised questions about the continuity of the robust stability margin, and the subsequent computational and interpretational difficulties raised by the presence of discontinuities. A 1993 paper in IEEE® Transactions on Automatic Control entitled "Continuity properties of the real/complex structured singular value" (vol 38, no 3, Packard and Pandey) describes in detail the consequences and interpretations of the regularization illustrated in this small example. It also extensively analyzes the regularization as applied to the interesting 2-parameter example from Barmish, et., al. A pdf file of the IEEE TAC paper is available at http://jagger.me.berkeley.edu/publications.php
Thanks to Barry Postma and John Glass for suggesting this example.