function [J,c,ceq,Wc,QHA_07,Cross,filename,THA_01]=f(w,p1,p2,p3)
%% Main Model code
% This code includes a steady state simulation of the whole mini LNG plant. 
% In this code all the models for the different units are solved
% simultaneously. The initialization comes from the files Main.m and
% Initialpoint.m. 

%%  Initialize variables

N_s=17;% Number of streams 
%Initialize variables
P    = zeros(1,N_s); % Pressures [pa]
T    = zeros(1,N_s); % Stream temperature [K]
n    = ones(1,N_s);  % Flows in [kmol/s]
h    = zeros(1,N_s); % Molar enthalpy [Kj/kmol]
VF   = zeros(1,N_s); % Vapor Fraction 
H    = zeros(1,N_s); % Total Enthalpy [Kj]
x    = zeros(N_s,5); % Liquid Composition
y    = zeros(N_s,5); % Vapor Composition
z    = zeros(N_s,5); % Total Composition
s    = zeros(N_s,5); % Liquid slack
beta = zeros(N_s,5); % Relaxation

 

% Set known compositions
z(1,:)  = [0.1074 0.3292 0.4096 0.1345 0.0193];
z(17,:) = [0.1074 0.3292 0.4096 0.1345 0.0193];


% Set known Temperatures
T(2) = 35 + 273.15;
T(6) = 273.15 - 75;
T(8) = 273.15 - 146.6;


% Set known vapor fractions
VF(1)  = 1;
VF(5)  = 1;
VF(4)  = 0;



% Define heat from the environment 
Q_env = 80.5852*p3;


%% Calculate the BOG requirements 

Tin_bog = -32 + 273.15+p1; %Inlet temperature [K]
m_bog=610*p2; %Gas mass flow [Kg/h] (Neskaa, 2010)
[ QHA_07] = BOGC( Tin_bog, m_bog );

%% Data extraction from the initial point vector 
%Define sets for different kinds of flashes 
TP  = [ 6 ];  %TP Flash
PH  = [ 3, 7, 11, 13, 14, 15]; %PH Flash
PHs = [ 16, 17];%PH Flash with slack variables
F = [TP, PH, PHs]; % Total set of flash calculations
FH= [PH,PHs]; % Set of all PH flash calculations

%Define important parameters
NC=5;

T(1) = w(1);
T1id = w(2);
a=2;
for i = F
    x(i,1:NC-1)=w(a+1:a+NC-1);
    x(i,NC)=1-sum(x(i,1:NC-1));
    y(i,1:NC-1)=w(a+NC:a+2*NC-2);
    y(i,NC)=1-sum(y(i,1:NC-1));
    VF(i)=w(a+2*NC-1);
    Zl(i)=w(a+2*NC);
    Zg(i)=w(a+2*NC+1);
    a=a+2*NC+1;
    if ismember(i,FH)
        T(i)=w(a+1); a=a+1;
        if ismember(i,PHs)
            s(i,:)=w(a+1:a+NC);
            beta(i,:)=w(a+NC+1:a+2*NC);
            a=a+2*NC;
        end
    end
end


T(10)= T(8);
T(9) = T(8);
THA_01 = T(3);
%% Optimization variables

% Set known flows 
N = 0.0434*w(a+1); %Refrigerant flow [Kmol/s]
n = N*n;

% Define split fraction
Sf = w(a+2);

% Set the pressures for the whole system
lp = [7 12 11 13 14 15 16 17];
hp = [1 2 3 4 5 6 8 9 10]; 
P(lp) = 3e5*w(a+3); 
P(hp) = 20e5*w(a+4);

% Define heat for HA-01
QHA_01 = 250*w(a+5);

%% Model Initial parameters

% The following section of the code includes the code for the model of the
% system, it sets all the equations of the system in terms of non linear 
% equality constraints for the optimization
n_e = 1; %Initial number of equality constrains
n_i = 1; %Initial number of inequality constrains
PH='PH';
TP='TP';
Slack='on ';
noSlack='off';
%% Compressor 

% Pressure relation calculation
rel = P(1)/P(17);
p = [0.0015   -0.0314    0.1911    0.4217]; %3rd degree polynomial. 
eta = polyval(p,rel);


% Isoentropic calculation
[Z(17),~,h(17),~,s17]= srks(z(17,:),T(17),P(17),2);    
[~,~,h1id,~,s1id]= srks(z(1,:),T1id,P(1),2);  
f{n_e}= (s17 - s1id)*10^-2; n_e=n_e+1;

% Real calculations
Wc =(h1id-h(17))/eta;           % Work calculation for 
h1real = h(17) + Wc;            % Real output enthalpy
[Z(1),~,h(1),~]  = srks(z(1,:),T(1),P(1),2);
Wc = Wc*n(17);
f{n_e}= (h1real - h(1))*10^-4; n_e=n_e+1;

%Assign variables
n(2)=n(1);
z(2,:)=z(1,:);
H(1) = h(1)*n(1);

%% MR Precooler

% Vapour only
z(2,:)=z(1,:);
y(2,:)=z(1,:);
[~,~,h(2),~]= srks(z(2,:),T(2),P(2),2);
H(2) = h(2)*n(2);
QMR = H(2)-H(1);
VF(2) = 0; 

%Assign values for the next stream
z(3,:) = z(2,:);
n(3) = n(2);

%% HA - 01
% Energy balance
H(3) = H(2)-QHA_01;
h(3) = H(3)/n(3);

% PH Flash without slack variables
[f{n_e},~,~] = flash(PH,noSlack,z(3,:),x(3,:),y(3,:),T(3),P(3),...
    VF(3),Zl(3),Zg(3),h(3));
n_e=n_e+1; %Increment constraints indeces

%Assign flows and compositions

%Vapour stream
z(5,:) = y(3,:);
y(5,:) = y(3,:);
n(5) = n(3)*VF(3);
T(5) = T(3);
[Z(5),~,h(5),~] = srks(z(5,:),T(5),P(5),2);
H(5)=h(5)*n(5);


%Liquid stream
z(4,:)=x(3,:);
x(4,:)=x(3,:);
n(4)=n(3)*(1-VF(3));
T(4)=T(3);
[Z(4),~,h(4),~] = srks(z(4,:),T(4),P(4),1);
H(4)=h(4)*n(4);

%% Valve FCV-01
% Isoenthalpic Valve
h(7) = h(4);
n(7) = n(4);
z(7,:) = z(4,:);

% PH flash without slack variables 
[f{n_e},~,~] = flash(PH,noSlack,z(7,:),x(7,:),y(7,:),T(7),P(7),...
    VF(7),Zl(7),Zg(7),h(7));
n_e=n_e+1; %Increment constraints indeces

H(7)=n(7)*h(7);


%% Heat exchanger HA-02 first side
%Calculate stream 6
z(6,:)=z(5,:);
n(6)=n(5);

% TP Flash without slack
[f{n_e},~,h(6)] = flash(TP,noSlack,z(6,:),x(6,:),y(6,:),T(6),P(6),...
    VF(6),Zl(6),Zg(6));
n_e=n_e+1; %Increment constraints indeces

% Energy Calculation
H(6)=h(6)*n(6);

%QHA_02 Duty 
QHA_02=H(5)-H(6);

%% Calculate HA-04 first side
%Calculate stream 7
z(8,:) = z(6,:);
x(8,:) = z(8,:);
n(8) = n(6);
VF(8) = 0; % Liquid only
[~,~,h(8),~] = srks(z(8,:),T(8),P(8),1);
H(8) = h(8)*n(8);
QHA_04 = H(6)-H(8);

%Split

%Stream 10
z(10,:) = z(8,:);
x(10,:) = x(8,:);
VF(10) = VF(8);
n(10) = n(8)*Sf;
h(10) = h(8);
H(10) = h(10)*n(10);
T(10) = T(8);

%Stream 9
z(9,:) = z(8,:);
x(9,:) = x(8,:);
VF(9) = VF(8);
n(9) = n(8)*(1-Sf);
h(9) = h(8);
H(9) = h(9)*n(9);
T(9) = T(8);

%% FCV-04 // FCV-05
%Assign values for stream 11
z(11,:)=z(9,:);
n(11)=n(9);
h(11)=h(9);

% PH flash without slack variables 
[f{n_e},~,~] = flash(PH,noSlack,z(11,:),x(11,:),y(11,:),T(11),P(11),...
    VF(11),Zl(11),Zg(11),h(11));
n_e=n_e+1; %Increment constraints indeces

% Total enthalpy 
H(11)=n(11)*h(11);

%Stream 12 - Properties are the same as for stream 11
z(12,:) = z(10,:);
n(12) = n(10);
h(12) = h(10);
T(12) = T(11);
x(12,:) = x(11,:);
y(12,:) = y(11,:);
VF(12) = VF(11);
H(12) = n(12)*h(12);

%% HA-04 Second side
n(13)=n(12);
H(13)=H(12)+QHA_04;
h(13)=H(13)/n(13);
z(13,:)=z(12,:);

%PH flash without slack variables
[f{n_e},~,~] = flash(PH,noSlack,z(13,:),x(13,:),y(13,:),T(13),...
    P(13),VF(13),Zl(13),Zg(13),h(13));
n_e=n_e+1; %Increment constraints indeces

H(13)=n(13)*h(13);

%% HA-07 Refrigerant side 
z(14,:)=z(11,:);
n(14)=n(11);
H(14)=H(11)-QHA_07;
h(14)=H(14)/n(14);

%PH flash without slack variables
[f{n_e},~,~] = flash(PH,noSlack,z(14,:),x(14,:),y(14,:),T(14),...
    P(14),VF(14),Zl(14),Zg(14),h(14));
n_e=n_e+1; %Increment constraints indeces


%% Mixer
%Balances
n(15)=n(7)+n(13)+n(14);
z(15,:)= (n(7).*z(7,:)+n(13).*z(13,:)+n(14).*z(14,:))/n(15);
H(15)=H(7)+H(13)+H(14);
h(15)=H(15)/n(15);

%PH flash without slack variables
[f{n_e},~,~] = flash(PH,noSlack,z(15,:),x(15,:),y(15,:),T(15),...
    P(15),VF(15),Zl(15),Zg(15),h(15));
n_e=n_e+1; %Increment constraints indeces


%% QHA_02 second side
z(16,:)=z(15,:);
n(16)=n(15);
H(16)=H(15)+QHA_02;
h(16)=H(16)/n(16);

%PH flash with slack variables
[f{n_e},fi{n_i},~,Pen1] = flash(PH,Slack,z(16,:),x(16,:),y(16,:),T(16),...
    P(16),VF(16),Zl(16),Zg(16),h(16),s(16,:),beta(16,:));
n_e=n_e+1; n_i=n_i+1; %Increment constraints indeces


%% Heat Transfer from the environment 
%It is calculated by carrying out an energy balance

H(17) = Q_env + H(16); % This corresponds to the heat from the environment
n(17) = n(16);
h(17) = H(17)/n(17);

%PH flash with slack variables
[f{n_e},fi{n_i},~,Pen2] = flash(PH,Slack,z(17,:),x(17,:),y(17,:),T(17),...
    P(17),VF(17),Zl(17),Zg(17),h(17),s(17,:),beta(17,:));
n_e=n_e+1; n_i=n_i+1; %Increment constraints indeces

x(1,:)=x(17,:);
y(1,:)=y(17,:);
VF(1)=VF(17);
%% Operational constratints

% Only vapour entering to the compressor
f{n_e} = 1-VF(17);n_e=n_e+1;
N0 = n_i;

% Constraints regarding pressure relation boundaries 2 <= rel <=10
fi{n_i} = rel - 10; n_i=n_i+1;% 31
fi{n_i} = 2 - rel;  n_i=n_i+1;% 32

% Boundaries for the splitflow
fi{n_i} = n(9)/(n(10)+n(9)) - 1;n_i=n_i+1; %33
fi{n_i} = 0 - n(9)/(n(10)+n(9));n_i=n_i+1; %34

% Constraints regarding pressure values

% Low pressure between 2 and 4 bar
fi{n_i} = 1.5/3 - w(a+3);n_i=n_i+1; % 35

% High pressure between 18 and 25 bar
fi{n_i} = w(a+4) - 25/20;n_i=n_i+1; % 36

% Stream 3 output T between -44 C and 35 C
fi{n_i} = T(3) - (35 + 273.15);n_i=n_i+1; % 37
fi{n_i} = (273.15 -44) - T(3);n_i=n_i+1; % 38

N1= n_i-1;
%% Temperature crosses
Cross=[];

%BOG Temperatures from BOG.m 
TBOGin = Tin_bog; %Inlet temperature [K]
TBOGout = -154 +273.15; %Outlet temperature [K]

% HA - 02
Cross(1) = -T(16) + T(5);
Cross(2) = -T(15) + T(6);

% HA - 04
Cross(3) = - T(13) + T(6);
Cross(4) = -T(12) + T(8);

% HA - 07
Cross(5) = - T(14) + TBOGin;
Cross(6) = - T(11) + TBOGout;

P_c=sum(abs(Cross));
%% Heat exchanger sizes
%Load calculated values
filename='HXSizing.mat';
load(filename)

%The following contraints make sure that the value for UA for each
%exchanger is constant 

%HA-02
[ UA_HA_02_c ] = hxsizing( QHA_02, T(5), T(6), T(15), T(16) );
f{n_e} = UA_HA_02-UA_HA_02_c;n_e=n_e+1;

%HA-04
[ UA_HA_04_c ] = hxsizing( QHA_04, T(6), T(8), T(12), T(13) );
f{n_e} = UA_HA_04-UA_HA_04_c;n_e=n_e+1;

%HA-07
[ UA_HA_07_c ] = hxsizing( -QHA_07,  T(11), T(14), TBOGin, TBOGout);
f{n_e} = (UA_HA_07-UA_HA_07_c);n_e=n_e+1;

%% Objective function 
W0 = 250; %Initial value for dimensionless work
M1 = 10^6;
M2 = 10^3;

J = (Wc/W0)+M1*(Pen1+Pen2)+M2*P_c;

%% Define the constraints vectors
const_i=[];
for i = 1:n_i-1
    const_i=[const_i;fi{i}];
end

const_e=[];
for i = 1:n_e-1 
    const_e=[const_e;f{i}];
end
c=[const_i]; ceq=[const_e];



