%**************************************************************
%    A PROGRAM FOR RUNNING THE HEAT INTEGRATED HP/LP COLUMNS
%
% - dynamic simulation of multieffect distillation columns
% - calls on multieffect.m 
% - inputs:  Xinitial 
% - outputs: Results_steady_state
%
% H.K.Engelien, 30/03/2001, edited by Magnus H. Nilsen 17/06/2002
%**************************************************************
close all
clear all
global  Y1 Y2 time P_HP P_LP Response 

% Number of stage in columns
NT = 36;                  % HP column
Nt = 48;                  % LP column
NS  = 5;                  % Sidestream stage LP col.
F_HP = 4320;
time = [ ];Response = [];te = [];xe = [];
delta_t=5;                %time between each online MPC-calc.(minutes)
proctime = delta_t/60;    %minutes

% Get nominal values, u0 and initial state vector from steady state case
load tuning27 Xss_I D_HP B_HP VD_HP LT_HP VB_HP D_LP B_LP VD_LP LT_LP VB_LP S xAD_HP xAD_LP P_top_LP xBB_LP 
start = [D_HP B_HP VD_HP LT_HP D_LP B_LP VD_LP LT_LP VB_LP VB_HP S];
Xss = Xss_I;

% Stepresponse model:Ystep=[y1/u1 .. yn/u1 y1/u2 .. yn/u2 .... y1/d1 .. yn/d1 ..]

% y1 = xAD_HP, y2 = xAD_LP, y3 = P_LP(condenser), y4 = T6_LP, y5 = xBB_LP
% u1 = LT_HP,  u2 = LT_LP,  u3 = VD_LP,           u4 = B_HP,  u5 = S_LP
% d1 = F_HP


load Ystep_5min Ystep_5min
Ystep = Ystep_5min;
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Controller - QDMC

yr    =[0.99;0.99;1.05e5;381.9;0.99];                        %Setpoints
ur    =[5468.8;3744.1;5187.5;2652.5;236.4];                  %Ideal values
u0    =[LT_HP;LT_LP;VD_LP;B_HP;S];                           %Input reference to deviation variable
d0    =F_HP;                                                 %Disturbance reference to deviation variable          
y0    =[xAD_HP;xAD_LP;P_top_LP;Xss(4*NT+3*Nt+6);xBB_LP];     %Measurement reference to deviation variable    
M     =180;  %180                                            %"Model horizon"(steps with delta_t)
N     =40;   %40                                             %Control horizon
P     =N+M-1;                                                %Prediction horizon
n1    =1;    %1                                              %Start constraint horizon
n2    =50;   %50                                             %End constraint horizon
gamma =[10000/0.001;1000/0.001;100/0.45e5;1/10;10/0.01];     %Dynamic opt. output weights
sigma =[1/5468.8;1/3744.1;1/5187.5;1/2652.5;10/236.4];       %Dynamic opt. input weights
lambda=[1/1093.8;1/748.8;1/1037.5;1/530.5;1/236.4];          %Dynamic opt. input change weights
q1    =[1/5468.8;1/3744.1;1/5187.5;1/2652.5;10/236.4];       %Steady-state opt. weights - ideal values
q2    =[1000/0.0001;10000/0.0001;1/0.45e5;1/10;10/0.01];     %Steady-state opt. weights - setpoints
q3a   =[1/0.01;1/0.01;1/0.45e5;1/10;1/0.01];                 %Steady-state opt. weights - upper output constraints
q3b   =[10000/0.001;10000/0.001;100000/00.0362e5;1/10;10/0.01];       %Steady-state opt. weights - lower output constraints
umin  =[0;0;0;0;0];                              %Min. input value
umax  =[10937.6;7488.2;10375;5305;472.8];        %Max. input value 
ymin  =[0.99;0.99;1.0138e5;361.9;0.989];         %Min. output value
ymax  =[1;1;2.5e5;391.9;1];                      %Max. output value
dumin =[-1093.8;-748.8;-1037.5;-530.5;-236.4];   %Min. change in input value
dumax =[ 1093.8; 748.8; 1037.5; 530.5; 236.4];   %Max. change in input value 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Simulation - closed loop

global Y    %measurements
linear=0;   %0:nonlinear, 1:linear
nntime=900; %sim. time(minutes)
ntime=nntime/delta_t;   %number of timesteps

%initial conditions
uk=u0;
duk=[0;0;0;0;0]; uk=uk+duk;
dk=d0;
dkold=d0;
ddk=0; dk=dk+ddk;
yk=y0;

%initialization
dUK=[];
T=0;
U=uk; D=d0; Yw=yk; 
Ur=ur; Umax=umax; Umin=umin;
Yr=yr; Ymax=ymax; Ymin=ymin;

%qdmc_init;
[YS,AS,ASD,H,E,E1,F,G,Gd,H1,Aeq1,A1,H2,Aeq2,A2,H3,Aeq3,A3,Iu,Iy,Ic,Is0,Is1,Is2]=qdmc_init(u0,d0,y0,Ystep,M,N,P,n1,n2,gamma,sigma,lambda,q1,q2,q3a,q3b);

ASp = AS; ASDp = ASD; dkold=dk; YSp = YS; dkp = dk; Up = U; Dp = D; Ywp = Yw; ukp = uk;
ykp = yk;

% Simulation/mpc - loop (START)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tcpu=cputime;

for i=1:ntime
   %Calculate input by qdmc    
   [duk,dUK,YS,YSD,YPD,umaxopt,uminopt,uropt,ymaxopt,yminopt,yropt,dkold]=qdmc(dUK,duk,dkold,uk,dk,yk,ur,yr,u0,d0,y0,M,N,P,n1,n2,YS,AS,ASD,H,E,E1,F,umin,umax,ymin,ymax,dumin,dumax,G,Gd,H1,Aeq1,A1,q1,H2,Aeq2,A2,q2,H3,Aeq3,A3,q3a,q3b,Iu,Iy,Ic,Is0,Is1,Is2);
   %Plot prediction

   %plotpred(dUK,YSD,YP,yropt,[],[],yk,uropt,[],[],uk)
   %plotpredakk(dUK,YSD,YPD,yk,yropt,yminopt,ymaxopt,uk,uropt,uminopt,umaxopt,ntime,i,delta_t)
   %pause
   %Update manipulated variables
   uk=uk+duk;
   %Data accumulation
   
   T=[T i-1];
   U=[U uk];, Yw=[Yw yk]; D=[D dk]; Up = [Up ukp]; Dp = [Dp dkp]; Ywp = [Ywp ykp]; 
   Ur=[Ur uropt]; Umax=[Umax umaxopt]; Umin=[Umin uminopt];
   Yr=[Yr yropt]; Ymax=[Ymax ymaxopt]; Ymin=[Ymin yminopt];
   start = [D_HP uk(4) VD_HP uk(1) D_LP B_LP uk(3) uk(2) VB_LP VB_HP uk(5) i delta_t dk];
   %pause
   %Implement input - process
   if linear==1
       ddkp=dkp-dkold;
       dkold=dkp;
       %duk (k-1)
       YSp=YS+AS*duk;
       ysp=YSp(Is0);
       YSp(Is1)=YSp(Is2);
       YSp=YSp+ASD*ddkp;
       pause
       ukp=uk+u0; dkp=dk+d0; ykp=yk+y0; %"absolute variables"
   else
      options = odeset('Maxstep',1/120,'InitialStep',1/1200);
      [t,x] = ode15s('multieffect',[0 proctime],Xss,options,start);
      t_final = size(x);
      %ud=[uk;dk];
      %options=odeset; options.Stats='on';
      %[t,x]=ode15s('oppgave_A',[0 1],xk,options,ud,p0); 
      xe =[xe;x];
      t = (i-1)*proctime + t;
      te = [te;t];
      X1final = [x(t_final(1),1:4*NT)]';            % Final states HP col
      X2final = [x(t_final(1),4*NT+1:4*NT+4*Nt)]';  % Final states LP col
      Xss = [X1final;X2final];
      %pause
      %yk=Y;
      %xk=x(size(x,1),:)';
  end
   D_HP=Y1(1);VD_HP=Y1(3);VB_HP=Y1(5);D_LP=Y2(1);B_LP=Y2(2);VB_LP=Y2(5);
   yk = [Y1(6) Y2(6) Y2(13) Y2(14) Y2(9)]'; 
   T=[T i];
   U=[U uk]; Yw=[Yw yk];  D=[D dk]; Up = [Up ukp]; Dp = [Dp dkp]; Ywp = [Ywp ykp]; 
   Ur=[Ur uropt]; Umax=[Umax umaxopt]; Umin=[Umin uminopt];
   Yr=[Yr yropt]; Ymax=[Ymax ymaxopt]; Ymin=[Ymin yminopt];
   ddk=zeros(1,1);
   %if i==6
   %    ddk=-21.6; dk=dk+ddk;
   %  yr(1)=0.985;yr(2)=0.985;
   %end
   %dk
   i
end
cputime-tcpu
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Simulation/mpc - loop (END)

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

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

% Plot results
T = T.*(delta_t/60);
plotresults(T,Yw,Yr,Ymin,Ymax,U,Ur,Umin,Umax)

Y1';                  % Vector for HP col : [D B VD LT VB xAD xBD xAB xBB F_HP]
Y2';                  % Vector for LP col : [D B VD LT VB xAD xBD xAB xBB xAS xBS S]             
X1final = [x(t_final(1),1:4*NT)]';            % Final states HP col
X2final = [x(t_final(1),4*NT+1:4*NT+4*Nt)]';  % Final states LP col

% Results
Xss_I = [x(t_final(1),:)]';
Xss = [X1final;X2final];
xC = 1 - x(t_final(1),4*NT+1:4*NT+Nt) - x(t_final(1),4*NT+Nt+1:4*NT+2*Nt);   % Final composition of xC in LP column

D_HP= Y1(1), B_HP= Y1(2), VD_HP= Y1(3), LT_HP= Y1(4)          
VB_HP= Y1(5), xAD_HP= Y1(6),xBD_HP= Y1(7), xAB_HP= Y1(8)
xBB_HP= Y1(9)

D_LP= Y2(1), B_LP= Y2(2), VD_LP= Y2(3), LT_LP= Y2(4), VB_LP= Y2(5)
xAD_LP= Y2(6), xBD_LP= Y2(7), xAB_LP= Y2(8), xBB_LP= Y2(9), xCB_LP =(1-xAB_LP-xBB_LP)
xAS_LP = Xss(NS), xBS_LP= Xss(Nt+NS), S=Y2(12)

%--------
D_tot = D_HP+D_LP
P_top_HP = P_HP(NT), Pbtm_LP = P_LP(1)
xAS_LP = Xss(4*NT+NS),xBS_LP= Xss(4*NT+Nt+NS)
xBB_LP= Xss(4*NT+Nt+1)
%------------
P_top_HP = P_HP(NT);P_top_LP = P_LP(Nt);Pbtm_LP = P_LP(1);

T_top_HP = x(t_final(1),4*NT); T_btm_LP = x(t_final(1),4*NT+3*Nt+1);

HvapA = 30384;HvapB = 36876;HvapC = 32396; % (kJ/kmol)
Hvap_top_HP = xAD_HP*HvapA + xBD_HP*HvapB +(1-xAD_HP-xBD_HP)*HvapC; % (kJ/kmol)
Hvap_btm_HP = xAB_HP*HvapA + xBB_HP*HvapB +(1-xAB_HP-xBB_HP)*HvapC; % (kJ/kmol) 
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 tuning30 Xss_I D_HP B_HP VD_HP LT_HP VB_HP xAD_HP xBD_HP xAB_HP xBB_HP D_LP B_LP VD_LP LT_LP VB_LP xAD_LP xBD_LP xAB_LP xBB_LP P_top_HP P_top_LP S 

%save steady_state_x0 Xss_I

%for i = 1:350
%    for i = 1:length(time)-1
%    if time(i) < time(i+1)
%    time(i) = time(i+1);
%    Response(i,:) = Response(i+1,:);
%end
%end
%end

% Plotting profiles and responses of the multi-effect columns
%**************************************************************
% HP column profiles
figure(1)
subplot(2,2,1); plot(x(t_final(1),1:NT),'b-')
title('Composition Profile, xA, in HP Column')
subplot(2,2,2); plot(x(t_final(1),2*NT+1:3*NT),'b-')      
title('Holdup Profile in HP Column')
subplot(2,2,3); plot(x(t_final(1),3*NT+1:4*NT),'b-')      
title('Temperature Profile in HP Column')
subplot(2,2,4); plot(P_HP/1e5,'r-')
title('Pressure Profile in HP Column')
hold;

% Bottom & top composition response & internal flow response
figure(2)
subplot(2,2,1); plot(te,xe(:,1),'r-')
title('Bottom Composition (xA) Response in HP Column')
subplot(2,2,2); plot(te,xe(:,NT),'r-')
title('Top Composition (xA) Response in HP Column')
subplot(2,2,3); plot(time,Response(:,5),'r-')
title('Boilup Response in HP Column')
subplot(2,2,4); plot(time,Response(:,4),'r-')
title('Reflux Response in HP Column')
hold

% Holdup response  and external flow response
figure(3)
subplot(2,2,1); plot(te,xe(:,2*NT+1),'r-')
title('Btm Holdup Response in HP Column')
subplot(2,2,2); plot(te,xe(:,3*NT),'r-')
title('Top Holdup Response in HP Column')
subplot(2,2,3); plot(time,Response(:,2),'r-')
title('Bottom flow Response in HP Column')
subplot(2,2,4); plot(time,Response(:,1),'r-')
title('Distillate Response in HP Column')
hold;

% LP column profiles
figure(4)
subplot(2,2,1); plot(x(t_final(1),4*NT+1:4*NT+Nt),'r-')
title('Composition Profile, xA, in LP Column')
subplot(2,2,2); plot(x(t_final(1),4*NT+2*Nt+1:4*NT+3*Nt),'r-')      
title('Holdup Profile in LP Column')
subplot(2,2,3); plot(x(t_final(1),4*NT+3*Nt+1:4*NT+4*Nt),'r-')      
title('Temperature Profile in LP Column')
subplot(2,2,4); plot(P_LP/1e5,'r-')
title('Pressure Profile in LP Column')
hold;

% Pressure and condensate flow response
figure(5)
subplot(2,2,1); plot(time,Response(:,12)/1e5,'r-')
title('Pressure Response HP Condenser'),ylabel('Pressure (bar)')
subplot(2,2,3); plot(time,Response(:,3),'r-')
title('Condensation (VD) in HP Column')
subplot(2,2,2); plot(time,Response(:,13)/1e5,'r-')
title('Pressure Response in LP Condenser'),ylabel('Pressure (bar)')
subplot(2,2,4); plot(time,Response(:,8),'r-')
title('Flow condensed (VD) in LP Column')
hold;

% Composition response (light component) & reflux and boilup response
figure(6)
subplot(2,2,1); plot(te,xe(:,4*NT+Nt+1),'r-')
title('Bottom Composition (xB) Response in LP Column')
subplot(2,2,2); plot(te,xe(:,4*NT+Nt),'r-')
title('Top Composition (xA) Response in LP Column')
subplot(2,2,3); plot(time,Response(:,10),'r-')
title('Boilup Response in LP Column')
subplot(2,2,4); plot(time,Response(:,9),'r-')
title('Reflux Response in LP Column')
hold;

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

% Temperature response
figure(8)
subplot(2,1,1);plot(te,xe(:,4*NT+3*Nt+6));
title('Temperature response in LP column, tray 6');
subplot(2,1,2); plot(time,Response(:,2));
title('Bottom flowrate, HP column')
hold;

% Ethanol composition in LP column
%figure(9)
%plot(xC(:,1:Nt),'r-')
%title('Compostion profile (xC) in LP Column')
%hold

% Composition response xA, xB in sidestream
figure(10)
subplot(2,2,1); plot(te,xe(:,4*NT+5))
title('Sidestream Composition (xA) Response in LP Column')
subplot(2,2,2); plot(te,xe(:,4*NT+Nt+5))
title('Sidestream Composition (xB) Response in LP Column')
subplot(2,2,3); plot(xC(:,1:Nt))
title('Compostion profile (xC) in LP Column')
subplot(2,2,4); plot(time,Response(:,9))
title('Reflux Response in LP Column')
hold;

%Sidestream response
figure(11)
subplot(2,1,1);plot(time, Response(:,14));
title('Sidestream Response in LP Column')
subplot(2,1,2); plot(te,xe(:,4*NT+Nt+1),'r-')
title('Bottom Composition (xB) Response in LP Column')
hold;

figure(12)
subplot(2,2,1);plot(te,xe(:,NT),'r-')
title('xAD-Response in HP Column')
subplot(2,2,2); plot(te,xe(:,4*NT+Nt),'r-')
title('xAD-Response in LP Column')
subplot(2,2,3); plot(te,xe(:,4*NT+Nt+1),'r-')
title('xBB-Response in LP Column')
subplot(2,2,4);plot(te,xe(:,4*NT+3*Nt+6));
title('T-Response in LP column, tray 6'), ylabel('T[K]');
hold;

figure(13)
subplot(2,1,1); plot(time,Response(:,12)/1e5,'r-')
title('Pressure Response HP Condenser'),ylabel('Pressure (bar)')
subplot(2,1,2); plot(time,Response(:,13)/1e5,'r-')
title('Pressure Response in LP Condenser'),ylabel('Pressure (bar)')
hold;