%****************************************************************
%            A PROGRAM FOR COMBINING THE HP/LP COLUMNS
%
% - calles on colamod1.m and colamod2.m, the models for the 
%   HP and LP columns respectively
% - The no. of stages, feed data and nominal flows for control
%   loops must be changed for a new case
%
% H.K.Engelien, 30/03/2001, edited by Magnus H. Nilsen 17/06/2002
% edited H.K.Engelien 06/02/02
%****************************************************************

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

global  Y1 Y2 time Response 

NT = 36;                          % No. stages in HP column
Nt = 48;                          % No. stages in LP column
NS = 5;                           % Sidestream stage in LP col.
proctime = start(13)/60;
X1 = X(1:4*NT);                   % States in HP column
X2 = X(4*NT+1:4*NT+4*Nt);         % States in LP column
k = start(12);

%################## HP COLUMN #################################

% Feed (HP column)
F_HP   = start(14);  %4320                   % Feedrate (kmol/h)
zAF_HP = 0.73;  %0.73                 % Feed composition, A(methanol)
zBF_HP = 0.23;  %0.25                 % Feed composition, B(water)	
zCF_HP = 0.04;	%0.02			 	 % Feed composition, C(ethanol)
qF_HP  = 1;                      % Feed liq. fraction
TF_HP  = 406;					 % Temperature in feed (K)

% Get actual values from state vector
MB_HP  = X1(2*NT+1);             % Actual reboiler holdup (kmol)
MD_HP  = X1(3*NT);          	 % Actual condenser holdup (kmol)
xAD_HP = X1(NT);            	 % Concentration in distillate, A
xBD_HP = X1(2*NT);               % Concentration in distillate, B
xAB_HP = X1(1);                  % Concentration in btm product, A
xBB_HP = X1(NT+1);				 % Concentration in btm product, B
TD_HP  = X1(4*NT);          	 % Temperature in condenser (K)
TB_HP  = X1(3*NT+1);             % Temperature in reboiler (K)

% Physical data (taken @ 8 bar, from HYSYS)
TBA = 401.25;               % Boiling point, comp. A (K)
TBB = 443.57;                  % Boiling point, comp. B (K)
TBC = 414.92;               % Boiling point, comp. C (K)				 
HvapA = 30384;              % Hvap for component A (kJ/kmol)
HvapB = 36876;              % Hvap for component B (kJ/kmol)
HvapC = 32396;			    % Hvap for component C (kJ/kmol)
p0A = 8e5;                  % Vap. pressure for pure liq. comp. A (Pa)
p0B = 8e5;                  % Vap. pressure for pure liq. comp. B( Pa)
p0C = 8e5;                  % Vap. pressure for pure liq. comp. C (Pa)
R = 8.314;                  % Universal gas constant (kJ/kmol*K)

% Pure component constants for Antoine equation (Gmehling et al, 1977-1990)
A1=8.08097; B1=1582.271; C1=239.726;  % (T-range: 15C - 84C, Component 1)
A2=8.07131; B2=1730.630; C2=233.426;  % (T-range: 1C - 100C, Component 2)
A3=8.11220; B3=1592.864; C3=226.184;  % (T-range: 20C - 93C, Component 3)

% Partial molar volumes of pure components (Gmehling et al, 1977-1990)
V1=40.73;                    % Component 1 (cm^3/mol)
V2=18.07;                    % Component 2 (cm^3/mol)
V3=58.68;                    % Component 3 (cm^3/mol)

% Binary interaction parameters Aij for Wilson equation (Gmehling et al. 1977-1990)
A11=0;                       % Pure component 1 (cal/mol)
A22=0;                       % Pure component 2 (cal/mol)
A33=0;                       % Pure component 3 (cal/mol)
A12=-102.7135;%19.2547;                 % Binary mixture 1-2 (cal/mol) 
A21=695.4909;%554.0494;                % Binary mixture 1-2 (cal/mol) 
A13=-240.152;                % Binary mixture 1-3 (cal/mol) 
A31=637.320;                 % Binary mixture 1-3 (cal/mol) 
A23=478.0720;%904.717;                 % Binary mixture 2-3 (cal/mol) 
A32=875.1566;%447.171;                 % Binary mixture 2-3 (cal/mol) 

% End of data that have to be changed for different mixtures 
%*******************************************************************************

% Gas constant
R=1.98721;                  % cal/(mol*K)
i=1;
xA = xAD_HP;
xB = xBD_HP;

% Binary interaction parameters for Wilson equation 
L11=V1/V1*exp(-A11./(R*TD_HP));            % Pure component 1 
L22=V2/V2*exp(-A22./(R*TD_HP));            % Pure component 2 
L33=V3/V3*exp(-A33./(R*TD_HP));            % Pure component 3 
L12=V2/V1*exp(-A12./(R*TD_HP));            % Binary mixture 1-2 
L21=V1/V2*exp(-A21./(R*TD_HP));            % Binary mixture 1-2 
L13=V3/V1*exp(-A13./(R*TD_HP));            % Binary mixture 1-3  
L31=V1/V3*exp(-A31./(R*TD_HP));            % Binary mixture 1-3  
L23=V3/V2*exp(-A23./(R*TD_HP));            % Binary mixture 2-3  
L32=V2/V3*exp(-A32./(R*TD_HP));            % Binary mixture 2-3 

%Vapor pressures by Antoine equation
P1s=10.^(A1-B1./(TD_HP-273.15+C1))*(1e5/760);    % Component 1 (Pa)
P2s=10.^(A2-B2./(TD_HP-273.15+C2))*(1e5/760);   % Component 2 (Pa)
P3s=10.^(A3-B3./(TD_HP-273.15+C3))*(1e5/760);    % Component 3 (Pa)

% Activity coefficients by Wilson equation
G1(i) = exp(1-log(xA(i).*L11(i)+xB(i).*L12(i)+(1-xA(i)-xB(i)).*L13(i))...
          - xA(i).*L11(i)./(xA(i).*L11(i)+xB(i).*L12(i)+(1-xA(i)-xB(i)).*L13(i))...
          - xB(i).*L21(i)./(xA(i).*L21(i)+xB(i).*L22(i)+(1-xA(i)-xB(i)).*L23(i))...    
          - (1-xA(i)-xB(i)).*L31(i)./(xA(i).*L31(i)+xB(i).*L32(i)+(1-xA(i)-xB(i)).*L33(i)));    % Component 1

G2(i) = exp(1-log(xA(i).*L21(i)+xB(i).*L22(i)+(1-xA(i)-xB(i)).*L23(i))...
          - xA(i).*L12(i)./(xA(i).*L11(i)+xB(i).*L12(i)+(1-xA(i)-xB(i)).*L13(i))...
          - xB(i).*L22(i)./(xA(i).*L21(i)+xB(i).*L22(i)+(1-xA(i)-xB(i)).*L23(i))...    
          - (1-xA(i)-xB(i)).*L32(i)./(xA(i).*L31(i)+xB(i).*L32(i)+(1-xA(i)-xB(i)).*L33(i)));    % Component 2

G3(i) = exp(1-log(xA(i).*L31(i)+xB(i).*L32(i)+(1-xA(i)-xB(i)).*L33(i))...
          - xA(i).*L13(i)./(xA(i).*L11(i)+xB(i).*L12(i)+(1-xA(i)-xB(i)).*L13(i))...
          - xB(i).*L23(i)./(xA(i).*L21(i)+xB(i).*L22(i)+(1-xA(i)-xB(i)).*L23(i))...    
          - (1-xA(i)-xB(i)).*L33(i)./(xA(i).*L31(i)+xB(i).*L32(i)+(1-xA(i)-xB(i)).*L33(i)));    % Component 3

% Total pressure in condenser (top of column)

P_HP = (xA(i).*G1(i).*P1s(i) + xB(i).*G2(i).*P2s(i) + (1-xA(i)-xB(i)).*G3(i).*P3s(i));

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

% Nominal flows, u0, in HP column (kmol/h)
D0_HP  = start(1);
B0_HP  = start(2);    
VD0_HP = start(3);   
LT0_HP = start(4);
VB0_HP = start(10);

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

% Bottom (reboiler) holdup control
KcB = -60; MBs = 3230.7; 
%B_HP = B0_HP + (MBs-MB_HP)*KcB;
VB_HP = VB0_HP + (MBs-MB_HP)*KcB;

% Distillate composition control
%KcXAD = -2444; TauI_xAD = 1;   %KcXAD = -844444 TauI_xAD = 1
%xADs_HP = 0.99;                   
%xAD_HP_error  = (xADs_HP-xAD_HP); Interr_xAD = X(4*NT+4*Nt+1);   
LT_HP = LT0_HP ;%+ xAD_HP_error*KcXAD + KcXAD/TauI_xAD*Interr_xAD;  

%tl = length(time);
%if tl == 0
%    tl = 1;
%elseif time(tl) >= 1
%    LT_HP = LT0_HP - 20;
%end


% Boilup in HP column (constant boilup)
%Hvap_btm_HP = xAB_HP*HvapA + xBB_HP*HvapB +(1-xAB_HP-xBB_HP)*HvapC; % (kJ/kmol) 
%QB = 66.1;                                % Fixed boilup MW
%VB_HP = QB*3600*10^6/1000/Hvap_btm_HP;    % Boilup (kmol/h)

%#################### END HP COLUMN ###############################

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

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

Hvap_top_HP = xAD_HP*HvapA + xBD_HP*HvapB +(1-xAD_HP-xBD_HP)*HvapC; % (kJ/kmol)    
VD_HP = Q/Hvap_top_HP;        % Condensation in HP column (kmol/h)

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

%################### LP COLUMN #####################################

% Feed (LP column - from btm stage HP column)
zAF_LP = xAB_HP;                % Feed composition, A
zBF_LP = xBB_HP;                % Feed composition, B
zCF_LP = 1 - xAB_HP -xBB_HP;    % Feed composition, C
qF_LP  = 1;                     % Feed liq. fraction

% Get actual values
MB_LP  = X2(2*Nt+1);            % Actual reboiler holdup (LP)
MD_LP  = X2(3*Nt);              % Actual condenser holdup (LP)
xAD_LP = X2(Nt);                % Concentration in distillate
xBD_LP = X2(2*Nt);              % Concentration in distillate
xAB_LP = X2(1);                 % Concentration in btm product
xBB_LP = X2(Nt+1);              % Concentration in btm product
TD_LP  = X2(4*Nt);              % Temperature in condenser (LP)
T_6    = X2(3*Nt+6);		    % Temperature on stage 6(LP)
xAS_LP = X2(NS);                % Sidestream composition xA
xBS_LP = X2(Nt+NS);             % Sidestream composition xB

% Thermodynamic data (taken @ 1.013 bar, from HYSYS))
TBA   = 337.626;                % Boiling point, A 
TBB   = 373.141;                % Boiling point, B 
TBC	  = 351.315;	            % Boiling point, C
HvapA = 35291;                  % Hvap for A (kJ/kmol)
HvapB = 40651;                  % Hvap for B (kJ/kmol)
HvapC = 38489;                  % Hvap for C (kJ/kmol)
p0A   = 1.01e5;                 % Vap. pressure for A (Pa)
p0B   = 1.01e5;                 % Vap. pressure for B (Pa)
p0C	  = 1.01e5;			        % Vap. pressure for C (Pa)
R     = 8.314;                  % Universal gas constant (kJ/kmol*K)

% Estimating approx. feed temperature to LP col at 1 bar
TF_LP = xAB_HP^2*15.155 - xAB_HP*51.521 + 373.54;  % Feed temperature at 1 bar.
%TF_LP = xAB_HP*TBA + xBB_HP*TBB +(1-xAB_HP-xBB_HP)*TBC;	

% Pure component constants for Antoine equation (Gmehling et al, 1977-1990)
A1=8.08097; B1=1582.271; C1=239.726;  % (T-range: 15C - 84C, Component 1)
A2=8.07131; B2=1730.630; C2=233.426;  % (T-range: 1C - 100C, Component 2)
A3=8.11220; B3=1592.864; C3=226.184;  % (T-range: 20C - 93C, Component 3)

% Partial molar volumes of pure components (Gmehling et al, 1977-1990)
V1=40.73;     % Component 1 (cm^3/mol)
V2=18.07;     % Component 2 (cm^3/mol)
V3=58.68;     % Component 3 (cm^3/mol)

% Binary interaction parameters Aij for Wilson equation (Gmehling et al. 1977-1990)
A11=0;                       % Pure component 1 (cal/mol)
A22=0;                       % Pure component 2 (cal/mol)
A33=0;                       % Pure component 3 (cal/mol)
A12=19.2547;                 % Binary mixture 1-2 (cal/mol) 
A21=554.0494;                % Binary mixture 1-2 (cal/mol) 
A13=-240.152;                % Binary mixture 1-3 (cal/mol) 
A31=637.320;                 % Binary mixture 1-3 (cal/mol) 
A23=904.717;                 % Binary mixture 2-3 (cal/mol) 
A32=447.171;                 % Binary mixture 2-3 (cal/mol) 

% End of data that have to be changed for different mixtures 
%*******************************************************************************

% Gas constant
R=1.98721;     % cal/(mol*K)
i=1;
xA = xAD_LP;
xB = xBD_LP;

% Binary interaction parameters for Wilson equation 
L11=V1/V1*exp(-A11./(R*TD_LP));            % Pure component 1 
L22=V2/V2*exp(-A22./(R*TD_LP));            % Pure component 2 
L33=V3/V3*exp(-A33./(R*TD_LP));            % Pure component 3 
L12=V2/V1*exp(-A12./(R*TD_LP));            % Binary mixture 1-2 
L21=V1/V2*exp(-A21./(R*TD_LP));            % Binary mixture 1-2 
L13=V3/V1*exp(-A13./(R*TD_LP));            % Binary mixture 1-3  
L31=V1/V3*exp(-A31./(R*TD_LP));            % Binary mixture 1-3  
L23=V3/V2*exp(-A23./(R*TD_LP));            % Binary mixture 2-3  
L32=V2/V3*exp(-A32./(R*TD_LP));            % Binary mixture 2-3 

%Vapor pressures by Antoine equation
P1s=10.^(A1-B1./(TD_LP-273.15+C1))*(1e5/760);     % Component 1 (Pa)
P2s=10.^(A2-B2./(TD_LP-273.15+C2))*(1e5/760);     % Component 2 (Pa)
P3s=10.^(A3-B3./(TD_LP-273.15+C3))*(1e5/760);     % Component 3 (Pa)

% Activity coefficients by Wilson equation
G1(i) = exp(1-log(xA(i).*L11(i)+xB(i).*L12(i)+(1-xA(i)-xB(i)).*L13(i))...
          - xA(i).*L11(i)./(xA(i).*L11(i)+xB(i).*L12(i)+(1-xA(i)-xB(i)).*L13(i))...
          - xB(i).*L21(i)./(xA(i).*L21(i)+xB(i).*L22(i)+(1-xA(i)-xB(i)).*L23(i))...    
          - (1-xA(i)-xB(i)).*L31(i)./(xA(i).*L31(i)+xB(i).*L32(i)+(1-xA(i)-xB(i)).*L33(i)));    % Component 1

G2(i) = exp(1-log(xA(i).*L21(i)+xB(i).*L22(i)+(1-xA(i)-xB(i)).*L23(i))...
          - xA(i).*L12(i)./(xA(i).*L11(i)+xB(i).*L12(i)+(1-xA(i)-xB(i)).*L13(i))...
          - xB(i).*L22(i)./(xA(i).*L21(i)+xB(i).*L22(i)+(1-xA(i)-xB(i)).*L23(i))...    
          - (1-xA(i)-xB(i)).*L32(i)./(xA(i).*L31(i)+xB(i).*L32(i)+(1-xA(i)-xB(i)).*L33(i)));    % Component 2

G3(i) = exp(1-log(xA(i).*L31(i)+xB(i).*L32(i)+(1-xA(i)-xB(i)).*L33(i))...
          - xA(i).*L13(i)./(xA(i).*L11(i)+xB(i).*L12(i)+(1-xA(i)-xB(i)).*L13(i))...
          - xB(i).*L23(i)./(xA(i).*L21(i)+xB(i).*L22(i)+(1-xA(i)-xB(i)).*L23(i))...    
          - (1-xA(i)-xB(i)).*L33(i)./(xA(i).*L31(i)+xB(i).*L32(i)+(1-xA(i)-xB(i)).*L33(i)));    % Component 3

% Total pressure in condenser
P_LP=(xA(i).*G1(i).*P1s(i) + xB(i).*G2(i).*P2s(i)+(1-xA(i)-xB(i)).*G3(i).*P3s(i));

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

% Nominal flows, u0, in LP column (kmol/h)
D0_LP  = start(5);
B0_LP =  start(6);
VD0_LP = start(7);
LT0_LP = start(8);
S0_LP = start(11);

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

% Bottom (reboiler) holdup control
KcB  = -20; MBs  = 1056.3;
B_LP = B0_LP + (MBs-MB_LP)*KcB; 

% Distillate composition control
%KcXAD = 2461; TauI_xAD = 1;  %KcXAD=246153 TauI_xAD=1
%xADs_LP = 0.99;                         
%xAD_LP_error  = (xADs_LP-xAD_LP); Interr_xAD = X(4*NT+4*Nt+2);
LT_LP = LT0_LP ;%+ KcXAD*xAD_LP_error + KcXAD/TauI_xAD*Interr_xAD; 

%tl = length(time);
%if tl == 0
%    tl = 1;
%elseif time(tl) >= 1
%    LT_LP = LT0_LP - 30;
%end

% Top pressure control
%KcP = -0.024; TauI_P =0.2; %KcP=-0.024 TauI_P=0.2
%Ps_LP  = 1.013e5;           
%P_LP_error = (Ps_LP-P_LP);Interr_P = X(4*NT+4*Nt+3);
VD_LP = VD0_LP ;%+ KcP*P_LP_error ;%+ KcP/TauI_P*Interr_P;

% Bottoms composition control
%KcXBB = 5735; TauI_xBB = 0.53; %KcXBB=57353 TauI_xBB=0.53
%xBBs_LP = 0.99;                         
%xBB_LP_error  = (xBBs_LP-xBB_LP); Interr_xBB = X(4*NT+4*Nt+5);
S = S0_LP ;%+ KcXBB*xBB_LP_error ;%+ KcXBB/TauI_xBB*Interr_xBB; 

% Temperature control on tray 6
%KcT = -0.068;	TauI_T = 1.33; %KcT=-0.068 Tau_T=1.33
%Ts = 381.9;     % K Ts=381.9 from optimisation_case7.mat
%T_error = (Ts - T_6); Interr_T = X(4*NT+4*Nt+4);
B_HP = B0_HP ;%+ KcT*T_error ;%+ KcT/TauI_T*Interr_T;

% Boilup in LP column (given by combined reboiler/condenser)
Hvap_btm_LP = xAB_LP*HvapA + xBB_LP*HvapB +(1-xAB_LP-xBB_LP)*HvapC; % (kJ/kmol)			
VB_LP = Q/Hvap_btm_LP;        % Boilup in LP column (kmol/h)

%##################### END LP COLUMN #############################

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

% 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)=zAF_HP;U1(7)=zBF_HP;U1(8)=zCF_HP;U1(9)=qF_HP;U1(10)=VD_HP;%U1(11)=TF_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)=zAF_LP;U2(7)=zBF_LP;U2(8)=zCF_LP;U2(9)=qF_LP;U2(10)=VD_LP;U2(11)=TF_LP;U2(12)=S;

xprime1 = colamod1(t,X1,U1,start);
xprime2 = colamod2(t,X2,U2,start);

t;
tm = (k-1)*proctime + t;
time = [time;tm];                     % Time vector

Y1 =[D_HP B_HP VD_HP LT_HP VB_HP xAD_HP xBD_HP xAB_HP xBB_HP F_HP];
Y2 =[D_LP B_LP VD_LP LT_LP VB_LP xAD_LP xBD_LP xAB_LP xBB_LP xAS_LP xBS_LP S P_LP T_6];

Response = [Response; D_HP B_HP VD_HP LT_HP VB_HP D_LP B_LP VD_LP LT_LP VB_LP F_HP P_HP P_LP S];

%IntegralError= [xAD_HP_error;xAD_LP_error;P_LP_error;T_error;xBB_LP_error];  % Integral error
xprime = [xprime1; xprime2]; 


