%************************************************************** 
%          RUN OPTIMISATION OF MULTIEFFECT COLUMNS 
%
%   - Runs the program optimisation.m
%   - Finds the nominal optimum (no disturbances) and the optimum
%     J(d,u) at disturbances d (disturbances in F, zF and qF)
%
% H.K.Engelien 01/07/2001
%**************************************************************


global Results Xfinal uopt

format short e

NT = 16; Nt = 26;

load Results_steady_state
u0 = [D_HP B_HP D_LP B_LP];  % Nominal values for control loops
Y1=[];
X1=[];
Ainit=1000;
%uxopt0 = [LT_HP;VB_HP;Ainit;Xss(1:3*NT);LT_LP;VD_LP;Xss(3*NT+1:3*NT+3*Nt)];

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

load Results_optimisation Y X
uxopt0=[Y(1,14); Y(1,30); Ainit; X(2:3*NT+1,1);Y(1,15);Y(1,31);X(3*NT+1+1:3*NT+3*Nt+1,1)];

% No disturbances (nominal operating point)
d=[zF_HP; qF_HP; F_HP];
optimisation(d,uxopt0,u0);
Y1=[Y1;0 Results];
Xfinal = [0;Xfinal];
X1 = [ X1 Xfinal];
i=0

'Results at nominal optimum (no disturbances):'
LT_HPopt=uopt(1), VB_HPopt=uopt(2), VD_HPopt=Results(28)
LT_LPopt=uopt(4), VB_LPopt=Results(29), VD_LPopt=uopt(5)
D_HPopt=Results(24),B_HPopt=Results(12),D_LPopt=Results(25)
B_LPopt=Results(26),xD_HPopt=Results(21), xD_LPopt=Results(22)
P_HPopt=Results(4),P_LPopt=Results(20), 
Q_LPopt=Results(19), Aopt = uopt(3), Jopt = Results(2)
Xfinal_opt=Xfinal;

uxopt0=[Y(2,14); Y(2,30); Ainit; X(2:3*NT+1,2);Y(2,15);Y(2,31);X(3*NT+1+1:3*NT+3*Nt+1,2)];

% Disturbance in Feedrate: F + 20%
d1=F_HP+0.2*F_HP; 
d=[zF_HP; qF_HP; d1];
optimisation(d,uxopt0,u0);
Y1=[Y1;1 Results];
Xfinal = [1;Xfinal];
X1 = [X1 Xfinal];
i=1

uxopt0=[Y(3,14); Y(3,30); Ainit; X(2:3*NT+1,3);Y(3,15);Y(3,31);X(3*NT+1+1:3*NT+3*Nt+1,3)];

% Disturbance in Feedrate: F - 20%
d2=F_HP-0.2*F_HP;       % Decrease feedrate by 10%
d=[zF_HP; qF_HP; d2];
optimisation(d,uxopt0,u0);
Y1=[Y1;2 Results];
Xfinal = [2;Xfinal];
X1=[X1 Xfinal];
i=2



% Results from optimisation
Y1                       % Results for each disturbances
%save Results_optimisation Y  X LT_HPopt VB_HPopt VD_HPopt LT_LPopt VB_LPopt VD_LPopt Xfinal_opt D_HPopt B_HPopt D_LPopt B_LPopt
