%Example 10.10

% START: LINEAR MODEL
% Distillation column model with 5 states and disturbances
A = [ -5.131e-3 0 0 0 0; 0 -7.366e-2 0 0 0; 0 0 -1.829e-1 0 0;
       0 0 0 -4.620e-1  9.895e-1; 0 0 0 -9.895e-1 -4.620e-1]
B = [-.629 .624; .055 -0.172; 0.030 -0.108; 
      -0.186 -0.139; -1.23 -0.056]
C= [-0.7223 -0.5170 0.3386 -0.1633e-1 0.1121;
    -0.8913 0.4728 0.9876 0.8425 0.2186]
D = [ 0 0; 0 0]
Bd = [-0.062 -0.067; 0.131 0.040; 0.022 -0.106; 
     -0.188 0.027; -0.045 0.014]
G = pck(A,B,C,D);
Gd = pck(A,Bd,C,D);
% Change the disturbance magnitudes from 30 to 20% for d1 and
% from 40 to 20% for d2 (this is because the scalings used in this
% example are somewhat different from that used in the paper by
% Hovd and Skogestad (Automatica, 1992) from which the model above
% was obtained).
Dmag = [0.2/0.3 0; 0 0.2/0.4];
Gd = mmult(Gd,Dmag);
[Ad,Bd,Cd,Dd]=unpck(Gd);
% END: LINEAR MODEL

zeros=tzero(A,B,C,D)

% Steady-state analysis
G0 = -C*inv(A)*B + D
Gd0 = -Cd*inv(Ad)*Bd + Dd
RGA0=vrga(G0)
PRGA0 = mmult(vdiag(vdiag(G0)),minv(G0))
CLDG0 = mmult(PRGA0,Gd0)
GinvGd0 = mmult(inv(G0),Gd0)


w = logspace(-3,1,41);
Gf = frsp(G,w); Gdf=frsp(Gd,w);
%Figure 10.9  
figure(1);clf;
vplot('liv,lm',Gdf,1,':');
axis([0.001,10,.1,20]);
text(.08,0.2,'11');
text(.117,0.30,'12');
text(0.15,0.4,'22');
text(0.5,0.4,'21');
xlabel('Frequency');ylabel('Magnitude');
title('Elements in Gd')
drawnow;
%vplot('liv,lm', mmult(minv(Gf),Gdf)); title('Ginv*Gd')

rgag=vrga(Gf);
gdiag=vdiag(vdiag(Gf));
prga = mmult(gdiag,minv(Gf));
cldg=mmult(prga,Gdf);

%Figure 10.10
figure(2);clf
vplot('liv,lm',cldg,1,':');
axis([0.001,10,.1,100]);
xlabel('Frequency');ylabel('Magnitude');
title('CLDG-elements');
text(0.005,20,'11'); text(0.005,0.4,'12');
text(0.005,5,'22'); text(0.02,25,'21');

%Figure 10.11
figure(3);clf
vplot('liv,lm',prga,1,':');
axis([0.001,10,.1,100]);
xlabel('Frequency');ylabel('Magnitude');
title('PRGA-elements');
text(4,1.8,'11=22'); text(0.1,1.3,'12');
text(0.6,2,'21');

k1 = nd2sys([3.76 1], [3.76 1.e-4],  0.261);
k2 = nd2sys([3.31 1], [3.31 1.e-4], -0.375);
K = daug(k1,k2);
GK = mmult(G,K);
S = minv(madd(eye(2),GK)); SGd = mmult(S,Gd);
Sf=frsp(S,w);
d = transp(vpck([1 0; 1 1],[0; 50]));
y = trsp(SGd,d,100,0.2);
%Figure 10.12
figure(4);clf
vplot(y);
xlabel('Time');
text(30,-0.25,'y2'); text(30., 0.17,'y1')
