%% Generating an uncertainty description for a parametrically uncertain transfer function
% Example 7.6 from Skogestad (2nd ed.)

k0 = 2.5; t0 = 0; T0 = 2.5;                     % Nominal values
kmin = 2; kmax = 3;                             % Lower and upper bounds 
tmin = 2; tmax = 3; 
Tmin = 2; Tmax = 3;

%% Using an array of LTI objects
freqRange = logspace(-2,1,20);

G0_tf = tf(k0,[T0 1],'iodelay',t0); 
G0 = frd(G0_tf,freqRange);

k = linspace(kmin, kmax, 10);                   % gridding the parameters 
t = linspace(tmin, tmax, 10); 
T = linspace(Tmin, Tmax, 10);

ix = 0;
for ii = 1:length(k)                            % building an array of LTI
    for jj = 1:length(t)
        for kk = 1:length(T)
            ix = ix+1;
            G_tf(:,:,ix) = tf(k(ii),[T(kk) 1],'iodelay',t(jj)); 
            G = frd(G_tf,freqRange); 
            M(:,:,ix) = (G(:,:,ix) - G0)/G0; 
        end
    end
end

W1 = tf([4 0.2],[4/2.5 1]);                     % designing the W filter(s) 
W = W1*tf([1 1.6 1],[1 1.4 1]);

figure(1), 
bode(G);                                        % bode for the plant

figure(2), 
bodemag(M);                                     % relative error 
hold on,                                        %
bodemag(W1,'r'),                                % and the filters
bodemag(W,'k')

figure(3)
iw = 20;                                        % pick the index of frequency         
w = freqRange(iw);                              % the picked frequency from grid
G0atw = evalfr(G0_tf,i*w);                      % freq. resp of G0 at some frequency 
Gatw = evalfr(G_tf(:,:,iw),i*w)                 % freq. response for samples
plot(G0atw,'ro')                                % plotting the point for nominal
hold on
plot(Gatw,'.r')                                 % plotting the samples
hold off

%% Using the functionality from Robust Control Toolbox    

k = ureal('k',k0,'Range',[kmin,kmax])
t = ureal('t',k0,'Range',[tmin,tmax])
T = ureal('T',k0,'Range',[Tmin,Tmax])

G = tf(k,[T 1],'InputDelay',tmin)
Gsamples = usample(G,1000);

w = 0.1
Gatw = evalfr(Gsamples,i*w)

plot(squeeze(Gatw),'.')