%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%                                              %%%%%
%%%%%      PERSTORP STEADY-STATE CASE STUDY        %%%%%
%%%%%                                              %%%%%
%%%%%          Daniel Greiner Edvardsen            %%%%%
%%%%%                 IKP, NTNU                    %%%%%
%%%%%               June 14, 2011                  %%%%%
%%%%%                                              %%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Info: This m-file calculates the hot and cold outlet 
% temperatures for a heat exchanger network with three 
% streams; on the upper stream there are two heat
% exchangers in series, on the middle stream there is 
% one heat exchanger and on the lower stream there is 
% also one heat exchanger. The conditions are taken from
% Perstorp in Sweden. 
% The optimum is found (optimal split) that 
% results in the maximum end temperature. The SOC split
% is also found and compared with the optimal split and
% the actual split at Perstorp.

clc;
clear all;
close all;

% Which year? Alternatives: 2008.mat, 2009.mat, 2010.mat. Remember to change 'year'.mat in perstorpdata.m as well.
load 2010.mat;

% Which data point? k=1: first hour, k=8760: last hour
k=4000;

%% DEFINING PARAMETERS
% Loading relevant temperatures, mass streams, heat capacities, and densities from 'perstorpdata.m' All units in SI.
[T0 Tp Th11in Th1in Th2in Th3in Th11 Th1 Th2 Th3 T11_ T11 T1 T2 T3 m0 mh11 mh1 mh2 mh3 m11 m1 m2 m3 cp0 cp11 cp1 cp2 cp3 rho0 rho11 rho1 rho2 rho3 w0 w11 w1 w2 w3 wh11 wh1 wh2 wh3] = perstorp_data(k);

%% UA AND SPLIT RATIOS

% Calculating UA-values from Perstorp Data
UA11 = w11*(T11-T0)/((Th11-T0)-(Th11in-T11))*log((Th11-T0)/(Th11in-T11));
UA1  = w1*(T1-T11)/((Th1-T11)-(Th1in-T1))*log((Th1-T11)/(Th1in-T1));
UA2  = w2*(T2-T0)/((Th2-T0)-(Th2in-T2))*log((Th2-T0)/(Th2in-T2));
UA3  = w3*(T3-T0)/((Th3-T0)-(Th3in-T3))*log((Th3-T0)/(Th3in-T3));

% Split ratios
ptm1 = m1/m0; ptm2 = m2/m0; ptm3 = m3/m0;
sum = ptm1+ptm2+ptm3;
pt1 = ptm1/sum; pt2 = ptm2/sum; pt3 = ptm3/sum;
P = [pt1 pt2 pt3];

%% ITERATION PROCEDURE

% Number of iterations
iter=100;

U1=1/iter:1/iter:1;
U2=1/iter:1/iter:1;

n = zeros(iter,iter);
T11=n; T1=n; T2=n; T3=n; Th11=n; Th1=n; Th2=n; Th3=n; Tmix=n; e11=n; eh11=n; e1=n; eh1=n; e2=n; eh2=n; e3=n; eh3=n; C11=n; C1=n; C2=n; C3=n; NTU11=n; NTU1=n; NTU2=n; NTU3=n; U=n; CON=n; CHEN=n; CHEN2=n; CON2=n;

for i=1:iter
for j=1:iter

u1 = i/iter;
u2 = j/iter;

if(u1+u2<1) % Assuring u3=1-u1-u2 is never = 0 (making NTU3 = inf) nor < 0

% Calculating outlet temperatures and info about HEs (Only u1 and u2 are changing)
[T HE] = tempCalc3(T0,w0,UA11,UA1,UA2,UA3,Th11in,Th1in,Th2in,Th3in,wh11,wh1,wh2,wh3,u1,u2);

T11(i,j)=T(1); T1(i,j)=T(2); T2(i,j)=T(3); T3(i,j)=T(4); Th11(i,j)=T(5); Th1(i,j)=T(6); Th2(i,j)=T(7); Th3(i,j)=T(8); Tmix(i,j)=T(9);e11(i,j)=HE(1); eh11(i,j)=HE(2); e1(i,j)=HE(3); eh1(i,j)=HE(4); e2(i,j)=HE(5); eh2(i,j)=HE(6); e3(i,j)=HE(7); eh3(i,j)=HE(8);C11(i,j)=HE(9); C1(i,j)=HE(10); C2(i,j)=HE(11); C3(i,j)=HE(12); NTU11(i,j)=HE(13); NTU1(i,j)=HE(14); NTU2(i,j)=HE(15); NTU3(i,j)=HE(16);

else
break;
end

end 
end

%% Results

[v p]=max(Tmix);
[v2 p2]= max(max(Tmix));

p1=p(p2);

% OPTIMAL (maximum) temperature
Tmax = Tmix(p1,p2);

% OPTIMAL splits
u1 = p1/iter;
u2 = p2/iter;

% Controlled variables
c1 = ((Th1in-T1)./(Th11in-T0)-1).*(T11-T0).^2./(Th1in-T11)+(T1-T0).^2./(Th1in-T11)-(T2-T0).^2./(Th2in-T0);
c2 = (T3-T0).^2./(Th3in-T0)-(T2-T0).^2./(Th2in-T0);
   
%% Plot Self-Optimizing solution
close all;
%How many contour lines? Default: 20
v=fix(Tmax)-20:1:fix(Tmax);
%How many labels? Defalt: Every 5th
v2=fix(Tmax)-20:5:fix(Tmax);

figure
contour(U1,U2,double(c1),[0 0],'b','linewidth',2); hold on;
contour(U1,U2,double(c2),[0 0],'r','linewidth',2); hold on;

[C,h] = contour(U1,U2,double(Tmix), v, 'HandleVisibility','off'); hold on; 
th = clabel(C,h,v2);
set(th,'BackgroundColor',[1 1 .6],'Edgecolor',[.7 .7 .7])

[C2,h2]=contour(U1,U2,double(Tmix), round(Tp*10)./10, 'linewidth',3,'HandleVisibility','off'); hold on;
th2 = clabel(C2,h2);
set(th2,'BackgroundColor',[0 1 0],'Edgecolor',[.7 .7 .7])

scatter(u2,u1,'k','filled');

xlabel('u_2'); ylabel('u_1');
legend('c_{1}=0' ,'c_{2}=0','T_{opt}');
grid on