%**************************************************************
%    A PROGRAM FOR RUNNING THE HEAT INTEGRATED HP/LP COLUMNS
%
% - dynamic simulation of multieffect distillation columns
% - calls on multieffect.m 
% - inputs:  Xinitial - a vector with the states of the two 
%            columns (composition, holdup and temperature)
%            u0 - a vector with the nominal values for the 
%            contol loops in both columns
% - outputs: Results_steady_state
%
% H.K.Engelien, 30/03/2001
%**************************************************************

global Ptop_HP Ptop_LP Y1 Y2 time P_HP P_LP Flow Xss 

% Number of stage in columns
NT = 16;                  % HP column
Nt = 26;                  % LP column

Ptop_HP = [ ];Ptop_LP = [ ];time = [ ];Flow = [];

%load Results_steady_state_new  D_HP B_HP VD_HP LT_HP VB_HP D_LP B_LP VD_LP LT_LP Xss
%u0 = [D_HP B_HP VD_HP LT_HP VB_HP D_LP B_LP VD_LP LT_LP];
load Results_optimisation
u0 = [D_HPopt B_HPopt VD_HPopt LT_HPopt VB_HPopt D_LPopt B_LPopt VD_LPopt LT_LPopt];
Xss=Xfinal_opt(2:3*NT+3*Nt+1);
Xss=[Xss;0;0;0;0;0];


options = [];
[t,x] = ode15s('multieffect',[0 100000],Xss,options,u0);


t_final = size(x);

Y1';                  % Vector for HP col : [D B VD LT VB xD xB]
Y2';                  % Vector for LP col : [D B VD LT VB xD xB]             
X1final = [x(t_final(1),1:3*NT)]';            % Final states HP col
X2final = [x(t_final(1),3*NT+1:3*NT+3*Nt)]';  % Final states LP col

% Display results on screen
Xss = [X1final;X2final];
D_HP= Y1(1), B_HP= Y1(2), VD_HP= Y1(3), LT_HP= Y1(4)
VB_HP= Y1(5), xD_HP= Y1(6),xB_HP= Y1(7)
D_LP= Y2(1), B_LP= Y2(2), VD_LP= Y2(3), LT_LP= Y2(4)
VB_LP= Y2(5), xD_LP= Y2(6), xB_LP= Y2(7)
c=size(Ptop_HP);
P_top_HP = Ptop_HP(c(1))
P_top_LP = Ptop_LP(c(1))
T_top_HP = x(t_final(1),3*NT)
T_btm_LP = x(t_final(1),3*NT+2*Nt+1)
Hvap_top_HP = 31710;          % From HYSYS at 5.2bar & 0.99 MeOH
Hvap_btm_HP = 35440;          % From HYSYS at 5.2bar & 0.99 MeOH
Q_top_HP = Hvap_top_HP*VD_HP*1000/(3600*1000000)    % MW
Q_btm_HP = Hvap_btm_HP*VB_HP*1000/(3600*1000000)    % MW

%save Results_steady_state_new Xss D_HP B_HP VD_HP LT_HP VB_HP xD_HP xB_HP D_LP B_LP VD_LP LT_LP VB_LP xD_LP xB_LP P_top_HP P_top_LP


% Plotting profiles and responses of the multi-effect columns
%**************************************************************

% Plotting HP column profiles:
%------------------------------------

figure(1)
subplot(2,2,1); plot(x(t_final(1),1:NT),'r-')
title('Composition Profile in HP Column')
subplot(2,2,2); plot(x(t_final(1),NT+1:2*NT),'r-')      
title('Holdup Profile in HP Column')
subplot(2,2,3); plot(x(t_final(1),2*NT+1:3*NT),'r-')      
title('Temperature Profile in HP Column')
subplot(2,2,4); plot(P_HP,'r-')
title('Pressure Profile in HP Columns')
hold;

% Plotting HP column responses over the whole simulation time:
%-------------------------------------------------------------
% Pressure and condensate flow response
figure(2)
subplot(2,1,1); plot(time,Ptop_HP/1e5)
title('Pressure Response in HP Condenser')
ylabel('Pressure (bar)')
subplot(2,1,2); plot(time,Flow(:,5))
title('Boilup (VB) in HP Column')
hold;

% Bottom & top composition response & internal flow response
figure(3)
subplot(2,2,1); plot(t,x(:,1))
title('Bottom Composition (LK) Response in HP Column')
subplot(2,2,2); plot(t,x(:,NT))
title('Top Composition (LK) Response in HP Column')
subplot(2,2,3); plot(time,Flow(:,5))
title('Boilup Response in HP Column')
subplot(2,2,4); plot(time,Flow(:,4))
title('Reflux Response in HP Column')
hold

% Holdup response  and external flow response
figure(4)
subplot(2,2,1); plot(t,x(:,NT+1))
title('Btm Holdup Response in HP Column')
subplot(2,2,2); plot(t,x(:,2*NT))
title('Top Holdup Response in HP Column')
subplot(2,2,3); plot(time,Flow(:,11))
title('Feed flow  Response in HP Column')
subplot(2,2,4); plot(time,Flow(:,1))
title('Distillate Response in HP Column')
hold;

%********************************************
% Plotting LP column profiles:
%------------------------------------

figure(5)
subplot(2,2,1); plot(x(t_final(1),3*NT+1:3*NT+Nt),'b-')
title('Composition Profile in LP Column')
subplot(2,2,2); plot(x(t_final(1),3*NT+Nt+1:3*NT+2*Nt),'b-')      
title('Holdup Profile in LP Column')
subplot(2,2,3); plot(x(t_final(1),3*NT+2*Nt+1:3*NT+3*Nt),'b-')      
title('Temperature Profile in LP Column')
subplot(2,2,4); plot(P_LP,'b-')
title('Pressure Profile in LP Column')
hold;

% Plotting LP column responses over the whole simulation time:
%--------------------------------------------------------------

% Pressure and condensate flow response
figure(6)
subplot(2,1,1); plot(time,Ptop_LP)
title('Pressure Response in LP Condenser')
subplot(2,1,2); plot(time,Flow(:,8))
title('Flow condensed (VD) in LP Column')
hold;

% Coposition response (light component) & reflux and boilup response
figure(7)
subplot(2,2,1); plot(t,x(:,3*NT+1))
title('Bottom Composition (LK) Response in LP Column')
subplot(2,2,2); plot(t,x(:,3*NT+Nt))
title('Top Composition (LK) Response in LP Column')
subplot(2,2,3); plot(time,Flow(:,10))
title('Boilup Response in LP Column')
subplot(2,2,4); plot(time,Flow(:,9))
title('Reflux Response in LP Column')
hold;

% Holdup response and external flow response
figure(8)
subplot(2,2,1); plot(t,x(:,3*NT+Nt+1))
title('Btm Holdup Response in LP Column')
subplot(2,2,2); plot(t,x(:,3*NT+2*Nt))
title('Top Holdup Response in LP Column')
subplot(2,2,3); plot(time,Flow(:,7))
title('Bottom Flow Response in LP Column')
subplot(2,2,4); plot(time,Flow(:,6))
title('Distillate Response in LP Column')
hold;

% Temperature response
figure(9)
subplot(2,1,1);plot(t,x(:,3*NT+2*Nt+6));
title('Temperature response in LP column, tray 6');
subplot(2,1,2); plot(time,Flow(:,2));
title('Bottom flowrate, HP column')





