%% Analysis of robust performance for a destillation column in LV configuration
% Here we show how robust performance analysis can be reformulated as a
% robust stability analysis upon introduction of a fictitious uncertainty.
% Then, of course, we need to invoke a tool for testing a robust stability with structured uncertainty. 

%% Entering a transfer function model of the plant (matrix numerator and scalar denominator)
% The model is found in Skogestad's book (the second edition) on page 322.

s = tf('s');
G0 = [87.8, -86.4; 108.2, -109.6]/(75*s+1); 

%% Design a controller 
% The best controller typically contains the inverse of the transfer
% function (matrix) of the plant. We add an integrator and a gain:

K = tf(0.7,[1 0]) * inv(tf(G0));

%% Model of uncertainty
% In this case we do not have any additional information, hence we follow
% the rule of thumb that there is always at least some 10% relative error
% at low to middle frequencies:

tau = 1.0;
r0 = 0.2;
rinf = 2.0;

%wi = (tau*s+r0)/(tau/rinf*s+1) 
%wi = ss(wi);  

wi = makeweight(r0,1/tau,rinf);

%%
% We now take advantage of having a possibility to define the two $\Delta$
% terms in Matlab and include them in the overal description of the
% uncertain system

Delta_i = [ultidyn('Delta_i_1',[1 1]), 0; 0, ultidyn('Delta_i_2',[1 1])];
G = G0*(eye(2,2)+wi*Delta_i);

%%
% This is actually not really needed, but it will make plotting the
% responses of uncertain systems easier and we will do such plotting at the end.

%% Performance specifications
% As usual, the required "shape" of the sensitivity function must be
% specified by means of a weighting filter

M = 2;
wb = 0.005;
A = 0;

wp = (s/M+wb)/(s+wb*A); wp = ss(wp);

%% Closed-loop transfer functions
% Now we close the feedback loop and find all the necessary closed-loop
% transfer functions

C = loopsens(G,K);
C0 = loopsens(G0,K); So = C0.So; Ti = C0.Ti;

%%
% Having the individual closed-loop transfer functions at hand, we can now
% build the lower LFT of the generalized plant (including the weights) and the
% feedback controller. 

N = [wi*Ti, wi*K*So; wp*So*G0, wp*So];

%% Computing the structured singular values
% The N matrix is what is now used for interconnection into an upper LFT
% with the structured uncertainty block. However, it is not only the
% original structured uncertainty $\Delta$ but we also introduce an
% auxiliary (artificial) uncertainty $\Delta_\mathrm{p}$, which captures the performance specs 

omega = logspace(-3,2,100);
uncBlckStruc = ones(2,2);
perfUncBlckStruc = [2 2];
fictUncBlckStruc = [uncBlckStruc; perfUncBlckStruc];

%%
% We now compute the structured singular values (SSV, or $\mu$) for both the $N_{11}$ term and the full $N$. 
% Note that in Matlab these are computed frequency by frequency. 

muWT = mussv(frd(wi*C0.Ti,omega),uncBlckStruc); 
muN = mussv(frd(N,omega),fictUncBlckStruc); 

%% Plotting the unstructured and the structured singular values
% Finally we can plot the singular and structured singular values. Here we
% decide to plot these for the weighted closed-loop transfer functions. We
% are then checking if the plots stay below 0 dB.

sigma(muWT(1)) 
hold on 
sigma(muN(1)) 
sigma(frd(wp*C0.So,omega)) 
hold off
legend('\mu(Wi Ti)', '\mu(N)','\sigma_1(Wp Si)')

title('');
ylabel('Magnitude')
xlabel('Frequency')
grid on

%%
% The conclusion is that while the nominal performance and robust stability
% are guaranteed, the controller fails to guarantee robust performance.