%% Model to simulate a steady state HEN
%  Optimal Design

% Topology to be investigated:

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

close all;
clear all;
clc;

%% Parameters

par.w0  = 95;   %[kW/K] w= miCpi
par.wh1 = 60;   %[kW/K] 
par.wh2 = 27;   %[kW/K]
par.wh3 = 65;   %[kW/K]

par.Th1 = 203;  %[degC]
par.Th2 = 255;  %[degC]
par.Th3 = 248;  %[degC]
par.T0 = 130;   %[degC]

par.DeltaTmin = 0.5; %[degC]
par.n = 0.65;
% par.c0 = 1.4; %[$/m2] 
par.sc.x = [200*ones(7,1);45;45;5000*ones(3,1);200*ones(3,1)];
par.sc.j = 200;

%% Optimization

% x0 = [ Tend T1   T2   T3 Th1out Th2out  Th3out    w1   w2 Q1    Q2     Q3       UA1     UA2     UA3]

x0   = [ 225  168  199  208  154    150    145     42.5   42.5 2896   2323   6808      3        2       9]';

x0 = x0./par.sc.x;
A = []; b = []; Aeq = []; Beq = [];
LB = 0.00*ones(15,1); UB = inf*ones(15,1);
options = optimset('Algorithm','interior-point','display','iter'...
                   ,'MaxFunEvals',9000,'TolCon',1e-10,'TolX',1e-10);

%% FIXED INLET TEMEPRATURE
% [x,J,exitflag] = fmincon(@(x)Object_21_HXD(x,par),x0,A,b,Aeq,Beq,LB,UB,@(x)HEN_Constraints_21_HXD(x,par),options);
% exitflag
% x=x.*par.sc.x; % Unscale variables
% 
% Tend = x(1); T1 = x(2); T2 = x(3); T3 = x(4); Th1out = x(5); Th2out = x(6); Th3out = x(7); w1 = x(8);
% w2 = x(9); Q1 = x(10); Q2 = x(11); Q3 = x(12); UA1 = x(13); UA2 = x(14); UA3 = x(15);
% T0 = par.T0; Th1 = par.Th1; Th2 = par.Th2; Th3 = par.Th3;
% 
% display(['   T1         T2         T3       Th1out    Th2out    Th3out   [degC]'])
% disp([T1 T2 T3   Th1out Th2out Th3out])
% 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     [Wm2/K]'])
% disp([UA1 UA2 UA3])
% 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(['  Q1            Q2         Q3'])
% disp([Q1 Q2 Q3])

%% RESULTS WITH VARIATION OF C0
% Defining the initial points 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 = [];
Tend_vec = [];
w1_vec = [];
w2_vec = [];
UA1_vec = [];
UA2_vec = [];
UA3_vec = [];
Th1out_vec = [];
Th2out_vec = [];
Th3out_vec = [];
Q1_vec = [];
Q2_vec = [];
Q3_vec = [];
exitflag_vec = [];
DeltaTHX1_vec = [];
DeltaTHX2_vec = [];
DeltaTHX3_vec = [];

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

Tend = x(1); T1 = x(2); T2 = x(3); T3 = x(4); Th1out = x(5); Th2out = x(6); Th3out = x(7); w1 = x(8);
w2 = x(9); Q1 = x(10); Q2 = x(11); Q3 = x(12); UA1 = x(13); UA2 = x(14); UA3 = x(15);


w0 = par.w0; 
wh1 = par.wh1; wh2 = par.wh2; wh3 = par.wh3;
Th1 = par.Th1; Th2 = par.Th2; Th3 = par.Th3;
T0 = par.T0;

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;

T1_vec(i)= T1;
T2_vec(i) = T2;
T3_vec(i) = T3;
Tend_vec(i) = Tend;
w1_vec(i) = w1;
w2_vec(i) = w2;
UA1_vec(i) = UA1;
UA2_vec(i) = UA2;
UA3_vec(i) = UA3;
exitflag_vec(i) = exitflag;
Th1out_vec(i) = Th1out;
Th2out_vec(i) = Th2out;
Th3out_vec(i) = Th3out;
Q1_vec(i) = Q1;
Q2_vec(i) = Q2;
Q3_vec(i) = Q3;
DeltaTHX1_vec(i) = DeltaTHX1;
DeltaTHX2_vec(i) = DeltaTHX2;
DeltaTHX3_vec(i) = DeltaTHX3;
c0_vec(i) = c0;

c0 = c0+DeltaC0;

end

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(3,1,1)
plot(c0_vec,T1_vec,'LineWidth',2)
xlabel('c0 cost factor')
ylabel('T1 [^{\circ}C]')
legend('T1')

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

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

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

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

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

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

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

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

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

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

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


