function xprime=multivessel(t,X,U) 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%  Program written by Stathis Skouras, May 2001  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Multivessel - This is a nonlinear model of a multivessel batch distillation column with a decanter (top vessel),
%    NT1 theoretical stages on the top section, a middle vessel, NT2 theoretical stages on the bottom 
%    section and a reboiler. Reboiler is considered as an additional stage. 
%    The number of real stages is NT1 + NT2 + 1. For practical reasons we consider the decanter 
%    and the middle vessel as 2 additional stages. So total number of 'stages' are NT=NT1+NT2+3. 
%    Stages numbered from the top to the bottom. 
%    First 'stage' is condenser, then NT1 real stages in top section, then middle vessel 'stage', 
%    then NT2 real stages in bottom stage and finally the reboiler which is also a real stage.
%    Model assumptions: 
%    1. Staged distillation column sections
%    2. Perfect mixing and equilibrium in all stage and vessels
%    3. No vapor holdup on the stages
%    4. Constant liquid holdup on all stages (i.e., neglecting liquid dynamics)
%    5. Constant molar vapor flows in all stages (i.e., simplified energy balances 
%    6. Constant pressure
%    7. Total condenser
%    8. Ideal behaviour in the vapor phase
%    9. VLE modelled by NRTL equation 
%      
%
% Inputs and controller parameters 
%            t    - time in [min].
%            X    - Contains initial composition for all the stages.
%                   The first NT molar fractions are for light component A
%                   The next NT molar fractions are for intermediate component B
%                   Next NT states are molar liquid holdups in all stages
%                   Next NT states are temperature estimations in all stages
%            U(1) - Boilup, VB
%            U(2) - Nominal liquid flows in the column, L0
%            U(3) - Proportional controller gain at the top section, Kc1
%            U(4) - Integral time, tI
%            U(5) - Level controller set-point, Ms (kmol) 
%            U(6) - Proportional controller gain at the bottom section, Kc2 (kmol/(hr*C)
%            U(7) - Temperature setpoint in the bottom section (C), Tsp2
%            U(8) - Minimum liquid holdup in condenser and middle vessel, Hmin
%            
% Outputs:   xprime - vector with time derivative of all the states 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% The following data need to be changed for a new column.
% Number of stages in the top section: 
NT1=25;
% Number of stages in the bottom section (without the reboiler): 
NT2=25;
% Total number of 'stages'  
NT=NT1+NT2+3;                           % We include condenser, middle vessel and reboiler 
    
% Splitting the states
xA=X(1:NT)';                            % Liquid composition of A from top to btm
xB=X(NT+1:2*NT)';                       % Liquid composition of B from top to btm
M=X(2*NT+1:3*NT)';                      % Liquid hold up from top to btm
T=X(3*NT+1:4*NT);                       % Temperature from top to btm
IE1=X(4*NT+1);

% Inputs and disturbances
VB=U(1);                                % Boilup
L0=U(2);                                % Nominal liquid flows in the column (achieved only at steady state)
Kc1=U(3);                               % Controller gain at the top section (Kmol/(hr*C))
tI=U(4);                                % Integration time (hr)
Ms=U(5);                                % Level controller set-point (kmol) 
Kc2=U(6);                               % Controller gain at the bottom section (kmol/(hr*C))
Tsp2=U(7);                              % Temperature setpoint at the bottom section (C)
Hmin=U(8);                              % Minimum liquid holdup in condenser (kmol)

% THE MODEL

% Vapor-liquid equilibria by NRTL
Y=NRTL(xA,xB,T,NT);                     % Calls 'NRTL' VLE model and returns vapor compositions, temperature and stage pressure

% Splitting the states
yA=Y(1,1:NT);                           % Vapor composition for component A
yB=Y(2,1:NT);                           % Vapor composition for component B 
T=Y(3,1:NT);                            % Temperature estimation in each stage
DP=Y(4,1:NT);                           % Pressure difference vector (for BP calculations)

% Total liquid composition in the decanter 
z=[xA(1);xB(1);1-xA(1)-xB(1)];          

% Liquid-Liquid Equilibria
[Xaq,Xorg,Laq] = LLE(z)                 % Returns composition of aqueous and organic phase and split ratio

% Moles in the aqueous phase
Maq = M(1)*Laq

% Moles in the organic phase
Morg = M(1)*(1-Laq)

if Maq<0 
  pause
end
if Morg<0 
    pause
end


% Liquid flows in top section (L-controller applied. Drag all of the organic phase out of the decanter)
i=1:NT1+1; 
error=Morg-Ms;
L(i) = (Kc1*error + Kc1/tI*IE1)*(1-exp(-Morg/Hmin));       
dIE1dt = Morg-Ms;

% Location of temperature measurements
Tm2=T(NT1+2+(NT2+1)/2);    % Temperature measurement in the middle of the bottom section (stage 13 of bottom section)

% Liquid flows in column sections
% Bottom section
i=NT1+2:NT;  
L(i) = (L0 + Kc2*(Tm2-Tsp2-273.15))*(1-exp(-M(NT1+2)/Hmin));      % Precaution so as middle vessel does not run empty  

% Vapor flows assuming constant molar flows
i=1:NT;    V(i)=VB*ones(1,NT);


% Material balances for the column 1. Overall material balances 2. Component material balances

% Total condenser (decanter)
dMdt = V(2) - L(1);
dMxAdt= V(2).*yA(2) - L(1).*Xorg(1); 
dMxBdt= V(2).*yB(2) - L(1).*Xorg(2); 

% Top column section (stage right after the decanter)
i=2;
dMdt(i) = L(i-1)           - L(i)        + V(i+1)          - V(i);                   %=0
dMxAdt(i)= L(i-1).*Xorg(1) - L(i).*xA(i) + V(i+1).*yA(i+1) - V(i).*yA(i);
dMxBdt(i)= L(i-1).*Xorg(2) - L(i).*xB(i) + V(i+1).*yB(i+1) - V(i).*yB(i);

% Top column section (not the last stage)
i=3:NT1;
dMdt(i) = L(i-1)           - L(i)        + V(i+1)          - V(i);                   %=0
dMxAdt(i)= L(i-1).*xA(i-1) - L(i).*xA(i) + V(i+1).*yA(i+1) - V(i).*yA(i);
dMxBdt(i)= L(i-1).*xB(i-1) - L(i).*xB(i) + V(i+1).*yB(i+1) - V(i).*yB(i);

% Last stage of top column section (Notice the vapour bypass)
i=NT1+1;
dMdt(NT1+1) = L(NT1)           - L(NT1+1)            + V(NT1+3)            - V(NT1+1);                   %=0
dMxAdt(NT1+1)= L(NT1).*xA(NT1) - L(NT1+1).*xA(NT1+1) + V(NT1+3).*yA(NT1+3) - V(NT1+1).*yA(NT1+1);
dMxBdt(NT1+1)= L(NT1).*xB(NT1) - L(NT1+1).*xB(NT1+1) + V(NT1+3).*yB(NT1+3) - V(NT1+1).*yB(NT1+1);

% Middle vessel
dMdt(NT1+2)= L(NT1+1) - L(NT1+2);
dMxAdt(NT1+2)= L(NT1+1).*xA(NT1+1) - L(NT1+2).*xA(NT1+2);
dMxBdt(NT1+2)= L(NT1+1).*xB(NT1+1) - L(NT1+2).*xB(NT1+2);

% Bottom column section 
i=NT1+3:NT-1;
dMdt(i) = L(i-1)         - L(i)       + V(i+1)         - V(i);                   %=0
dMxAdt(i)= L(i-1).*xA(i-1) - L(i).*xA(i) + V(i+1).*yA(i+1) - V(i).*yA(i);
dMxBdt(i)= L(i-1).*xB(i-1) - L(i).*xB(i) + V(i+1).*yB(i+1) - V(i).*yB(i);

% Reboiler (assumed to be an equilibrium stage)
dMdt(NT) = L(NT-1)      - V(NT);
dMxAdt(NT)= L(NT-1)*xA(NT-1) - V(NT)*yA(NT);
dMxBdt(NT)= L(NT-1)*xB(NT-1) - V(NT)*yB(NT);

% Compute the derivative for the mole fractions from d(Mx) = x dM + M dx
i=1:NT;
dxAdt(i) = (dMxAdt(i) - xA(i).*dMdt(i))./M(i);
dxBdt(i) = (dMxBdt(i) - xB(i).*dMdt(i))./M(i);


% Output
xprime=[dxAdt';dxBdt';dMdt';DP';dIE1dt];