%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Section 3.7.2 DISTILLATION PROCESS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%
%   The example process is from:
%   Skogestad, Morari and Doyle,``Robust Control of Ill-Conditioned
%   Plants: High-Purity Distillation'', IEEE Atomat. Control 33,
%   1092-1105 (1988).
%

%Plant (3.82): distilation column, G(s)=1/(75s+1)[87 -86.4;108.2 -109.6]
G0 = [87.8 -86.4; 108.2 -109.6];
[U,sval,V]=svd(G0);
dyn = nd2sys(1,[75 1]); Dyn = daug(dyn,dyn); G = mmult(Dyn,G0);

Lambda=rga(G0);
disp(sprintf('\nRGA(0)=[%6.1f%6.1f]\n       [%6.1f%6.1f]',Lambda'));

% Inverse-based controller (3.84)
dynk = nd2sys([75 1],[1 1.e-6],0.7); Dynk = daug(dynk,dynk);
Kinv = mmult(Dynk,minv(G0));
I2=eye(2); K = Kinv;

% Nominal
GK = mmult(G,K);
%% May perform minimal realization to get rid of the extra states
GK = sysbal(GK);
S = minv(madd(I2,GK)); T = msub(I2,S);

% With 20% INPUT uncertainty
Unc = [1.2 0; 0 0.8];
GKu = mmult(G,Unc,K);
Su = minv(madd(I2,GKu)); Tu = msub(I2,Su);
hinfnorm(S), hinfnorm(Su) % Peak of S; 1 nominally and 14.2 with uncertainty
hinfnorm(T), hinfnorm(Tu) % Peak of T; 1 nominally and 14.2 with uncertainty

% TIME  simulation
kr=nd2sys(1,[5 1]); % 5 min filter on reference change
Kr=daug(kr,kr); Tr = mmult(T,Kr);
y = trsp(Tr,[1;0],100,.1);
u = trsp(mmult(K,S,Kr),[1;0],100,.1);
Tru = mmult(Tu,Kr);
yu = trsp(Tru,[1;0],100,.1);
uu = trsp(mmult(K,Su,Kr),[1;0],100,.1);

subplot(211); vplot(y,yu,'--');  %Figure 3.12
axis([0 100 -0.2 2.7]);text(77,1.2,'y1'),text(77,0.2,'y2');title('OUTPUTS');
subplot(212);vplot(u,uu,'--'),title('INPUTS');figure(1);xlabel('TIME (min)');
