function p=optimisation(d,uxopt0,u0)
%**************************************************************
% OPTIMISATION

% Main program for self-optimising control 
% of multi-effect distillation columns
% H.K. Engelien, 02/04/2001

% The program has the following subroutines:

%   cola_lv_fun     - the object function
%   cola_lv_nonlcon - generates the constraints 
%   colamod1        - the model for the HP column
%   colamod2        - the model for the LP column
%   Xss_multieffect - initial values for colamod1 and colamod2
%
%**************************************************************
global Y P_HP P_LP  T Results Xfinal uopt S

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

% Feed specifications - nominal:
%-------------------------------
zF_HP = d(1);        % Feed composition light component p(1)
qF_HP = d(2);;       % Feed liquid fraction p(2)
F_HP =  d(3);        % Feedrate (kmol/h) p(3)

% Prices for object function:
%----------------------------
PD =  5.13;              % Distillate price p(4)
PV =  0.5;               % Boilup price p(5)

% Max/min values for inequality constraints:
%--------------------------------------------

% HP Column
xD_HPmin = 0.99;     % Min light component in HP column, p(6)
P_HPmax = 10e5;      % Max pressure HP column (Pa), p(7)
VB_HPmin = 0;        % Min reboil flow (kmol/h), p(8) 
LT_HPmin = 0;        % Min reflux flow (kmol/h), p(9)
VD_HPmin = 0;        % Min. flow condensated (kmol/h) p(10)

% LP Column
xD_LPmin = 0.99;     % Min light component in LP column, p(11)
P_LPmin = 1.01e5;    % Min pressure LP column (Pa), p(12)
VB_LPmin = 0;        % Min reboil flow (kmol/h), p(13) 
LT_LPmin = 0;        % Min reflux flow  (kmol/h), p(14)
VD_LPmin = 0;        % Min. flow condensated (kmol/h) p(15)

pnom = [zF_HP; qF_HP; F_HP; PD; PV; xD_HPmin; P_HPmax; VB_HPmin;...
        LT_HPmin; VD_HPmin; xD_LPmin; P_LPmin; VB_LPmin;... 
        LT_LPmin; VD_LPmin];

p = pnom;

%***************************************************************
%                   Steady state optimisation
%***************************************************************

% Options
options        =optimset;
options.Display='iter';     % Displays the interations

% Bounds
LB=0;%[];
UB=[];


% Optimisation (run optimisation twice to ensure convergence) 

[uxopt,fval,exitflag,output,lambda,grad]=fmincon('cola_lv_fun',uxopt0,[],[],[],[],LB,UB,'cola_lv_nonlcon',options,p,u0);
[uxopt,fval,exitflag,output,lambda,grad]=fmincon('cola_lv_fun',uxopt,[],[],[],[],LB,UB,'cola_lv_nonlcon',options,p,u0);

exitflag         % If algorithm has converted then exitflag >0

%-------------------------------------------------------------------
%                           Results
%-------------------------------------------------------------------

% The states [x M T] in the columns :
xopt_HP=  [uxopt(4:51)];       % HP column
xopt_LP = [uxopt(54:131)];     % LP column

popt=p; 
yopt=Y';         % [LT VB VD xB xD B D J] for both columns 

% Evaluating contraints at optimum
[C, Ceq] = cola_lv_nonlcon(uxopt,p,u0);        % Scaled C and Ceq

%A = [5e5 1000 1000 1000 1000 1000 1000];  %0.01 0.01 5e5 10
A = [ 1 1 1 1 1 1 1];
S1 = diag(A);
C_unscaled = S1*C;

S2 = eye(NT*3+Nt*3+4)*100;
Ceq_unscaled = S2*Ceq;

% Value of objective function
Jopt  = -fval*10;                     % Objective (mol/s) (scale back)

uopt  = [uxopt(1:3);uxopt(52:53)];      % Optimum [LT_HP VB_HP A LT_LP VD_LP]

% Concentration in top and btm of columns
xD = Y(1,[5 11]);                      % [HP LP]
xB = Y(1,[4 10]);                      % [HP LP]

% Column flowrates (kmol/h)
D_HP = Y(7);
B_HP = Y(6);
D_LP = Y(13);
B_LP = Y(12);

%Column pressures (Pa), temperatures (K), and duty (kJ/h)
P_btm_HP = P_HP(1);
P_top_HP = P_HP(NT);
P_btm_LP = P_LP(1);
P_top_LP = P_LP(Nt);
T_top_HP = T(1);
T_btm_LP = T(2);

Hvap_btm_HP = 35440;         % From HYSYS at 5.2bar & 0.53 MeOH 
QB_HP = uopt(2)*Hvap_btm_HP;  % (kJ/h)

VD_HP = Y(3);
VB_LP = Y(15);
VD_LP = Y(9);

Hvap_btm_LP = 40650;          % Assume pure water at 1 bar
QB_LP = uopt(5)*Hvap_btm_LP;   % (kJ/h)



Xfinal = [xopt_HP;xopt_LP];


Results =[exitflag Jopt QB_HP P_top_HP (P_btm_HP-P_top_HP) xB(1) xB(2) T_btm_LP...
          xopt_LP(2*Nt+2) xopt_LP(2*Nt+4) xopt_LP(2*Nt+6) B_HP...
          uopt(1) uopt(4) QB_HP/F_HP QB_HP/uopt(1) QB_HP/uopt(4)... 
          B_HP/F_HP  QB_LP P_top_LP xD(1) xD(2) uopt(3) D_HP D_LP B_LP uopt(3) VD_HP VB_LP VD_LP];

