% HEN_Constraints function 4:1 HEN for simulations with the Jaeschke temp

% Nonlinear constraints for optimizing a HEN
% Includes mass, energy and steady state balances and the Jaeschke temp


function [Cineq, Res] = HEN_Constraints_41_DJT(x,par)

% Defining state variables
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);

% Defining parameteres
w0 = par.w0; 
wh1 = par.wh1; wh2 = par.wh2; wh3 = par.wh3; wh4 = par.wh4; wh5 = par.wh5;
Th1 = par.Th1; Th2 = par.Th2; Th3 = par.Th3; Th4 = par.Th4; Th5 = par.Th5;
T0 = par.T0;
UA1 = par.UA1; UA2 = par.UA2; UA3 = par.UA3; UA4 = par.UA4; UA5 = par.UA5;
DeltaTmin = par.DeltaTmin;
P1 = par.P1; P2 = par.P2; P3 = par.P3; P4 = par.P4; P5 = par.P5;



%% INEQUALITY CONSTRAINTS

% HX1
Cineq1 = -(Th1-T1-DeltaTmin); % HOT SIDE HX1
Cineq2 = -(Th1out-T0-DeltaTmin); % COLD SIDE HX1

% HX2
Cineq3 = -(Th2-T2-DeltaTmin); % HOT SIDE HX2
Cineq4 = -(Th2out-T1-DeltaTmin); % COLD SIDE HX2

% HX3
Cineq5 = -(Th3-T3-DeltaTmin); % HOT SIDE HX3
Cineq6 = -(Th3out-T2-DeltaTmin); % COLD SIDE HX3

% HX4
Cineq7 = -(Th4-T4-DeltaTmin); % HOT SIDE HX4
Cineq8 = -(Th4out-T3-DeltaTmin); % COLD SIDE HX4

% HX 5
Cineq9 = -(Th5-T5-DeltaTmin); % HOT SIDE HX5
Cineq10 = -(Th5out-T0-DeltaTmin); % COLD SIDE HX5

Cineq = [Cineq1;Cineq2;Cineq3;Cineq4;Cineq5;Cineq6;Cineq7;Cineq8;...
        Cineq9;Cineq10];
Cineq = [];



%% MODEL EQUATIONS

% % AMTD
% DeltaT1 = 0.5*((Th1out-T0)+(Th1-T1));
% DeltaT2 = 0.5*((Th2out-T1)+(Th2-T2));
% DeltaT3 = 0.5*((Th3out-T2)+(Th3-T3));
% DeltaT4 = 0.5*((Th4out-T3)+(Th4-T4));
% DeltaT5 = 0.5*((Th5out-T0)+(Th5-T5));

%UNDERWOOD APPROXIMATION
 DeltaT1 = ((((Th1out-T0)^1/3)+((Th1-T1)^1/3))/2)^3;
 DeltaT2 = ((((Th2out-T1)^1/3)+((Th2-T2)^1/3))/2)^3;
 DeltaT3 = ((((Th3out-T2)^1/3)+((Th3-T3)^1/3))/2)^3;
 DeltaT4 = ((((Th4out-T3)^1/3)+((Th4-T4)^1/3))/2)^3;
 DeltaT5 = ((((Th5out-T0)^1/3)+((Th5-T5)^1/3))/2)^3;

 
 
%% JAESCHKE TEMPERATURES
% Upper path
JT11 = P1*(T1-T0)^2/(Th1-T0);
JT12 = P2*((T2-T1)*(T2+T1-2*T0-JT11))/(Th2-T1);
JT13 = P3*((T3-T2)*(T3+T2-2*T0-JT12))/(Th3-T2);
JT14 = P4*((T4-T3)*(T4+T3-2*T0-JT13))/(Th4-T3);
% Lower path
JT21 = P5*(T5-T0)^2/(Th5-T0);



%% EQUALITY CONSTRAINTS
Res = [ % Upper path, 1st HX
         Q1-(w1*(T1-T0));                % Cold Stream, w1
         Q1+(par.wh1*(Th1out-Th1));      % Hot Stream, wh1   
         Q1-(UA1*DeltaT1);               % HX Design Equation
         
             
        % Upper path, 2nd HX
         Q2-(w1*(T2-T1));                % Cold Stream, w1
         Q2+(par.wh2*(Th2out-Th2));      % Hot Stream, wh2
         Q2-(UA2*DeltaT2);               % HX Design Equation
         
        
        % Upper path, 3rd HX
         Q3-(w1*(T3-T2));                % Cold Stream, w1
         Q3+(par.wh3*(Th3out-Th3));      % Hot Stream, wh3
         Q3-(UA3*DeltaT3);               % HX Design equation
         
        % Lower path, 4th HX
         Q4-(w1*(T4-T3));                % Cold stream, w2
         Q4+(par.wh4*(Th4out-Th4));      % Hot stream, wh4
         Q4-(UA4*DeltaT4);               % HX design equation
         
         % Lower path, 5th HX
         Q5-(w2*(T5-T0));                % Cold stream, w2
         Q5+(par.wh5*(Th5out-Th5));      % Hot stream, wh4
         Q5-(UA5*DeltaT5);               % HX design equation
         
        % Mass balance
         w1+w2-w0;

        % Energy balance
         (w0*Tend)-(w1*T4)-(w2*T5);
     
        % Jaeschke temperature
          (JT11+JT12+JT13+JT14)-JT21];
      
end