% Optimal operation of a 1:1 HEN and 
% operation using the Jaeschke temperature.
% Simulations are based on the NTU-method


% Topology to be investigated

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%               1                %
%        -------0-------         %
%   ----|               |----    %
%        -------0-------         %
%               2                %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


clc;
clear all;
close all;

% Defining parameters

% Cases evaluated
% Vector parameters: [T0 w0 wh1 wh2 Th1in Th2in UA1 UA2]

caseI        = [130 100  50 50  203 248  10 30];
caseII       = [130 100 50 50  203 248  31.1 93.9];
caseIII      = [130 50 100 100  203 248  10 30];
caseIV       = [130 100 50 50 203 248 100 300];
caseV        = [130 100 400 100 203 248 1000 100];
caseVI       = [130 100 400 100 203 248 1000 1000];

% Select case
casesel = caseI;

% Operation parameters
T0      = casesel(1); % Feed stream temperature [degC] 	
w0      = casesel(2); % [kW/K]

% Utility parameters
wh1     = casesel(3); % Hot stream 1 Heat Capacity rate [kW/K]
wh2 	= casesel(4); % Hot stream 2 Heat Capacity rate [kW/K]
Th1in	= casesel(5); % Hot stream 1 Temperature [degC]
Th2in 	= casesel(6); % Hot stream 2 Temperature [degC]

% Design parameters
UA1 	= casesel(7);       % [kW/K]
UA2 	= casesel(8);       % [kW/K]

% Number of iterations
N=1000;

n = zeros(N,1);
T1=n; T2=n; Th1=n; Th2=n; Tmix=n; e1=n; eh1=n; e2=n; eh2=n; 
C1=n; C2=n; NTU1=n; NTU2=n; U=n;

% Calculating HX based on the NTU-method for all splits ranging [0,1]:
for i=1:N
    
    u = i/N;
    U(i)=u;

%     Calculating outlet temperatures and info about HEs 
%     (only u is changing)
    [T HE] = TempCalc(T0,w0,UA1,UA2,Th1in,wh1,Th2in,wh2,u);
     
    T1(i)=T(1); T2(i)=T(2); Th1(i)=T(3); Th2(i)=T(4); 
    Tmix(i)=T(5); e1(i)=HE(1); eh1(i)=HE(2); e2(i)=HE(3); 
    eh2(i)=HE(4); C1(i)=HE(5); C2(i)=HE(6); NTU1(i)=HE(7); 
    NTU2(i)=HE(8);
    
end


% RESULTS

% Finding optimal split
[Tmixm,nr]=max(Tmix);

split=U(nr);
T1m=T1(nr);
Th1m=Th1(nr);
T2m=T2(nr);
Th2m=Th2(nr);
Tmixm
split

% Finding the self-optimizing split

% Jaeschke Temperature for HX1 and HX2
JT = (T1-T0).^2./(Th1in-T0) - (T2-T0).^2./(Th2in-T0);

[JTmin,nr2]=min(abs(JT));

JT_opt=JT(nr);
JTsplit=U(nr2);
T1JT=T1(nr2);
Th1JT=Th1(nr2);
T2JT=T2(nr2);
Th2JT=Th2(nr2);
JTmin;
JTTmax=Tmix(nr2);
JTTmax
JTsplit

% Jaeschke temperature in the presence of measurement errors, max

JTTmax_vec = [];
TempLoss = [];

nT0 = 0;
nTh1 = 0;
nTh2 = 0;
nT1 = 0;
nT2 = 0;
 
M = 1000;

% Simulating HX with measurement errors, with given Measurement errors
% (data file Measurement_Errors.m)
for j=1:M;
    
    load Measurement_Errors
    
    nT0 = noise(1,j);
    nTh1 = noise(2,j);
    nTh2 = noise(3,j);
    nT1 = noise(4,j);
    nT2 = noise(5,j);

    % Implementing the noise in the control variable   
    JT_noise = ((T1+nT1)-(T0+nT0)).^2./((Th1in+nTh1)-(T0+nT0)) - ((T2+nT2)-(T0+nT0)).^2./((Th2in+nTh2)-(T0+nT0));

    [JTmin_noise,nr3] = min(abs(JT_noise));
    JT_noise_split = U(nr3);
    JTnoiseTmax = Tmix(nr3);

    JTTmax_vec(j) = JTnoiseTmax;
    TempLoss(j) = Tmixm-JTnoiseTmax;

    noise(:,j) = [nT0, nTh1, nTh2, nT1, nT2]';

end

% Worst case loss and avergae loss
WCloss = max(TempLoss);
AVGloss = sum(TempLoss)/M;
WCloss
AVGloss

% % Calculating temperature difference on each side of each HX

dTcold1=Th1-T0;
dThot1=Th1in-T1; 

dTcold2=Th2-T0; 
dThot2=Th2in-T2;

% Calculating errors from AMTD approximation
[eU1 eU2 eAM1 eAM2] = ErrorCalc(dTcold1, dThot1, dTcold2, dThot2);

% Calculating the AMTD approximation valid range..
theta1 = dThot1./dTcold1;
theta2 = dThot2./dTcold2;


% PLOTTING THE RESULTS

% Temperature and control variable profile with split u
% return
h = figure;
% return
% figure(1)
y1start = 160; y1end = 210; y1step = 10;
y2start = -60; y2end = 60; y2step = 60;

split = [split split];
JTs = [JTsplit JTsplit];
y11 = [y1start y1end];
y22 = [y2start y2end];

[AX,H1,H2] = plotyy(U,JT,U,Tmix);
set(get(AX(2),'Ylabel'),'String',...
    'T_{end} [ \circC]','fontsize',12) 
set(get(AX(1),'Ylabel'),'String',...
    'Controlled variable, JT [ \circC]','fontsize',12)
axis(AX(2),[0 1 y1start y1end]);
axis(AX(1),[0 1 y2start y2end]);
set(AX(2),'YLim',[y1start y1end])
set(AX(2),'YTick',y1start:y1step:y1end)
set(AX(1),'YLim',[y2start y2end])
set(AX(1),'YTick',y2start:y2step:y2end)
set(H1,'linewidth',2)
set(H2,'linewidth',2)
xlabel('Split, u','fontsize',12);
hold on; 
H3 = plot(JTs,y22,'Color','k','LineStyle','--','LineWidth',2);
hold on
H4 = plot(split,y22,'Color','r','LineStyle','--','LineWidth',2);

set(H3,'parent',AX(1));
% hold on;
grid on;
print(h,'-depsc','CaseIId_optCalc.eps');


% Validity of AMTD approximation
UB = [1.4 1.4]; % Upper AMTD limit
LB = [(1/1.4) (1/1.4)]; % Lower AMTD limit
s = [0 1];

k = figure;
% figure(6);
plot(U,theta2,U,theta1,'LineWidth',2); 
xlabel('Split, u','fontsize',12);
ylabel('\theta_{1}/\theta_{2}','fontsize',12);
% legend('HX_{1,2}','HX_{1,1}','fontsize',12);
axis([0 1 0 2]);
% Using hline.m to include upper and lower bounds:
hline([1/1.4 1.4],{'m','m'},{'AMTD LB','AMTD UB'})
hold on; 
plot(splitline,solid,'Color','r','LineStyle','--','LineWidth',2);
legend('HX_{1,2}','HX_{1,1}','Optimal split','fontsize',11);
print(k,'-depsc','AMTD_CaseIIb.eps');




