%Example 6.7 Feedback control of distillation process.
%

%Model (6.83):
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);

%1. Inverse-based control
dynk = nd2sys([75 1],[1 1.e-6],0.7); Dynk = daug(dynk,dynk); 
Kinv = mmult(Dynk,minv(G0));

% Inverse (Kinv) Check with RGA-bound
w = logspace(-2,2,81);
k=0.7; t = nd2sys(k, [1 k]); s = msub(1,t); tf=frsp(t,w); sf=frsp(s,w);
rgamaxrow=norm(vrga(G0),'inf') % 69.1
wit=vabs(mmult(0.2,tf)); witrga = mmult(rgamaxrow, wit, minv(madd(1,wit)));
smaxrga = mmult(vabs(sf), madd(1,  witrga));
%peak value is 6.81 (at about w=0.7), which is less than 14.2:
Sinv_rga=pkvnorm(smaxrga)

%Actual peak value, with uncertainty Ei=diag(0.2, -0.2).
L=mmult(G,daug(1.2,0.8),Kinv);
Sinv=minv(madd(eye(2),L));
Sinv_act=hinfnorm(Sinv) % 14.2

%Compare with worst-case performance computed with skrpmu.

%       Uncertainty weight
Wi = daug(0.2,0.2);
systemnames = 'G Kinv Wi';
inputvar = '[di(2); r(2)]';
outputvar = '[Wi; r-G]';
input_to_G = '[Kinv+di]';
input_to_Kinv = '[r-G]';
input_to_Wi = '[Kinv]';
sysoutname = 'Pmsku';
cleanupsysic = 'yes';
sysic;

omega=linspace(0.6,0.75,10); 
Mjw = frsp(Pmsku, omega);
[skm, bnds, dvec, sens, pvec, gvec] = skrpmu(Mjw, [1 1;1 1; 2 2]);
Sinv_wc=pkvnorm(skm) % 14.5

% Alternative command in new version of mu-tools: 
%   [delwc,mulow,muup]=wcperf(Mjw,[1 1;1 1; 2 2],1)
% which also gives a peak value of 14.5.
% This implies that the actual worst-case peak of S' is 14.5, so the 
% uncertainty E_= diag(0.2, -0.2), which yields a peak of 14.2 
% (see p. 242), is NOT  quite the worst case.

%2. Diagonal feedback controller.
k2=2.4e-2; tau=75;
dynk=nd2sys([tau 1],[1 1e-6],k2);
Kdiag=mmult(daug(dynk,dynk),[1 0;0 -1]);

%Norminal S and T:
S=minv(madd(eye(2),mmult(G,Kdiag)));
T=msub(eye(2),S);
gok=1;wi=0.2;
Ss=veval('max',vsvd(frsp(S,w)));
Ts=veval('max',vsvd(frsp(T,w)));

%Upper bound of uncertain ||S||, (6.77)
Subd=mmult(Ss,minv(msub(1,mmult(gok*wi,Ts))));
Sdiag_ubd=hinfnorm(Subd) %1.26

%Actual peak value, with uncertainty Ei=diag(0.2, -0.2).
L=mmult(G,daug(1.2,0.8),Kdiag);
Sdiag=minv(madd(eye(2),L));
Sdiag_act=hinfnorm(Sdiag) % 1.05

%Compare with worst-case performance computed with skrpmu.

%       Uncertainty weight
Wi = daug(0.2,0.2);
systemnames = 'G Kdiag Wi';
inputvar = '[di(2); r(2)]';
outputvar = '[Wi; r-G]';
input_to_G = '[Kdiag+di]';
input_to_Kdiag = '[r-G]';
input_to_Wi = '[Kdiag]';
sysoutname = 'Pmsku1';
cleanupsysic = 'yes';
sysic;

omega=linspace(0.9,1.1,10); 
Mjw1 = frsp(Pmsku1, omega);
[skm1, bnds, dvec, sens, pvec, gvec] = skrpmu(Mjw1, [1 1;1 1; 2 2]);
Sdiag_wc=pkvnorm(skm1) % 1.05

%Result:
clc
disp('1. Inverse-based controller:');
disp(sprintf('   Lower bound of ||S|| estimated using RGA:%6.2f',Sinv_rga));
disp(sprintf('   Actual ||S|| with E_I=diag(0.2,-0.2):    %6.2f',Sinv_act(1)));
disp(sprintf('   Worst-case performance using skewed-mu:  %6.2f',Sinv_wc));
disp('2. Diagonal controller:');
disp(sprintf('   Upper bound of ||S|| using (6.77):       %6.2f',Sdiag_ubd));
disp(sprintf('   Actual ||S|| with E_I=diag(0.2,-0.2):   %7.2f',Sdiag_act(1)));
disp(sprintf('   Worst-case performance using skewed-mu:  %6.2f',Sdiag_wc));

