% The following simulations confirm the frequency plot in Fig.10.7.
% Results by: Kjetil Havre, July 1996

clear all
close all

% 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 = [-6.292e-2 6.236e-2; 5.481e-3 -1.719e-2; 3.041e-3 -1.078e-2;
     -1.856e-2 -1.393e-2; -1.229e-1 -5.608e-3]

C= [-7.223 -5.170 3.386 -1.633e-1 1.121;
    -8.913 4.728 9.876 8.425 2.186]

D = [ 0 0; 0 0]

Bd = [-9.364e-3  -1.333e-2 ; 1.960e-2  8.018e-3 ; 3.266e-3 -2.116e-2;
      -2.827e-2  5.319e-3; -6.784e-3  2.719e-3]

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
Dmag = [0.2/0.3 0; 0 0.2/0.4];
Gd = mmult(Gd,Dmag);

% Simple PI controller
Ti = 2; Kp = -1;
K2 = nd2sys([Ti 1], [Ti 1e-6], Kp);

G22 = sel(G,2,2); w = logspace(-3, 1, 201);
L22 = mmult(G22,K2);
L22jw = frsp(L22,w);
vplot( 'liv,lm', L22jw, 1, ':');

Kf = nd2sys([1], [3 1], 0.54);
%
% Note that there is minus error in Eq. (10.35)
%

[A,B,C,D] = unpck(G);
[Ad,Bd,Cd,Dd] = unpck(Gd);
Gs = pck(A,[B Bd], C, [D Dd]);

%
% Remark:
%  If one forms the combined system consisting of G and Gd  
%  so that the A matrix is duplicated, that is, the states 
%  are taken into account twice it is diffucult to remove 
%  extra states after the partially controlled system is 
%  formed, resulting in numerical problems.
%

systemnames     = 'Gs K2 Kf';
   inputvar     = '[r2(1); d(2)]';
   outputvar    = '[Gs; r2-Gs(2); Kf; K2]';
   input_to_Gs  = '[Kf; K2; d]';
   input_to_K2  = '[r2-Gs(2)]';
   input_to_Kf  = '[d(1)]';
   sysoutname   = 'P';
   cleanupsysic = 'yes';
sysic;

%P = sminreal(P, 1e-1);
if isstable(P)
   P = sysbal(P);
else
   disp('Partiall controlled system is unstable');
end

spoles(P)

% Reference change
y1 = trsp(P, [1 0 0]', 50, 0.1);
figure(1); vplot(sel(y1,[1 2], 1)); 
ylabel('Compositions %'); xlabel('Time [min]')
figure(2); vplot(sel(y1,[4 5],1)); 
ylabel('Controls'); xlabel('Time [min]')

% Feed flowrate disturbance
y2 = trsp(P, [0 1 0]', 50, 0.1);
figure(1); vplot(sel(y2,[1 2],1)); 
ylabel('Compositions %'); xlabel('Time [min]')
figure(2); vplot(sel(y2,[4 5],1)); 
ylabel('Controls'); xlabel('Time [min]')

% Feed composition disturbance
y3 = trsp(P, [0 0 1]', 50, 0.1);
figure(1); vplot(sel(y3,[1 2],1)); 
ylabel('Compositions %'); xlabel('Time [min]')
figure(2); vplot(sel(y3,[4 5],1)); 
ylabel('Controls'); xlabel('Time [min]')

