% Shell Model relating y1,y2,y7 to u1 u2 u3 l1 and l2.

% Created by B. Joseph April 2001

%

e1=-1.0;e2=-1.0;e3=-1.0;e4=1.0;e5=1.0; %  Fractional error in gains

% Uses kij  to represent gains, tij to represent time constants and dij to represent delays


% u1=top draw, u2=side draw, u3=botttoms duty, u4= reflux duty, u5= upper reflux duty


% y1=top end point, y2=side end point, y3= bottoms reflux temp


% gij represents transfer function relating input i and output j.


%


% Defining the transfer function of the process
delt = 0;
num = 4.05+2.11*e1;
den = [50 1];
delay = 27;  
g11 = poly2tfd(num,den,delt,delay);
%
%
num = 5.39+3.29*e1;
den = [50 1];
delay = 18;  
g12 = poly2tfd(num,den,delt,delay);
%
%
num = 4.38+3.11*e1;
den = [33 1];
delay = 20;  
g13= poly2tfd(num,den,delt,delay);
%
%
num = 1.77+.39*e2;
den = [60 1];
delay = 28;
g21 = poly2tfd(num,den,delt,delay);
%
%
num = 5.72+.57*e2;
den = [60 1];
delay = 14;  
g22= poly2tfd(num,den,delt,delay);
%
%
num = 4.42+3.11*e2;
den = [44 1];
delay = 22;
g23 = poly2tfd(num,den,delt,delay);
%
%
num = 5.88+.59*e3;
den = [50 1];
delay = 27;  
g31 = poly2tfd(num,den,delt,delay);
%
%
num = 6.90+0.89*e3;
den = [40 1];
delay = 15;
g32 = poly2tfd(num,den,delt,delay);
%
%
num = 7.2+1.33*e3;
den = [19 1];
delay = 0;  
g33= poly2tfd(num,den,delt,delay);
%
%
num = 1.2+.12*e4;
den = [45 1];
delay = 27;  
g41 = poly2tfd(num,den,delt,delay);
%
%
num = 1.52+.13*e4;
den = [25 1];
delay = 15;  
g42= poly2tfd(num,den,delt,delay);
%
%
num = 1.14+.18*e4;
den = [27 1];
delay = 0;
g43 = poly2tfd(num,den,delt,delay);
%
%
num = 1.44+.16*e5;
den = [40 1];
delay = 27;  
g51 = poly2tfd(num,den,delt,delay);
%
%
num = 1.83+.13*e5;
den = [20 1];
delay = 15;  
g52 = poly2tfd(num,den,delt,delay);
%
%
num = 1.26+.18*e5;
den = [32 1];
delay = 0;  
g53 = poly2tfd(num,den,delt,delay);
%
%
tfinal = 250; 
delt = 5;
nout = 3;
plant = tfd2step(tfinal,delt,nout,g11,g12,g13,g21,g22,g23,g31,g32,g33);
%pause
dplant = tfd2step(tfinal,delt,nout,g41,g42,g43,g51,g52,g53);

%Model definition
e1=0;e2=0;e3=0;e4=0;e5=0;
delt = 0;
num = 4.05+2.11*e1;
den = [50 1];
delay = 27;  
g11 = poly2tfd(num,den,delt,delay);
%
%
num = 5.39+3.29*e1;
den = [50 1];
delay = 18;  
g12 = poly2tfd(num,den,delt,delay);
%
%
num = 4.38+3.11*e1;
den = [33 1];
delay = 20;  
g13= poly2tfd(num,den,delt,delay);
%
%
num = 1.77+.39*e2;
den = [60 1];
delay = 28;
g21 = poly2tfd(num,den,delt,delay);
%
%
num = 5.72+.57*e2;
den = [60 1];
delay = 14;  
g22= poly2tfd(num,den,delt,delay);
%
%
num = 4.42+3.11*e2;
den = [44 1];
delay = 22;
g23 = poly2tfd(num,den,delt,delay);
%
%
num = 5.88+.59*e3;
den = [50 1];
delay = 27;  
g31 = poly2tfd(num,den,delt,delay);
%
%
num = 6.90+0.89*e3;
den = [40 1];
delay = 15;
g32 = poly2tfd(num,den,delt,delay);
%
%
num = 7.2+1.33*e3;
den = [19 1];
delay = 0;  
g33= poly2tfd(num,den,delt,delay);
%
%
num = 1.2+.12*e4;
den = [45 1];
delay = 27;  
g41 = poly2tfd(num,den,delt,delay);
%
%
num = 1.52+.13*e4;
den = [25 1];
delay = 15;  
g42= poly2tfd(num,den,delt,delay);
%
%
num = 1.14+.18*e4;
den = [27 1];
delay = 0;
g43 = poly2tfd(num,den,delt,delay);
%
%
num = 1.44+.16*e5;
den = [40 1];
delay = 27;  
g51 = poly2tfd(num,den,delt,delay);
%
%
num = 1.83+.13*e5;
den = [20 1];
delay = 15;  
g52 = poly2tfd(num,den,delt,delay);
%
%
num = 1.26+.18*e5;
den = [32 1];
delay = 0;  
g53 = poly2tfd(num,den,delt,delay);
%
%
tfinal = 250; 
delt = 5;
nout = 3;
model = tfd2step(tfinal,delt,nout,g11,g12,g13,g21,g22,g23,g31,g32,g33);
%pause
dmodel = tfd2step(tfinal,delt,nout,g41,g42,g43,g51,g52,g53);
%
%  Calculate the MPC controller gain matrix for
%    No plant/model mismatch,
%    Output Weight = 1, Input Weight = 0
%    Input Horizon = 2, Output Horizon = 20
%    model = plant;
ywt = [1 1 1]; uwt = [0 0 0];
M = 2; P = 50;
Kmpc = mpccon(model,ywt,uwt,M,P);
%
%  Simulate and plot response for unmeasured and measured
%  step disturbance through dplant.
%
tend = 400;
r = [0 0 0]; 
dstep = [-0.5 -0.5];
tfilter=[0.0 0.0 0.0];
usat=[];
u0=[0.0 0.0 0.0]
d0=[0.0 0.0]
%[y1,u1] = mpcsim(plant,model,Kmpc,tend,r,usat,tfilter,dplant,dmodel,dstep);
ulim=[-0.5 -0.5 -0.5 0.5 0.5 0.5 0.05 0.05 0.05];
ylim=[-0.5 -Inf -0.5 0.5 Inf Inf];
[y1,u1] = cmpc(plant,model,ywt,uwt,M,P,tend,r,...
   ulim,ylim,tfilter,dplant,dmodel,dstep);
plotall(y1,u1,delt)