%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%                                              %%%%%
%%%%%               MONGSTAD 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 seven
% streams in parallel. The conditions are taken from
% Statoil Mongstad. The SOC split is found and compared
% to the real RTO split. 

clear all;
clc;

%% Defining all parameters

global T0;
global m0
global Tcold;
global Thot;

T0  = 132.909;
%% Initial flows (RTO)

m0A = 24.08753; % Stream A
m0B = 48.49833; % Stream B
m0C = 17.84011; % Stream C
m0D = 35.61500; % Stream D
m0E = 34.68861; % Stream E
m0F = 65.57833; % Stream F
m0G = 28.73583; % Stream G

% Total flow minus m0G
m0=m0A+m0B+m0C+m0D+m0E+m0F;

% Initial splits (without stream G)
u0=[m0A m0B m0C m0D m0E m0F]./m0;

%% RTO results (0 ('zero') indicates RTO)
[c,ceq0]=constraints(u0);
Tcold0=Tcold;
Thot0=Thot;

T1A0 = Tcold0(1); T2B0 = Tcold0(3); T2C0 = Tcold0(5);
T2D0 = Tcold0(7); T1E0 = Tcold0(8); T2F0 = Tcold0(10);

Tend0 = u0(1)*T1A0+u0(2)*T2B0+u0(3)*T2C0+u0(4)*T2D0+u0(5)*T1E0+u0(6)*T2F0;

%% Running fmincon to find optimal solution
options = optimset('Display','iter');

% Function to be minimized, J = 1
% Results: fmincon finds the 'u' that fulfils the constraints

u = fmincon(@(x)1,u0,[],[],[],[],0,1,@constraints,options);

% Calculating (intermediate) end temperature

Th1A = Thot(1); Th1B = Thot(2); Th2B = Thot(3); Th1C = Thot(4); Th2C = Thot(5); Th1D = Thot(6); Th2D = Thot(7); Th1E = Thot(8); Th1F = Thot(9); Th2F = Thot(10);
T1A = Tcold(1); T1B = Tcold(2); T2B = Tcold(3); T1C = Tcold(4); T2C = Tcold(5); T1D = Tcold(6); T2D = Tcold(7); T1E = Tcold(8); T1F = Tcold(9); T2F = Tcold(10);

Tend = u(1)*T1A+u(2)*T2B+u(3)*T2C+u(4)*T2D+u(5)*T1E+u(6)*T2F;

% New mass streams
m=u.*m0;
m(7) = m0G;

% Splits including stream G
unew = m./sum(m); % SOC splits
u0new=[m0A m0B m0C m0D m0E m0F m0G]./sum(m); % RTO splits

%% Temperature after O1

% Temperature after heat exchanger O1 (E153)
wc      = 615377.76;
n       = 1;
wh      = 186285.21;
Thin    = 282.608;
UA      = 144051.68;

[T Th] = tempCalc(Tend,Thin,wc,wh,UA,n);

% Temperature after mixing with stream G
cp0     = 2406.224915;
wc1G    = m0G*cp0;
n1G     = 1;
wh1G    = 289894.73;
Th1inG  = 167.487;
UA1G    = 96449.75;

[T1G Th1G] = tempCalc(T0,Th1inG,wc1G,wh1G,UA1G,n1G);

wc2G    = m0G*2622.608824;
n2G     = 1;
wh2G    = 317501.17;
Th2inG  = 262.586;
UA2G    = 125079.36;

[T2G Th2G] = tempCalc(T1G,Th2inG,wc2G,wh2G,UA2G,n2G);

TendFinal = (1-unew(7))*T + unew(7)*T2G;

Tfinal  = [T1A T2B T2C T2D T1E T2F T2G Tend TendFinal]';

%% RTO temperature

% Temperature after heat exchanger O1 (E153)

[Trto Thrto] = tempCalc(Tend0,Thin,wc,wh,UA,n);

% Temperature after mixing with stream G
TendFinal0 = (1-u0new(7))*Trto + u0new(7)*T2G;

Tfinal0 = [T1A0 T2B0 T2C0 T2D0 T1E0 T2F0 T2G Tend0 TendFinal0]';