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

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

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 = 40;   %[kW/degC]
par.wh6 = 35;   %[kW/degC]
par.wh7 = 30;   %[kW/degC]

% Hot stream inlet temperature
par.Th1 = 190;  %[degC]
par.Th2 = 203;  %[degC]
par.Th3 = 220;  %[degC]
par.Th4 = 235;  %[degC]
par.Th5 = 240;  %[degC]
par.Th6 = 245;  %[degC]
par.Th7 = 225;  %[degC]

% Cold stream inlet temperture
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]
par.UA6 = 8;    %[kWm2/degC]
par.UA7 = 11;   %[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]
par.P6 = 1.7;   %[$/kW]
par.P7 = 1.5;   %[$/kW]

% Scaling vector
par.sc.x = [200*ones(15,1);100;100;500*ones(7,1)];

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


%% OPTIMAL OPERATION

% Guessing outlet variables
% x0 = [Tend T1 T2 T3 T4 T5 T6 T7 Th1 Th2 Th3 Th4 Th5 Th6 Th7 w1 w2 ...
%       Q1 Q2 Q3 Q4 Q5 Q6 Q7]
x0 = [148 131 133 138 138 140 145 150 188 198 200 215 190 230 200 50 50 ...
      59 137 297 333 200 250 300]';

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

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

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

% RESULTS
% Outlet temperatures
Tend = x(1);
T1 = x(2); T2 = x(3); T3 = x(4); T4 = x(5); T5 = x(6); T6 = x(7); T7 = x(8);
Th1out = x(9); Th2out = x(10); Th3out = x(11); Th4out = x(12);
Th5out = x(13); Th6out = x(14); Th7out = x(15);
% Split
w1 = x(16); w2 = x(17);
% Heat transfer
Q1 = x(18); Q2 = x(19); Q3 = x(20); Q4 = x(21); Q5 = x(22);
Q6 = x(23); Q7 = x(24);
% 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_hot6 = Th6 - T6;
DeltaT_hot7 = Th7 - T7;
DeltaT_cold1 = Th1out - T0;
DeltaT_cold2 = Th2out - T1;
DeltaT_cold3 = Th3out - T2;
DeltaT_cold4 = Th4out - T3;
DeltaT_cold5 = Th5out - T4;
DeltaT_cold6 = Th6out - T5;
DeltaT_cold7 = Th7out - T0;

% Displaying the results
display(['  Tend [degC] = '])
disp(Tend)
display(['   T1          T2         T3        T4        T5        T6        T7   [degC]'])
disp([T1 T2 T3 T4 T5 T6 T7])
display(['    Th1out    Th2out    Th3out    Th4out    Th5out    Th6out    Th7out   [degC]'])
disp([Th1out Th2out Th3out Th4out Th5out Th6out Th7out])
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       HX6       HX7   '])
disp([DeltaT_hot1 DeltaT_hot2 DeltaT_hot3 DeltaT_hot4 DeltaT_hot5 DeltaT_hot6 DeltaT_hot7])
display([' DeltaT cold side '])
display(['    HX1        HX2       HX3       HX4       HX5       HX6       HX7   '])
disp([DeltaT_cold1 DeltaT_cold2 DeltaT_cold3 DeltaT_cold4 DeltaT_cold5 DeltaT_cold6 DeltaT_cold7])


%% OPERATION USING THE JAESCHKE TEMPERATURE

% Guessing outlet variables
% x0 = [Tend T1 T2 T3 T4 T5 T6 T7 Th1 Th2 Th3 Th4 Th5 Th6 Th7 w1 w2 ...
%       Q1 Q2 Q3 Q4 Q5 Q6 Q7]
x0 = [148 131 133 138 138 140 145 150 188 198 200 215 190 230 200 50 50 ...
      59 137 297 333 200 250 300]';
  
% Defining parameters  
Th1 = par.Th1; Th2 = par.Th2; Th3 = par.Th3; Th4 = par.Th4; Th5 = par.Th5;
Th6 = par.Th6; Th7 = par.Th7;
T0 = par.T0;


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

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

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


% 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); T6_DJT = xDJT(7); T7_DJT = xDJT(8);
Th1out_DJT = xDJT(9); Th2out_DJT = xDJT(10); Th3out_DJT = xDJT(11); Th4out_DJT = xDJT(12);
Th5out_DJT = xDJT(13); Th6out_DJT = xDJT(14); Th7out_DJT = xDJT(15);
% Split
w1_DJT = xDJT(16); w2_DJT = xDJT(17);
% Heat transfer
Q1_DJT = xDJT(18); Q2_DJT = xDJT(19); Q3_DJT = xDJT(20); Q4_DJT = xDJT(21); Q5_DJT = xDJT(22);
Q6_DJT = xDJT(23); Q7_DJT = xDJT(24);
% 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_hot6_DJT = Th6 - T6_DJT;
DeltaT_hot7_DJT = Th7 - T7_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 - T4_DJT;
DeltaT_cold6_DJT = Th6out_DJT - T5_DJT;
DeltaT_cold7_DJT = Th7out_DJT - T0;

% Displaying the results
display(['  Tend [degC] = '])
disp(Tend)
display(['   T1          T2         T3        T4        T5        T6        T7   [degC]'])
disp([T1_DJT T2_DJT T3_DJT T4_DJT T5_DJT T6_DJT T7_DJT])
display(['    Th1out    Th2out    Th3out    Th4out    Th5out    Th6out    Th7out   [degC]'])
disp([Th1out_DJT Th2out_DJT Th3out_DJT Th4out_DJT Th5out_DJT Th6out_DJT Th7out_DJT])
display(['      w1        w2'])
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       HX6       HX7   '])
disp([DeltaT_hot1_DJT DeltaT_hot2_DJT DeltaT_hot3_DJT DeltaT_hot4_DJT DeltaT_hot5_DJT DeltaT_hot6_DJT DeltaT_hot7_DJT])
display([' DeltaT cold side '])
display(['    HX1        HX2       HX3       HX4       HX5       HX6       HX7   '])
disp([DeltaT_cold1_DJT DeltaT_cold2_DJT DeltaT_cold3_DJT DeltaT_cold4_DJT DeltaT_cold5_DJT DeltaT_cold6_DJT DeltaT_cold7_DJT])






