%% Model to simulate a steady state 4:1 HEN
% Topology to be investigated:

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


close all;
clear all;
clc;

%% Parameters

% Heat Capacity rates
par.w0  = 100;   %[kW/degC] w= miCpi
par.wh1 = 50;    %[kW/degC] 
par.wh2 = 30;    %[kW/degC]
par.wh3 = 15;    %[kW/degC]
par.wh4 = 25;    %[kW/degC]
par.wh5 = 70;    %[kW/degC]

% Hot streams inlet temperature
par.Th1 = 190;  %[degC]
par.Th2 = 203;  %[degC]
par.Th3 = 220;  %[degC]
par.Th4 = 235;  %[degC]
par.Th5 = 210;  %[degC]

% Cold stream inlet temperature
par.T0 = 130;   %[degC]

% UA values for each heat exchanger
par.UA1 = 5;    %[kWm2/degC]
par.UA2 = 7;    %[kWm2/degC]
par.UA3 = 10;   %[kWm2/degC]
par.UA4 = 12;   %[kWm2/degC]
par.UA5 = 9;    %[kWm2/degC]

% Operating prices for each heat exchanger
par.P1 = 1;     %[$/kW]
par.P2 = 1.2;   %[$/kW]
par.P3 = 1.3;   %[$/kW]
par.P4 = 1.5;   %[$/kW]
par.P5 = 1.4;   %[$/kW]

%Inequality constraint
par.DeltaTmin = 0.5; %[degC]

% Scaling vector
par.sc.x = [200*ones(11,1);100;100;1000*ones(5,1)];
par.sc.j = 200;

% Defining parameters
Th1 = par.Th1; Th2 = par.Th2; Th3 = par.Th3; Th4 = par.Th4; Th5 = par.Th5;
T0 = par.T0;

%% OPTIMAL OPERATION  

% Guessing outlet variables
% x0 = [Tend T1 T2 T3 T4 T5 Th1out Th2out Th3out Th4out Th5out w1 w2 ...
%      [Q1 Q2 Q3 Q4 Q5]

x0   = [138 131 133 138 138 140 188 198 200 215 190 60 40 ...
        59 137 297 333 200]';
% x0   = [207 160 176 187 204 215 146 166 169 185 174 71 29 ...
%           1.9224e+03 778.4439 581.1345 921.1994 3.3767e+03]';



% Scaling variables
% x0 = x0./par.sc.x;

% Minimizing cost function based on equality constraints 
% using fmincon
A = []; b = []; Aeq = []; Beq = [];
LB = 0*ones(23,1); UB = inf*ones(23,1);

options = optimset('Algorithm','interior-point','display','iter',...
    'MaxFunEvals',9000,'TolCon',1e-12,'TolX',1e-12);

options = optimset('Algorithm','active-set','display','iter',...
    'MaxFunEvals',9000,'TolCon',1e-11,'TolX',1e-11);

options = optimset('display','iter',...
    'MaxFunEvals',9000,'TolCon',1e-10,'TolX',1e-10);

[x,J,exitflag] = fmincon(@(x)Object_41(x,par),x0,A,b,Aeq,Beq,...
    LB,UB,@(x)HEN_Constraints_41(x,par),options);
exitflag

% Unscaling variables
% x = x.*par.sc.x; 


% RESULTS
% Outlet temperatures
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);
% Split
w1 = x(12); w2 = x(13);
% Heat transfer
Q1 = x(14); Q2 = x(15); Q3 = x(16); Q4 = x(17); Q5 = x(18);
% Split ratio
w1_rat = w1/par.w0;
w2_rat = w2/par.w0;
% Delta Ts
DeltaT_hot1 = Th1 - T1;
DeltaT_hot2 = Th2 - T2;
DeltaT_hot3 = Th3 - T3;
DeltaT_hot4 = Th4 - T4;
DeltaT_hot5 = Th5 - T5;
DeltaT_cold1 = Th1out - T0;
DeltaT_cold2 = Th2out - T1;
DeltaT_cold3 = Th3out - T2;
DeltaT_cold4 = Th4out - T3;
DeltaT_cold5 = Th5out - T0;

% Displaying the results
display(['  Tend [degC] = '])
disp(Tend)
display(['   T1          T2         T3        T4        T5   [degC]'])
disp([T1 T2 T3 T4 T5])
display(['    Th1out    Th2out    Th3out    Th4out    Th5out   [degC]'])
disp([Th1out Th2out Th3out Th4out Th5out])
display(['      w1        w2'])
disp([w1 w2])
display(['  w1 ratio   w2 ratio [%]'])
disp([w1_rat w2_rat])
display([' DeltaT hot side '])
display(['    HX1        HX2       HX3       HX4       HX5   '])
disp([DeltaT_hot1 DeltaT_hot2 DeltaT_hot3 DeltaT_hot4 DeltaT_hot5])
display([' DeltaT cold side '])
display(['    HX1        HX2       HX3       HX4       HX5   '])
disp([DeltaT_cold1 DeltaT_cold2 DeltaT_cold3 DeltaT_cold4 DeltaT_cold5])


%% OPERATION USING THE JAESCHKE TEMPERATURE

% Guessing outlet variables
% x0 = [Tend T1 T2 T3 T4 T5 Th1out Th2out Th3out Th4out Th5out w1 w2...
%      [Q1 Q2 Q3 Q4 Q5]
x0   = [138 131 133 138 138 140 188 198 200 215 190 60 40 ...
        59 137 297 333 200]';
    
% Scaling variables    
% x0 = x0./par.sc.x;

% Defining parameters
Th1 = par.Th1; Th2 = par.Th2; Th3 = par.Th3; Th4 = par.Th4; Th5 = par.Th5;
T0 = par.T0;


% Minimizing cost function based on equality constraints and Jaeschke temp 
% using fmincon
A = []; b = []; Aeq = []; Beq = [];
LB = 0*ones(23,1); UB = inf*ones(23,1);

options = optimset('Algorithm','interior-point','display','iter',...
    'MaxFunEvals',9000,'TolCon',1e-12,'TolX',1e-12);

options = optimset('Algorithm','active-set','display','iter',...
    'MaxFunEvals',9000,'TolCon',1e-11,'TolX',1e-11);

options = optimset('display','iter',...
    'MaxFunEvals',9000,'TolCon',1e-10,'TolX',1e-10);

[xDJT,J,exitflag] = fmincon(@(x)Object_41(x,par),x0,A,b,Aeq,Beq,...
    LB,UB,@(x)HEN_Constraints_41_DJT(x,par),options);
exitflag

%Unscaling variables
% xDJT = xDJT.*par.sc.x;


% RESULTS
% Outlet temperatures
Tend_DJT = xDJT(1);
T1_DJT = xDJT(2); T2_DJT = xDJT(3); T3_DJT = xDJT(4); T4_DJT = xDJT(5);
T5_DJT = xDJT(6);
Th1out_DJT = xDJT(7); Th2out_DJT = xDJT(8); Th3out_DJT = xDJT(9);
Th4out_DJT = xDJT(10); Th5out_DJT = xDJT(11);
% Split
w1_DJT = xDJT(12); w2_DJT = xDJT(13);
% Heat transfer
Q1_DJT = xDJT(14); Q2_DJT = xDJT(15); Q3_DJT = xDJT(16); Q4_DJT = xDJT(17);
% Split ratio
w1_rat_DJT = w1_DJT/par.w0;
w2_rat_DJT = w2_DJT/par.w0;
% Delta Ts
DeltaT_hot1_DJT = Th1 - T1_DJT;
DeltaT_hot2_DJT = Th2 - T2_DJT;
DeltaT_hot3_DJT = Th3 - T3_DJT;
DeltaT_hot4_DJT = Th4 - T4_DJT;
DeltaT_hot5_DJT = Th5 - T5_DJT;
DeltaT_cold1_DJT = Th1out_DJT - T0;
DeltaT_cold2_DJT = Th2out_DJT - T1_DJT;
DeltaT_cold3_DJT = Th3out_DJT - T2_DJT;
DeltaT_cold4_DJT = Th4out_DJT - T3_DJT;
DeltaT_cold5_DJT = Th5out_DJT - T0;

% Displaying the results
display(['  Tend DJT [degC] = '])
disp(Tend_DJT)
display(['   T1 DJT   T2 DJT     T3 DJT    T4 DJT    T5 DJT  [degC]'])
disp([T1_DJT T2_DJT T3_DJT T4_DJT T5_DJT])
display(['Th1out DJT  Th2out DJT Th3out DJT Th4out DJT Th5out DJT  [degC]'])
disp([Th1out_DJT Th2out_DJT Th3out_DJT Th4out_DJT Th5out_DJT])
display(['    w1 DJT    w2 DJT'])
disp([w1_DJT w2_DJT])
display(['  w1 ratio   w2 ratio [%]'])
disp([w1_rat_DJT w2_rat_DJT])
display([' DeltaT hot side '])
display(['    HX1        HX2       HX3       HX4       HX5   '])
disp([DeltaT_hot1_DJT DeltaT_hot2_DJT DeltaT_hot3_DJT DeltaT_hot4_DJT DeltaT_hot5_DJT])
display([' DeltaT cold side '])
display(['    HX1        HX2       HX3       HX4       HX5   '])
disp([DeltaT_cold1_DJT DeltaT_cold2_DJT DeltaT_cold3_DJT DeltaT_cold4_DJT DeltaT_cold5_DJT])