clear;

s = tf('s');

k = ureal('k', 2.5, 'Range', [2, 3]);
tau = ureal('tau', 2.5, 'Range', [2, 3]);
theta = ureal('theta', 2.5, 'Range', [2, 3]);

% delay is approximated using pade:
% e^(-theta*s) = e^(-theta*s)/2      1 - theta*s/2
%                -------------  =  --------------
%                e^(theta*s)/2       1 + theta*s/2

G = k/(tau*s+1) * (1-theta*s)/2  / ((1+theta*s)/2);
G_nom = k.NominalValue/(tau.NominalValue*s+1);

% G = (1+I_m)G_nom => I_m = (G-G_nom)/G_nom
I_m = (G-G_nom)/G_nom;

% get wcgain
w = logspace(-3, 3, 40);
[wcg, wcu, info] = wcgain(I_m, w);

% fit worst-case magnitude to obtain a closed form of the multiplicative
% uncertainty
data = frd(info.Bounds(:, 2), w);
% use a second-order model
I_m_analytic = tf(fitmagfrd(data, 2));

% For MIMO systems, the division in the expression for I_m doesn't make
% sense
% This is another way to compute it:

% generate some instances of the uss
systems = usample(G, 20);
data = frd(systems, w);

% cover the frequency response data using an order 2 multiplicative
% uncertainty
% ucover returns [G, Info] where G is the uncertain system given by
% G(s) = [I+W1(s) Delta(s)] G_nom(s)
% W1(s) can be recovered from Info.W1
[G2, Info] = ucover(data, G_nom, 2);
I_m_analytic_2 = tf(Info.W1);

sigma(I_m, w); hold;
semilogx(w, 20*log10(info.Bounds(:, 2)), 'green');
sigma(I_m_analytic, 'red');
sigma(I_m_analytic_2, 'yellow');

I_m_analytic