%% 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/degC] w= miCpi
par.wh1 = 50;    %[kW/degC] 
par.wh2 = 30;    %[kW/degC]
par.wh3 = 15;    %[kW/degC]
par.wh4 = 70;    %[kW/degC]
par.wh5 = 20;

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

par.DeltaTmin = 0.00001; %[degC]
par.n = 0.65;
% par.c0 = 4; %[$/m2] 
par.sc.x = [200*ones(11,1);100;100;1000*ones(5,1);5*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   = [175 158 173 184 185 198 158 174 180 144 186 60 120 1564 854 590 5268 1233 1.52 1.90 1.59 9.72 2.18]';

x0 = x0./par.sc.x;
A = []; b = []; Aeq = []; Beq = [];
LB = 0.00*ones(23,1); UB = inf*ones(23,1);

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


% [x,J,exitflag] = fmincon(@(x)Object_32_HXD(x,par),x0,A,b,Aeq,Beq,...
%                  LB,UB,@(x)HEN_Constraints_32_HXD(x,par),options);
% exitflag
% x=x.*par.sc.x; % Unscale variables

% %% RESULTS
% 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); UA1 = x(19); UA2 = x(20); UA3 = x(21);
% UA4 = x(22); UA5 = x(23); % T0 = par.T0; Th1 = par.Th1;
% Th2 = par.Th2; Th3 = par.Th3; Th4 = par.Th4; Th5 = par.Th5;
% 
% 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])

%% VARIATION AND TRENDS
% Defining the initial point and step
c0 = 1;
co_end = 5;
DeltaC0 = 0.1;
n = (co_end-c0)/DeltaC0;
c0_vec = [];

% Utilizing the resulting variables
T1_vec = [];
T2_vec = [];
T3_vec = [];
T4_vec = [];
T5_vec = [];
Tend_vec = [];
w1_vec = [];
w2_vec = [];
UA1_vec = [];
UA2_vec = [];
UA3_vec = [];
UA4_vec = [];
UA5_vec = [];
Th1out_vec = [];
Th2out_vec = [];
Th3out_vec = [];
Th4out_vec = [];
Th5out_vec = [];
Q1_vec = [];
Q2_vec = [];
Q3_vec = [];
Q4_vec = [];
Q5_vec = [];
exitflag_vec = [];
DeltaTHX1_vec = [];
DeltaTHX2_vec = [];
DeltaTHX3_vec = [];
DeltaTHX4_vec = [];
DeltaTHX5_vec = [];

for i=1:n;
    
    [x,J,exitflag] = fmincon(@(x)Object_32_HXD(x,par,c0),...
        x0,A,b,Aeq,Beq,LB,UB,@(x)HEN_Constraints_32_HXD(x,par),options);
        x=x.*par.sc.x; % Unscale variables
exitflag

% States
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); UA1 = x(19); UA2 = x(20); 
UA3 = x(21); UA4 = x(22); UA5 = x(23); T0 = par.T0;
Th1 = par.Th1; Th2 = par.Th2; Th3 = par.Th3; Th4 = par.Th4; Th5 = par.Th5;

% Parameters
w0 = par.w0; 
wh1 = par.wh1; wh2 = par.wh2; wh3 = par.wh3;

% Calculating the DeltaT's for each HX
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;

Thot4 = Th4-T4;
Tcold4 = Th4out-T0;
DeltaTHX4 = Thot4/Tcold4;

Thot5 = Th5-T5;
Tcold5 = Th5out-T4;
DeltaTHX5 = Thot5/Tcold5;

% Inserting the calulcated value in the solution matrices
T1_vec(i)= T1;
T2_vec(i) = T2;
T3_vec(i) = T3;
T4_vec(i) = T4;
T5_vec(i) = T5;
Tend_vec(i) = Tend;
w1_vec(i) = w1;
w2_vec(i) = w2;
UA1_vec(i) = UA1;
UA2_vec(i) = UA2;
UA3_vec(i) = UA3;
UA4_vec(i) = UA4;
UA5_vec(i) = UA5,
exitflag_vec(i) = exitflag;
Th1out_vec(i) = Th1out;
Th2out_vec(i) = Th2out;
Th3out_vec(i) = Th3out;
Th4out_vec(i) = Th4out;
Th5out_vec(i) = Th5out;
Q1_vec(i) = Q1;
Q2_vec(i) = Q2;
Q3_vec(i) = Q3;
Q4_vec(i) = Q4;
Q5_vec(i) = Q5;
DeltaTHX1_vec(i) = DeltaTHX1;
DeltaTHX2_vec(i) = DeltaTHX2;
DeltaTHX3_vec(i) = DeltaTHX3;
DeltaTHX4_vec(i) = DeltaTHX4;
DeltaTHX5_vec(i) = DeltaTHX5;
c0_vec(i) = c0;

c0 = c0+DeltaC0;

end

% Plotting the results
figure(1)
plot(c0_vec,Tend_vec,'LineWidth',2)
xlabel('c0 cost factor')
ylabel('T_{end} [^{\circ}C]')
legend('T_{end}')

figure(2)
plot(c0_vec,exitflag_vec)
xlabel('c0')
ylabel('exitflag')

figure(3)
plot(c0_vec,w1_vec,'b','LineWidth',2)
hold on
plot(c0_vec,w2_vec,'r','LineWidth',2)
xlabel('c0 cost factor')
ylabel('w_{i} = mC_{p}')
legend('w_{1}','w_{2}')

figure(4)
subplot(5,1,1)
plot(c0_vec,T1_vec,'LineWidth',2)
xlabel('c0 cost factor')
ylabel('T1 [^{\circ}C]')
legend('T1')

subplot(5,1,2)
plot(c0_vec,T2_vec,'LineWidth',2)
xlabel('c0 cost factor')
ylabel('T2 [^{\circ}C]')
legend('T2')

subplot(5,1,3)
plot(c0_vec,T3_vec,'LineWidth',2)
xlabel('c0 cost factor')
ylabel('T3 [^{\circ}C]')
legend('T3')

subplot(5,1,4)
plot(c0_vec,T4_vec,'LineWidth',2)
xlabel('c0 cost factor')
ylabel('T4 [^{\circ}C]')
legend('T4')

subplot(5,1,5)
plot(c0_vec,T5_vec,'LineWidth',2)
xlabel('c0 cost factor')
ylabel('T5 [^{\circ}C]')
legend('T5')

figure(5)
subplot(5,1,1)
plot(c0_vec,Th1out_vec,'LineWidth',2)
xlabel('c0 cost factor')
ylabel('Th1out [^{\circ}C]')
legend('Th1out')

subplot(5,1,2)
plot(c0_vec,Th2out_vec,'LineWidth',2)
xlabel('c0 cost factor')
ylabel('Th2out [^{\circ}C]')
legend('Th2out')

subplot(5,1,3)
plot(c0_vec,Th3out_vec,'LineWidth',2)
xlabel('c0 cost factor')
ylabel('Th3out [^{\circ}C]')
legend('Th3out')

subplot(5,1,4)
plot(c0_vec,Th4out_vec,'LineWidth',2)
xlabel('c0 cost factor')
ylabel('Th4out [^{\circ}C]')
legend('Th4out')

subplot(5,1,5)
plot(c0_vec,Th5out_vec,'LineWidth',2)
xlabel('c0 cost factor')
ylabel('Th5out [^{\circ}C]')
legend('Th5out')

figure(6)
subplot(5,1,1)
plot(c0_vec,DeltaTHX1_vec)
xlabel('c0 cost factor')
ylabel('DeltaT HX1')
legend('DeltaT HX1')

subplot(5,1,1)
plot(c0_vec,DeltaTHX2_vec)
xlabel('c0 cost factor')
ylabel('DeltaT HX2')
legend('DeltaT HX2')

figure(7)
subplot(5,1,1)
plot(c0_vec,UA1_vec,'LineWidth',2)
xlabel('c0 cost factor')
ylabel('UA1 [kW/[^{\circ}C]]')
legend('UA1')

subplot(5,1,2)
plot(c0_vec,UA2_vec,'LineWidth',2)
xlabel('c0 cost factor')
ylabel('UA2 [kW/[^{\circ}C]]')
legend('UA2')

subplot(5,1,3)
plot(c0_vec,UA3_vec,'LineWidth',2)
xlabel('c0 cost factor')
ylabel('UA3 [kW/[^{\circ}C]]')
legend('UA3')

subplot(5,1,4)
plot(c0_vec,UA4_vec,'LineWidth',2)
xlabel('c0 cost factor')
ylabel('UA4 [kW/[^{\circ}C]]')
legend('UA4')

subplot(5,1,5)
plot(c0_vec,UA5_vec,'LineWidth',2)
xlabel('c0 cost factor')
ylabel('UA5 [kW/[^{\circ}C]]')
legend('UA5')

% figure(8)
% plot(c0_vec,UA4_vec,'LineWidth',2)
% axis([1 5 0 10])
% xlabel('c0 cost factor')
% ylabel('UA4 [kW/[^{\circ}C]]')
% legend('UA4')
% 
% figure(9)
% plot(c0_vec,UA1_vec,'LineWidth',2)
% axis([1 5 0 2])
% xlabel('c0 cost factor')
% ylabel('UA1 [kW/[^{\circ}C]]')
% legend('UA1')



