%% Model to simulate a steady state HEN
%  Optimal Design

% Topology to be investigated:

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%            1     2     3           %
%        ----0-----0-----0----       %
%   ----|                     |----  %
%        -------0------0------       %
%               4      5             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

close all;
clear all;
clc;

%% Parameters

par.w0  = 150;   %[kW/K] w= miCpi

par.wh1 = 50;   %[kW/K] 
par.wh2 = 30;   %[kW/K]
par.wh3 = 15;   %[kW/K]
par.wh4 = 70;
par.wh5 = 20;

par.Th1 = 205;  %[degC]
par.Th2 = 203;  %[degC]
par.Th3 = 220;  %[degC]
par.Th4 = 220;
par.Th5 = 248;

par.T0 = 130;   %[degC]

par.UA1 = 0.06;
par.UA2 = 0.31;
par.UA3 = 0.29;
par.UA4 = 2.29;
par.UA5 = 0.74;

par.DeltaTmin = 0.2; %[degC]
par.n = 0.65;
par.c0 = 4; %[$/m2] 
par.sc.x = [200*ones(11,1);100;100;1000*ones(5,1)];
% par.sc.j = 200;

%% Optimization


% x0 = [  Tend      T1     T2      T3     T4     T5        Th1out    Th2out    Th3out     Th4out     Th5out     w1    w2       Q1    Q2     Q3    Q4      Q5     UA1  UA2  UA3  UA4   UA5]

% x0   = [  182      150    164     174     173    186         165      176      182         145        170       53   124     1321   828   562   5190    1534]';


% 180/4
% x0   = [  175      153    168     178     172    185         165      177      182         144        173       80   80     1550   790   518   5102    1519]';

% 150/4
% x0   = [  180      139    161     176     169    181         183      178      188         154        176       33   117     314   737   472   4574    1422]';

% 150/2
x0   = [  187      151    167     179     178    192         168      176      182         149        180       75   75     1077   783  556   4911    1357]';


%% Optimization

x0 = x0./par.sc.x;
% x0 = 0.5*ones(size(x0))

A = []; b = []; Aeq = []; Beq = []; LB = 0.00*ones(18,1); UB = inf*ones(18,1);

options = optimset('display','iter','MaxFunEvals',9000,'TolCon',1e-10,'TolX',1e-10);
% options = optimset('display','iter','MaxFunEvals',9000,'TolCon',1e-6,'TolX',1e-6);
options = optimset('Algorithm','interior-point','display','iter','MaxFunEvals',9000,'TolCon',1e-12,'TolX',1e-12);
% options=[]

%% FIXED INLET TEMEPRATURE
[x,J,exitflag] = fmincon(@(x)Object_32_HXD_DJT(x,par),x0,A,b,Aeq,Beq,LB,UB,@(x)HEN_Constraints_32_HXD_DJT(x,par),options);
exitflag
x=x.*par.sc.x; % Unscale variables

Tend = x(1); T1 = x(2); T2 = x(3); T3 = x(4); T4 = x(5); T5 = x(6); Th1out = x(7); Th2out = x(8); Th3out = x(9); Th4out = x(10); Th5out = x(11);
w1 = x(12); w2 = x(13); 
Q1 = x(14); Q2 = x(15); Q3 = x(16); Q4 = x(17); Q5 = x(18);

T0 = par.T0; Th1 = par.Th1; Th2 = par.Th2; Th3 = par.Th3; Th4 = par.Th4; Th5 = par.Th5;
UA1 = par.UA1; UA2 = par.UA2; UA3 = par.UA3; UA4 = par.UA4; UA5 = par.UA5;

display(['   T1         T2         T3         T4         T5    Th1out    Th2out    Th3out    Th4out    Th5out [degC]'])
disp([T1 T2 T3 T4 T5   Th1out Th2out Th3out Th4out Th5out])
display(['  Tend [degC] = '])
disp(Tend)
w1_ratio = w1/par.w0;
w2_ratio = w2/par.w0;
display(['w1 ratio    w2 ratio      [J/K]'])
disp([w1_ratio w2_ratio])
display(['      w1        w2'])
disp([w1 w2])
display(['     UA1       UA2       UA3       UA4       UA5      [Wm2/K]'])
disp([UA1 UA2 UA3 UA4 UA5])
display(['DeltaTmin'])
display(['    Hot1      Cold1       Hot2    Cold2       Hot3    Cold3'])
Thot1 = Th1-T1;
Tcold1 = Th1out-T0;
Thot2 = Th2-T2;
Tcold2 = Th2out-T1;
Thot3 = Th3-T3;
Tcold3 = Th3out-T0;
disp([Thot1 Tcold1 Thot2 Tcold2 Thot3 Tcold3])
display(['    Hot4      Cold4       Hot5    Cold5'])
Thot4 = Th4-T4;
Tcold4 = Th4out-T0;
Thot5 = Th5-T5;
Tcold5 = Th5out-T4;
disp([Thot4 Tcold4 Thot5 Tcold5])
display(['  Q1            Q2        Q3        Q4       Q5'])
disp([Q1 Q2 Q3 Q4 Q5])


%% JAESCHKE TEMPERATURE

% x0 = [  Tend    T1        T2        T3     Th1out    Th2out      Th3out     w1       w2   Q1     Q2
x0   = [  182      150    164     174     173    186         165      176      182         145        170       53   124     1321   828   562   5190    1534]';
% x0   = [  207     169      204    216      163        176       155     60.45    69.54      2399   2120   5989]';

% x0 = 0.5*ones(size(x0));
x0 = x0./par.sc.x;

A = []; B = []; Aeq = []; Beq = []; LB = 0*ones(18,1); UB = inf*ones(18,1);
options = optimset('display','iter','MaxFunEvals',9000,'TolCon',1e-6,'TolX',1e-6);
% options = optimset('display','iter','MaxFunEvals',9000,'TolCon',1e-6,'TolX',1e-6);
options = optimset('Algorithm','interior-point','display','iter','MaxFunEvals',9000);

[xDJT,fval,exitflag] = fmincon(@(x)Object_32_HXD_DJT(x,par),x0,A,B,Aeq,Beq,LB,UB,@(x)HEN_Constraints_32_HXD_DJT(x,par),options);
exitflag
xDJT = xDJT.*par.sc.x;

TendDJT = xDJT(1); T1DJT = xDJT(2); T2DJT = xDJT(3); T3DJT = xDJT(4);
T4DJT = xDJT(5); T5DJT = xDJT(6);

Th1outDJT= xDJT(7); Th2outDJT = xDJT(8); Th3outDJT = xDJT(9); 
Th4outDJT = xDJT(10); Th5outDJT = xDJT(11);

w1DJT = xDJT(12); w2DJT = xDJT(13);
Q1DJT = xDJT(14); Q2DJT = xDJT(15); Q3DJT = xDJT(16);
Q4DJT = xDJT(17); Q5DJT = xDJT(18);

T0DJT = par.T0; Th1DJT = par.Th1; Th2DJT = par.Th2; Th3DJT = par.Th3;
Th4DJT = par.Th4; Th5DJT = par.Th5;
UA1DJT = par.UA1; UA2DJT = par.UA2; UA3DJT = par.UA3;
UA4DJT = par.UA4; UA5DJT = par.UA5;

display(['   T1_DJT         T2_DJT    T3_DJT      T4_DJT      T5_DJT      Th1out_DJT      Th2out_DJT   Th3out_DJT      Th4out_DJT      Th5out_DJT   [degC]'])
disp([T1DJT T2DJT T3DJT T4DJT T5DJT Th1outDJT Th2outDJT Th3outDJT Th4outDJT Th5outDJT])
display(['  Tend_DJT [degC] = '])
disp(TendDJT)
w1_ratio = w1/par.w0;
w2_ratio = w2/par.w0;
% display(['w1 ratio    w2 ratio      [J/K]'])
% disp([w1_ratio w2_ratio])
display(['      w1_DJT        w2_DJT'])
disp([w1DJT w2DJT])
display(['     UA1_DJT       UA2_DJT    UA3_DJT    UA4_DJT    UA5_DJT    [Wm2/K]'])

disp([UA1DJT UA2DJT UA3DJT UA4DJT UA5DJT])
display(['DeltaTmin_DJT'])
display(['    Hot1      Cold1       Hot2    Cold2     Hot3    Cold3'])
Thot1DJT = Th1DJT-T1DJT;
Tcold1DJT = Th1outDJT-T0DJT;
Thot2DJT = Th2DJT-T2DJT;
Tcold2DJT = Th2outDJT-T0DJT;
Thot3DJT = Th3DJT-T3DJT;
Tcold3DJT = Th3outDJT-T0;
disp([Thot1DJT Tcold1DJT Thot2DJT Tcold2DJT Thot3DJT Tcold3DJT])
display(['  Q1_DJT       Q2_DJT      Q3_DJT'])
disp([Q1DJT Q2DJT Q3DJT])



% %% VARIATION AND TRENDS
% % Defining the initial points
% c0 = 1;
% co_end = 5;
% DeltaC0 = 0.1;
% n = (co_end-c0)/DeltaC0;
% c0_vec = [];
% 
% % T0 = 120;
% % T0_end = 140;
% % DeltaT0 = 1;
% % n = (T0_end-T0)/DeltaT0; % # of points
% % T0_vec = [];
% 
% % w0 = 90;
% % w0_end = 150;
% % Deltaw0 = 1;
% % n = (w0_end-w0)/Deltaw0;
% % w0_vec = [];
% 
% % Utilizing the resulting variables
% T1_vec = [];
% T2_vec = [];
% T3_vec = [];
% Tend_vec = [];
% w1_vec = [];
% w2_vec = [];
% UA1_vec = [];
% UA2_vec = [];
% UA3_vec = [];
% Th1out_vec = [];
% Th2out_vec = [];
% Th3out_vec = [];
% Q1_vec = [];
% Q2_vec = [];
% Q3_vec = [];
% exitflag_vec = [];
% DeltaTHX1_vec = [];
% DeltaTHX2_vec = [];
% DeltaTHX3_vec = [];
% 
% for i=1:n;
%     [x,J,exitflag] = fmincon(@(x)Object_21_HXD(x,par,c0),x0,A,b,Aeq,Beq,LB,UB,@(x)HEN_Constraints_21_HXD(x,par),options);
% %     x0 = x;
%     x=x.*par.sc.x; % Unscale variables
% exitflag
% 
% Tend = x(1); T1 = x(2); T2 = x(3); T3 = x(4); Th1out = x(5); Th2out = x(6); Th3out = x(7); w1 = x(8);
% w2 = x(9); Q1 = x(10); Q2 = x(11); Q3 = x(12); UA1 = x(13); UA2 = x(14); UA3 = x(15);
% 
% 
% w0 = par.w0; 
% wh1 = par.wh1; wh2 = par.wh2; wh3 = par.wh3;
% Th1 = par.Th1; Th2 = par.Th2; Th3 = par.Th3;
% T0 = par.T0;
% 
% Thot1 = Th1-T1;
% Tcold1 = Th1out-T0;
% DeltaTHX1 = Thot1/Tcold1;
% 
% Thot2 = Th2-T2;
% Tcold2 = Th2out-T1;
% DeltaTHX2 = Thot2/Tcold2;
% 
% Thot3 = Th3-T3;
% Tcold3 = Th3out-T0;
% DeltaTHX3 = Thot3/Tcold3;
% 
% T1_vec(i)= T1;
% T2_vec(i) = T2;
% T3_vec(i) = T3;
% Tend_vec(i) = Tend;
% w1_vec(i) = w1;
% w2_vec(i) = w2;
% UA1_vec(i) = UA1;
% UA2_vec(i) = UA2;
% UA3_vec(i) = UA3;
% exitflag_vec(i) = exitflag;
% Th1out_vec(i) = Th1out;
% Th2out_vec(i) = Th2out;
% Th3out_vec(i) = Th3out;
% Q1_vec(i) = Q1;
% Q2_vec(i) = Q2;
% Q3_vec(i) = Q3;
% DeltaTHX1_vec(i) = DeltaTHX1;
% DeltaTHX2_vec(i) = DeltaTHX2;
% DeltaTHX3_vec(i) = DeltaTHX3;
% 
% c0_vec(i) = c0;
% % T0_vec(i) = T0;
% % w0_vec(i) = w0;
% 
% c0 = c0+DeltaC0;
% % T0 = T0+DeltaT0;
% % w0 = w0+Deltaw0;
% 
% end
% 
% % disp(Tend_vec(1))
% % disp(T1_vec(1))
% % disp(T2_vec(1))
% % disp(Th1out_vec(1))
% % disp(Th2out_vec(1))
% 
% disp(UA1_vec);
% disp(UA2_vec);
% 
% figure(1)
% plot(c0_vec,Tend_vec,'--')
% xlabel('c0 cost factor')
% ylabel('T_{end}')
% legend('T_{end}')
% 
% figure(2)
% plot(c0_vec,exitflag_vec,'--')
% xlabel('c0')
% ylabel('exitflag')
% 
% figure(3)
% plot(c0_vec,w1_vec,'x')
% hold on
% plot(c0_vec,w2_vec,'o')
% xlabel('c0 cost factor')
% ylabel('w_{i} = mC_{p}')
% legend('w_{1}','w_{2}')
% 
% 
% figure
% subplot(3,3,1)
% plot(c0_vec,T1_vec,'--')
% % axis([120 140 210 230])
% xlabel('c0 cost factor')
% ylabel('T1 [^{\circ}C]')
% legend('T1')
% 
% subplot(3,3,2)
% plot(c0_vec,T2_vec,'--')
% xlabel('c0 cost factor')
% ylabel('T2 [^{\circ}C]')
% legend('T2')
% 
% subplot(3,3,3)
% plot(c0_vec,T3_vec,'--')
% xlabel('c0 cost factor')
% ylabel('T3 [^{\circ}C]')
% legend('T3')
% 
% subplot(3,3,4)
% plot(c0_vec,Th1out_vec,'--')
% xlabel('c0 cost factor')
% ylabel('Th1out [^{\circ}C]')
% legend('Th1out')
% 
% subplot(3,3,5)
% plot(c0_vec,Th2out_vec,'--')
% % axis([120 140 190 210])
% xlabel('c0 cost factor')
% ylabel('Th2out [^{\circ}C]')
% legend('Th2out')
% 
% subplot(3,3,6)
% plot(c0_vec,Th3out_vec,'--')
% % axis([120 140 220 250])
% xlabel('c0 cost factor')
% ylabel('Th3out [^{\circ}C]')
% legend('Th3out')
% 
% subplot(3,3,7)
% plot(c0_vec,DeltaTHX1_vec,'--')
% xlabel('c0 cost factor')
% ylabel('DeltaT HX1')
% legend('DeltaT HX1')
% 
% subplot(3,3,8)
% plot(c0_vec,DeltaTHX2_vec,'--')
% xlabel('c0 cost factor')
% ylabel('DeltaT HX2')
% legend('DeltaT HX2')
% 
% subplot(3,3,9)
% plot(c0_vec,DeltaTHX3_vec,'--')
% xlabel('c0 cost factor')
% ylabel('DeltaT HX3')
% legend('DeltaT HX3')


%% STOP %%




% 
% % subplot(3,3,7)
% % plot(c0_vec,Th1out_vec,'--')
% % xlabel('c0 (Pua/Pt)')
% % ylabel('Th1_out')
% % legend('Th1_out')
% % 
% % subplot(4,2,8)
% % plot(c0_vec,Th2out_vec,'--')
% % xlabel('c0 (Pua/Pt)')
% % ylabel('Th2_out')
% % legend('Th2_out')
% % 
% % figure(2)
% % plot(c0_vec,Q1_vec,'x')
% % hold on
% % plot(c0_vec,Q2_vec,'o')
% % xlabel('c0 (Pua/Pt)')
% % ylabel('Q')
% % legend('Q1', 'Q2')
% % 
% % figure(3)
% % plot(c0_vec,exitflag_vec,'x')
% % axis([1 5 -5 5])
% % xlabel('c0 (Pua/Pt')
% % ylabel('exitflag')
% % legend('exitflag')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% figure
% subplot(4,2,1)
% plot(T0_vec,Tend_vec,'--')
% % axis([120 140 210 230])
% xlabel('w0 [J/K]')
% ylabel('Tend [^{\circ}C]')
% legend('Tend')
% 
% subplot(4,2,2)
% plot(w0_vec,w1_vec,'o')
% hold on
% plot(w0_vec,w2_vec,'x')
% xlabel('w0 [J/K]')
% ylabel('mC_p parallell streams [J(K)]')
% legend('w_1 [J/K]','w_2 [J/K]')
% 
% subplot(4,2,3)
% plot(w0_vec,UA1_vec,'o')
% xlabel('w0 [J/K]')
% ylabel('UA1')
% legend('UA1')
% 
% subplot(4,2,4)
% plot(w0_vec,UA2_vec,'x')
% xlabel('w0 [J/K]')
% ylabel('UA2')
% legend('UA2')
% 
% subplot(4,2,5)
% plot(w0_vec,T1_vec,'x')
% % axis([120 140 190 210])
% xlabel('w0 [J/K]')
% ylabel('T1')
% legend('T1')
% 
% subplot(4,2,6)
% plot(w0_vec,T2_vec,'x')
% % axis([120 140 220 250])
% xlabel('w0 [J/K]')
% ylabel('T2')
% legend('T2')
% 
% subplot(4,2,7)
% plot(w0_vec,Th1out_vec,'x')
% xlabel('w0 [J/K]')
% ylabel('Th1_out')
% legend('Th1_out')
% 
% subplot(4,2,8)
% plot(w0_vec,Th2out_vec,'o')
% xlabel('w0 [J/K]')
% ylabel('Th2_out')
% legend('Th2_out')
% 
% figure(2)
% plot(w0_vec,Q1_vec,'x')
% hold on
% plot(w0_vec,Q2_vec,'o')
% xlabel('w0 [J/K]')
% ylabel('Q')
% legend('Q1', 'Q2')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% figure(1)
% plot(T0_vec,Tend_vec)
% axis([120 140 215 235])
% xlabel('Inlet temperature T0 [^{\circ}C]')
% ylabel('Outlet temperature, Tend [^{\circ}C]')
% legend('Tend')
% 
% figure(2)
% plot(T0_vec,w1_vec,'o')
% hold on
% plot(T0_vec,w2_vec,'--')
% axis([120 140 0 100])
% xlabel('Inlet temperature T0 [^{\circ}C]')
% ylabel('Heat capasities parallell streams (split) [J/K])')
% legend('w1 [J/K]','w2 [J/K]')
% 
% figure(3)
% plot(T0_vec,T1_vec,'o')
% hold on
% plot(w0_vec,T2_vec,'--')
% plot(w0_vec,Tend_vec)
% % axis([120 140 190 255])
% xlabel('Inlet temperature T0 [^{\circ}C]')
% ylabel('Temperature [^{\circ}C]')
% legend('T1','T2','Tend')
% 
% figure(4)
% plot(w0_vec,Thot1_vec,'o')
% hold on
% plot(w0_vec,Tcold1_vec,'--')
% plot(w0_vec,Thot2_vec,'.')
% plot(w0_vec,Tcold2_vec,'-')
% % axis([120 140 0 50])
% xlabel('Inlet temperature T0 [^{\circ}C]')
% ylabel('DeltaTmin')
% legend('Hot side HX1','Cold side HX1','Hot side HX2','Cold side HX3')
% 
% figure(5)
% plot(w0_vec,UA1_vec,'-')
% xlabel('Heat capacity inlet cold stream (w0) [J/K]')
% ylabel('Heat Exchanger Duty UA [J/K]')
% legend('UA1')
% 
% figure(6)
% plot(w0_vec,UA2_vec,'-')
% axis([120 140 1300 1650])
% xlabel('Heat capacity inlet cold stream (w0) [J/K]')
% ylabel('Heat Exchanger Duty UA [J/K]')
% legend('UA2')




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
% % %JAESCHKE TEMPERATURES
% % x0 = x0./par.sc.x;
% % A = []; B = []; Aeq = []; Beq = []; LB = 0*ones(11,1); UB = inf*ones(11,1);
% % options = optimset('display','iter','MaxFunEvals',2000);
% % 
% % [xDJT,J,exitflag] = fmincon(@(x)Object_11_HXD(x,par),x0,A,B,Aeq,Beq,LB,UB,@(x)HEN_Constraints_11_HXD_DJT(x,par),options);
% % exitflag
% % xDJT = xDJT.*par.sc.x; % Unscale variables
% % 
% % TendDJT = xDJT(1); T1DJT = xDJT(2); T2DJT = xDJT(3); Th1outDJT = xDJT(4);
% % Th2outDJT = xDJT(5); w1DJT = xDJT(6); w2DJT = xDJT(7); Q1DJT = xDJT(8); 
% % Q2DJT = xDJT(9); UA1DJT = xDJT(10); UA2DJT = xDJT(11); 
% % Th1DJT = par.Th1; Th2DJT = par.Th2; T0 = par.T0;
% % 
% % display(['   T1(DJT)    T2(DJT)      [degC]'])
% % disp([T1DJT T2DJT])
% % display(['  TendDJT [degC] = '])
% % disp(TendDJT)
% % w1DJT_ratio = w1DJT/par.w0;
% % w2DJT_ratio = w2DJT/par.w0;
% % display(['w1DJT ratio   w2DJT ratio      [J/K]'])
% % disp([w1DJT_ratio w2DJT_ratio])
% % display(['     UA1DJT       UA2DJT     [Wm2/K]'])
% % disp([UA1DJT UA2DJT])
% % display(['DeltaTmin DJT'])
% % display(['    Hot1      Cold1       Hot2    Cold2'])
% % Thot1DJT = Th1DJT-T1DJT;
% % Tcold1DJT = Th1outDJT-T0;
% % Thot2DJT = Th2DJT-T2DJT;
% % Tcold2DJT = Th2outDJT-T0;
% % disp([Thot1DJT Tcold1DJT Thot2DJT Tcold2DJT])
% display(['Required heating of cold streams (test):'])
% disp([Q1_test_cold Q2_test_cold]')
% display(['Required heat loss for hot streams (test):'])
% disp([Q1_test_hot Q2_test_hot]')


