function f=cola_lv_fun(x,p,u0)

%**************************************************************
%                   OBJECT FUNCTION

% Subroutine for self optimisation program optmisation.m
% for optimisation of multi-effect distilaltion columns
% Contains the object function
%
%**************************************************************
global Y S

%x = S*x;         % Scaling back to "original" values

% Number of stages in columns:
NT = 16;              % HP column
Nt = 26;              % LP column

% Inputs, disturbances and states:
LT_HP = x(1);                  % Reflux in HP column
VB_HP = x(2);                  % Boilup in HP column
A     = x(3);                  % Condensed flow, HP column
X1 = x(4:3*NT+3);              % States HP column

LT_LP = x(3*NT+4);             % Reflux in LP column
VD_LP = x(3*NT+5);             % Condensed flow, LP column
X2 = x(3*NT+6:3*NT+5+3*Nt);    % States in LP column

% Prices for object function (from main program):
PD = p(4);            % Distillate price
PV = p(5);            % Boilup price

% P-controllers for control of reboiler and condenser hold-up:
KcB_HP = -60; KcD_HP = -60;         % Controller gains (HP col)
KcB_LP = -60; KcD_LP = -60;         % Controller gains (HP col)

MDs_HP = 100; MBs_HP = 100;     % Nominal holdups (these are rather small)
MDs_LP = 100; MBs_LP = 100;     % Nominal holdups (these are rather small)
D0_HP = u0(1); B0_HP = u0(2);     % Nominal flows HP column (kmol/h)
D0_LP = u0(3); B0_LP = u0(4);     % Nominal flows LP column (kmol/h)

MB_HP = X1(NT+1);                    % Actual reboiler holdup (HP)
MD_HP = X1(2*NT);                    % Actual condenser holdup (HP)
D_HP = D0_HP + (MDs_HP-MD_HP)*KcD_HP;   % Distillate flow (HP)
B_HP = B0_HP + (MDs_HP-MB_HP)*KcD_HP;   % Bottoms flow (HP)

MB_LP = X2(Nt+1);                    % Actual reboiler holdup (LP)
MD_LP = X2(2*Nt);                    % Actual condenser holdup (LP)
D_LP = D0_LP + (MDs_LP-MD_LP)*KcD_LP;   % Distillate flow (LP)
B_LP = B0_LP + (MDs_LP-MB_LP)*KcD_LP;   % Bottoms flow (LP)

% HP/LP Heat Exchanger Calculations:
%------------------------------------
Hvap_top_HP = 31710;           % From HYSYS at 5.2bar & 0.99 MeOH
Hvap_btm_LP = 40650;           % Assume pure water at 1 bar

U=3000;                        % Overall h.t. coefficient (W/m^2 C)
DeltaT = X1(3*NT)-X2(2*Nt+1);   % (Tcond-Treb)
Q_HP = (U*A*DeltaT)*3600/1000;  % Duty in reboiler/condenser (kJ/h)

VD_HP = Q_HP/Hvap_top_HP;        % Condensed flow HP column (kmol/h)
VB_LP = Q_HP/Hvap_btm_LP;        % Boilup LP column (kmol/h)


% Reboiler duty in HP column
Hvap_btm_HP = 35440;         % From HYSYS at 5.2bar & 0.53 MeOH 
Q_HP = VB_HP*Hvap_btm_HP;     % (kJ/h)
Q_HP_M= Q_HP/(1000*3600);     % MJ/s

% Object function:
%-------------------

wr = 0.6488;                 % Relative cost (mole/MJ)
J = (D_HP + D_LP)*(1000/3600) - wr*Q_HP_M;   % Objective (mol/s)
f = -J/10;        % Scale the objective fn by 1000 

xB_HP = X1(1);
xD_HP = X1(NT);
xB_LP = X2(1);
xD_LP = X2(Nt);

Y = [LT_HP VB_HP VD_HP xB_HP xD_HP B_HP D_HP LT_LP VD_LP xB_LP xD_LP B_LP D_LP J VB_LP];
