% Copyright 1986-2004 The MathWorks, Inc. 

%% 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 the Robust Control Toolbox.
% To estimate, e.g.,  the robust stability margin over the entire 
% frequency range, |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 dense enough frequency grid.  But 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 will 
% fail to detect worst-case perturbations and give over-optimistic    
% stability margins.
%
% This demo shows how to modify the uncertainty description to avoid 
% such discontinuities and get reliable estimates of the margin of
% uncertainty.

%% How Discontinuities Can Hide Robustness Issues
% For example, 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 |robuststab|, 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

%%
% Looking at the |mussv| bounds computed by |robuststab|, 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')

%%
% Note that making the grid denser won't help. Only if you add f=1 to the
% grid will you find out true margin.

f = 1;
stabmarg = robuststab(frd(sys,f))

%% How to Modify Your Uncertainty Model to Eliminate Discontinuitities 
% The example above shows that the robust stability margin can be a discontinuous
% function of frequency, in other words, have jumps.
% One way to eliminate such jumps is to add 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 mass/spring/damper example, model |c| as a |ureal| with 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('clti',[1 1],'Bound',0.05);
sysreg = tf(1,[m cc k]);

%%
% Now recompute the robust stability margin:

[stabmarg,du,report,info] = robuststab(frd(sysreg,omega));
stabmarg
report

%%
% 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. 

%% 
% Plot the MU bounds (reciprocal of robust stability margin) 
% as a function of frequency. Now a peak clearly appears around 
% 1 rad/sec

semilogx(info.MussvBnds)
ylim([0 1])
title('Regularized \mu, 40 frequency points, 5% added dynamics');
legend('Upper bound','Lower bound')


%% Increasing the Frequency Density for Sharper Estimates.
% The sharpness of the peak of the |mussv| plot as well as its shape
% suggests that perhaps a
% denser frequency grid is necessary.  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

%%
% 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')

%%
% For the modified uncertainty model, the robust stability margin is 1, 
% 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, it makes sense
% to put more faith in the value that allows for a few percents of unmodeled
% dynamics. Either way, the stability margin for modified problem is most
% trustworthy.

%% 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) which raised questions 
% about the continuity of the robust stability margin, and the subsequent
% computational and interpretational difficulties raised by the presence
% of discontinuities.  In 1993, a
% paper appeared in IEEE Transactions on Automatic Control entitled
% "Continuity properties of the real/complex structured singular value"
% (vol 38, no 3, Packard and Pandey).  This paper describes in detail the
% consequences and interpretations of the regularization illustrated in
% this small example.  It
% also has a extensive analysis of 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 jagger.me.berkeley.edu/papers.
%
% Thanks to Barry Postma and John Glass for suggesting this example.

