%**************************************************************
%            A PROGRAM FOR COMBINING THE HP/LP COLUMNS
%
% - calls on colamod1.m and colamod2.m, the models for the 
%   HP and LP columns, respectively
% - The no. of stages and feed data must be changed for a 
%   new case
%
% H.K.Engelien, 30/03/2001
%**************************************************************

function xprime = multieffect(t,X,options,u0)

global Ptop_HP Ptop_LP  Y1 Y2 time Flow 

NT = 16;                            % No. stages in HP column
Nt = 26;                            % No. stages in LP column

X1 = X(1:3*NT);                     % States in HP column
X2 = X(3*NT+1:3*NT+3*Nt);           % States in LP column

%*********************** HP COLUMN ****************************

% Feed (HP column)
F_HP0  = 4320;              % Feedrate (kmol/h) 4320
zF_HP = 0.73;              % Feed composition light component
qF_HP = 1;                 % Feed liq. fraction

% Get actual values
MB_HP = X1(NT+1);          % Actual reboiler holdup (kmol)
MD_HP = X1(2*NT);          % Actual condenser holdup (kmol)
xD_HP = X1(NT);            % Concentration in distillate
xB_HP = X1(1);             % Concentration in btm product
TD_HP = X1(3*NT);          % Temperature in condenser (K)
TB_HP = X1(2*NT+1);

% Physical data (taken @ 5.2 bar, from HYSYS)
TBL   = 385.85;             % Boiling point light component (K)
TBH   = 426.45;             % Boiling point heavy component (K)
HvapL = 31630;              % Hvap for light component (kJ/kmol)
HvapH = 37870;              % Hvap for heavy component (kJ/kmol)
p0L   = 5.2e5;              % Vap. pressure for pure light liq. comp.(Pa)
p0H   = 5.2e5;              % Vap. pressure for pure heavy liq. comp.(Pa)
R     = 8.314;              % Universal gas constant (kJ/kmol*K)

% Calculating top pressure for cooling control
%pL = p0L*exp(-HvapL/R*(1/TD_HP - 1/TBL));  % Vap. pres. light comp.
%pH = p0H*exp(-HvapH/R*(1/TD_HP - 1/TBH));  % Vap. pres. heavy comp.
%P_HP =  xD_HP*pL + (1-xD_HP)*pH;           % Stage pres.(Raoult's law)

% Calculating btm pressure for cooling control
pL = p0L*exp(-HvapL/R*(1/TB_HP - 1/TBL));  % Vap. pres. light comp.
pH = p0H*exp(-HvapH/R*(1/TB_HP - 1/TBH));  % Vap. pres. heavy comp.
P_HP =  xB_HP*pL + (1-xB_HP)*pH;           % Stage pres.(Raoult's law)


% PI-controllers for Control of HP Column:
%-------------------------------------------------------------
 
D0_HP  = u0(1);                    % Nominal distillate flow (kmol/h) 
B0_HP  = u0(2);                    % Nominal bottom flow (kmol/h)
VD0_HP = u0(3);                    % Nominal flow condensated (kmol/h)
LT0_HP = u0(4);                    % Nominal reflux (kmol/h)
VB0_HP = u0(5);                    % Nominal boilup (kmol/h)

% Top (condenser) holdup control
KcD  = -60;                       
MDs  = 102.6; 
D_HP = D0_HP + (MDs-MD_HP)*KcD;  

% Bottom (reboiler) holdup control
KcB  = 4;                         
MBs  = 97.4;         
%B_HP = B0_HP + (MDs-MB_HP)*KcB;  % Uses btm. flow to control holdup
F_HP = F_HP0 + (MDs-MB_HP)*KcB;   % Uses feedrate to control holdup

% Distillate composition control
KcXD = 937500; TauI_xD = 0.0000005;     
xDs_HP = 0.99;
xD_HP_error  = (xDs_HP-xD_HP); Interr_xD = X(3*NT+3*Nt+1);
LT_HP = LT0_HP + xD_HP_error*KcXD+1/TauI_xD*Interr_xD;  

% Pressure (btm) control
KcP = 0.028;    TauI_P = 0.067;
Ps_HP  = 5.79e5;           
P_HP_error = (Ps_HP-P_HP); Interr_P = X(3*NT+3*Nt+2);
VB_HP = VB0_HP + P_HP_error*KcP + 1/TauI_P*Interr_P;

% Btm composition control
%KcXB = -210460;  TauI_xB = -0.055;      
%xBs_HP = 0.53;
%xB_HP_error = (xBs_HP-xB_HP);Interr_xB = X(3*NT+3*Nt+2);
%VB_HP = VB0_HP + xB_HP_error*KcXB + 1/TauI_xB*Interr_xB; 

%---------------------------------------------------------------------

% Integrated condenser/reboiler

Hvap_top_HP = 31710;          % From HYSYS at 5.2bar & 0.99 MeOH
Hvap_btm_LP = 40650;          % Assume pure water at 1 bar

% Fixed area of reboiler/condenser (1000 m^2)
U=3000;                   % Overall h.t. coefficient (W/m^2 C)
A=1000;                   % Heat exchanger area (m^2)
DeltaT=TD_HP-X2(2*Nt+1);  % Reboiler temperature (LP) (K)
Q=U*A*DeltaT*3600/1000;   % Reboiler/condenser duty (kJ/h)

VD_HP = Q/Hvap_top_HP;        % Condensation in HP column (kmol/h)
VB_LP = Q/Hvap_btm_LP;        % Boilup in LP column (kmol/h)

%----------------------------------------------------------------------

%************************** END HP COLUMN **************************


%*************************** LP COLUMN *****************************


% Feed (LP column - from btm stage HP column)
%F_LP  = B_HP;             % Feedrate (kmol/h)
zF_LP = xB_HP;             % Feed composition light component
qF_LP = 1;                 % Feed liq. fraction

% Get actual values
MB_LP = X2(Nt+1);          % Actual reboiler holdup (LP)
MD_LP = X2(2*Nt);          % Actual condenser holdup (LP)
xD_LP = X2(Nt);            % Concentration in distillate
xB_LP = X2(1);             % Concentration in btm product
TD_LP = X2(3*Nt);          % Temperature in condenser (LP)
T_6   = X2(2*Nt+6);        % Temperature on tray 6 (LP)

% Thermodynamic data (taken @ 1.013 bar, from HYSYS))
TBL   = 337.63;             % Boiling point light component (K)
TBH   = 373.15;             % Boiling point heavy component (K)
HvapL = 35290;              % Hvap for light component (kJ/kmol)
HvapH = 40650;              % Hvap for heavy component (kJ/kmol)
p0L   = 1.01e5;             % Vap. pressure for pure light liq. comp.(Pa)
p0H   = 1.01e5;             % Vap. pressure for pure heavy liq. comp.(Pa)
R     = 8.314;              % Univarsal gas constant (kJ/kmol*K)

% Calculating top pressure for cooling control
pL = p0L*exp(-HvapL/R*(1/TD_LP - 1/TBL));  % Vap. pres. light comp.
pH = p0H*exp(-HvapH/R*(1/TD_LP - 1/TBH));  % Vap. pres. heavy comp.
P_LP =  xD_LP*pL + (1-xD_LP)*pH;           % Stage pres.(Raoult's law)


% PI-controllers for Control of LP Column:
%-------------------------------------------------------------

D0_LP  = u0(6);                 % Nominal distillate flow (kmol/h)
B0_LP  = u0(7);                 % Nominal bottom flow (kmol/h)
VD0_LP = u0(8);                 % Nominal amount condensed (kmol/h)
LT0_LP = u0(9);                 % Nominal reflux (kmol/h)

% Top (condenser) holdup control
KcD  = -60;          
MDs  = 100;
D_LP = D0_LP + (MDs-MD_LP)*KcD; 

% Bottom (reboiler) holdup control
KcB  = -60;              
MBs  = 102.6;
B_LP = B0_LP + (MDs-MB_LP)*KcB; 

% Distillate composition control
KcXD   = 5000; TauI_xD = 0.001;
xDs_LP = 0.99;                
xD_LP_error  = (xDs_LP-xD_LP); Interr_xD = X(3*NT+3*Nt+3);
LT_LP = LT0_LP + KcXD*xD_LP_error + 1/TauI_xD*Interr_xD; 

% Top pressure control
KcP    = -0.05; TauI_P = -0.067;
Ps_LP  = 1.01e5;           
P_LP_error = (Ps_LP-P_LP); Interr_P = X(3*NT+3*Nt+4);
VD_LP = VD0_LP + KcP*P_LP_error + TauI_P*Interr_P;   

% Temperature control on tray 6 (using btm. flow from HP col.)
KcT = -3; TauI_T = -0.5;
Ts  = 380.83;      % K
T_error = (Ts - T_6); Interr_T = X(3*NT+3*Nt+5);
B_HP = B0_HP + KcT*T_error + TauI_T*Interr_T;

%---------------------------------------------------------------------
% Store all inputs and disturbances (goes into colamod1.m)
U1(1)=LT_HP;U1(2)=VB_HP;U1(3)=D_HP;U1(4)=B_HP;U1(5)=F_HP;
U1(6)=zF_HP;U1(7)=qF_HP;U1(8)=VD_HP;

% Store all inputs and disturbances (goes into colamod2.m)
U2(1)=LT_LP;U2(2)=VB_LP;U2(3)=D_LP;U2(4)=B_LP;U2(5)=B_HP;
U2(6)=zF_LP;U2(7)=qF_LP;U2(8)=VD_LP;

xprime1 = colamod1(t,X1,U1);
xprime2 = colamod2(t,X2,U2);
t
time = [time t];                     % Time vector

Ptop_HP = [Ptop_HP; P_HP];           % Top or btm. pressure response (HP)
Ptop_LP = [Ptop_LP; P_LP];           % Top  pressure response (LP)

Y1 =[D_HP B_HP VD_HP LT_HP VB_HP xD_HP xB_HP];
Y2 =[D_LP B_LP VD_LP LT_LP VB_LP xD_LP xB_LP];

Flow = [Flow; D_HP B_HP VD_HP LT_HP VB_HP D_LP B_LP VD_LP LT_LP VB_LP F_HP];

%************************** END LP COLUMN ****************************

IntegralError= [xD_HP_error;P_HP_error;xD_LP_error; P_LP_error; T_error];  % Integral error


xprime = [xprime1; xprime2;IntegralError]; 


