%**************************************************************
%        A PROGRAM FOR COMBINING THE HP/LP COLUMNS (MODIFIED)
%
%   - modified version od multieffect.m to be run using broyden1.m
%     or secant_cont.m
%   - 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
%**************************************************************

function f = multieffect_mod(ux,d,n,cs,m)

global  Y1 Y2 P_HP P_LP

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

% Inputs 
X = ux(1:3*NT+3*Nt);                 % States in the two columns
u = ux(3*NT+3*Nt+1:3*NT+3*Nt+8);     % Manipulated variables,u
d = d(1:3);

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

t=0;

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

% Feed (HP column)
F_HP  = d(1);              % Feedrate (kmol/h) 4320
zF_HP = d(2);              % Feed composition light component
qF_HP = d(3);                 % 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)

% Steady state variables for HP Column:
%-------------------------------------------------------------
 
D_HP  = u(1);         % Distillate flowrate HP column (kmol/h)
B_HP  = u(2);         % Bottom flowrate HP column (kmol/h)
LT_HP = u(3);         % Reflux HP column (kmol/h)
VB_HP = u(4);         % Boilup HP column (kmol/h)

MBs_HP = m(1);
MDs_HP = m(2);
xDs_HP = 0.99;

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

% Integrated condenser/reboiler

A = 1000;
U = 3000;                     % Overall h.t. coefficient (W/m^2 C)
DeltaT =TD_HP-X2(2*Nt+1);     % (Tcond - Treb) (K)
QC_HP=A*U*DeltaT*(3600/1000);  % Calculated area in exchanger (m^2)

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

VD_HP = QC_HP/Hvap_top_HP;        % Duty (kJ/h)
VB_LP = QC_HP/Hvap_btm_LP;        % Boilup LP column (kmol)


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

% 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;

xprime1 = colamod1(t,X1,U1);

Y1 =[D_HP B_HP VD_HP LT_HP VB_HP xD_HP xB_HP];

%************************** 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)


% Steady state variables for HP Column:
%-------------------------------------------------------------

D_LP  = u(5);              % Distillate flow LP column (kmol/h)
B_LP  = u(6);              % Bottoms flow LP column (kmol/h)
VD_LP = u(7);              % Condensation LP column (kmol/h)
LT_LP = u(8);              % Reflux LP column (kmol/h)

MBs_LP = m(3);
MDs_LP = m(4);
xDs_LP = 0.99;
Ps_LP  = 1.01e5;           % Top pressure set point

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

% 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)=F_LP;
U2(6)=zF_LP;U2(7)=qF_LP;U2(8)=VD_LP;

xprime2 = colamod2(t,X2,U2);

Y2 =[D_LP B_LP VD_LP LT_LP VB_LP xD_LP xB_LP];

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

% Controlled variables, c:
Hvap_btm_HP = 35440;              % From HYSYS at 5.2bar & 0.53 MeOH 
QB_HP = VB_HP*Hvap_btm_HP;        % Duty in HP column (kJ/h)
DeltaP_HP = P_HP(1)-P_HP(NT);     % Pressure drop in HP column
TB_LP = X2(2*Nt+1);               % Btm temperature LP column
T2_LP = X2(2*Nt+2);               % temerature tray 2, LP column
T4_LP = X2(2*Nt+4);               % temerature tray 4, LP column
T6_LP = X2(2*Nt+6);               % temerature tray 6, LP column
TB_HP = X1(2*NT+1);               % Btm. temp. in HP column

c= [QB_HP P_HP(NT) DeltaP_HP xB_LP xB_HP TB_LP T2_LP T4_LP T6_LP TB_HP B_HP...
    LT_LP LT_HP (QB_HP/F_HP) (QB_HP/LT_LP) (QB_HP/LT_HP) (B_HP/F_HP) ];
c = c(n);

Scale = [1000 1000 1 0.1 0.1 10 10 10 10 10 10 10 10 10 10 10 10];
S = Scale(n);

xprime = [xprime1; xprime2];

f = [xprime;(MD_HP-MDs_HP); (MB_HP-MBs_HP); (xD_HP-xDs_HP);...
     (MD_LP-MDs_LP);(MB_LP-MBs_LP); (P_LP(Nt)-Ps_LP); (xD_LP-xDs_LP);(c-cs)/S];
