%% Mu-synthesis for a distillation column in LV configuration

%% Linearized model of a distillation column in LV configuration
% Taken from Skogestad's textbook (the second edition), page 100, 
% section 3.7.2.

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

%% Weighting filters, both for the uncertainty and for the performance specs 

wp = (s/2+0.04)/(s+0.04*1e-3);                              % weighting filter for the performance specs
wi = (s+0.2)/(0.5*s+1);                                     % weighting filter for 20% relative (multiplicative) uncertainty
delta1 = ultidyn('delta1',1); delta2 = ultidyn('delta2',1); % delta blocks for both inputs
G = G0*([1,0;0,1]+wi*[delta1,0;0,delta2]);                  % uncertain system with two multiplicative uncertainties at the input

%% Building the generalized plant P
% First, set the number of measured and control variables:

nmeas = 2; % the number of measured variables
ncon = 2;  % the number of control variables

%%
% Now, build the generalized plant P

I = eye(2,2); 
P = [wp*I; I]*[I G];

%%
% Note however, that P is now still containing the uncertainty. If we were to proceed "manually", 
% we would have to pull the uncertainty out of P. In order words, we would have to express the model 
% in the form of an upper LFT of the P (void of uncertainty) and the uncertainty block. However, 
% Robust Control Toolbox can do this procedure automatically. Hence we can keep the uncertainty in P. 

%% mu-synthesis through DK iterations 
% First we set some parameters for the solver. This is not crucial, though.

options = dkitopt;
options.DisplayWhileAutoIter='on';
options.AutoIter = 'on';
options.NumberOfAutoIterations=15;
%options.AutoScalingOrder=3;

%%
% and now we are able to call the solver:

[K,clp,bnd,dkinfo] = dksyn(P,nmeas,ncon,options); % DK iterace 

%% Robustness analysis
% There is also some functionality in Robust Control Toolbox for such
% analysis (for example, now already depreceated robustperf, or its
% successor robgain) but we do the analysis "manually" by checking the
% (structured) singular values of the relevant closed-loop transfer functions. 

%%
% We first form all the possible closed-loop transfer functions. In
% particular, however, we are interested in S and T:

CL = loopsens(G,-K); % note that the computed controller considers a positive feedback, 
                     % that is why the '-'.
S = CL.So; T = CL.To;

%%
% Now we plot the sensitivity and complementary sensitivity functions for
% samples of the uncertain system.

figure(1)
sigma(S,1/wp); 
title('Sensitivity function of the uncertain system'); xlabel('Frequency'); ylabel('Gain');
grid on;

figure(2)
step(T)
title('Unit step response of the uncertain system'); xlabel('Time'); ylabel('Output');
grid on;

%%
% Obviously, both from the plots and by noticing that the peak in the achieved structured singular value
% was above 1, we can conclude that the conditions of robust stability and
% performance WERE NOT SATISFIED. Nonetheless, note that even stating these
% conditions is a matter of engineering iteration. With such a minor
% breaking of the conditions (peak in mu is 1.073 after 15 iterations), 
% we can actually conclude that the design was successful.