%************************************************************** 
%          CALCULATION OF LOSS
%
%   - Runs the program optimisation.m
%   - Finds the loss when no disturbances (=0) and the loss
%     at disturbances d1 and d2(F+20% and F-20%)
%
% 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

% Specify set points for the controlled variables, cs
% (the nominal values are the values of c at the optimum)

Q_HP = Y(1,4);         % Heat duty HP column
P_top_HP  = Y(1,5);    % Pressure in top of HP column
DeltaP_HP = Y(1,6);    % Pressure drop HP column
xB_HP = Y(1,7);        % Btm. conc. of MeOH, HP column
xB_LP = Y(1,8);        % Btm. conc. of water, LP column
TB_LP = Y(1,9);        % Temperature in btm LP column
T2_LP = Y(1,10);       % Temperature on tray 2, LP column
T4_LP = Y(1,11);       % Temperature on tray 4, LP column
T6_LP = Y(1,12);       % Temperature on tray 6, LP column
B_HP  = Y(1,13);       % Btm flow in HP column
LT_HP = Y(1,14);       % Reflux in HP column
LT_LP = Y(1,15);       % Reflux in LP column
Q_HP_F = Y(1,16);      % Ratio btw duty and feed in HP column
Q_HP_LT_HP = Y(1,17);  % Ratio btw duty and reflux in HP column
Q_HP_LT_LP = Y(1,18);  % Ratio btw duty (HP) and reflux in LP col
B_HP_F = Y(1,19);      % Ratio btw btm flow and feed in HP column
TB_HP  = X(1+2*NT+1,1); % Temperature in btm. of HP column

cs0 = [Q_HP P_top_HP DeltaP_HP  xB_LP xB_HP TB_LP T2_LP T4_LP...
      T6_LP TB_HP B_HP LT_LP LT_HP Q_HP_F Q_HP_LT_LP Q_HP_LT_HP B_HP_F ];
  
Jopt=[Y(:,3)];
J=[]; 
Exitflag=[];
L=[];
for i=17:17
    cs= cs0(i);                 % Constant controlled variable
    n=i
    Jd=[];
    E=[];
    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,cs,n);
    E=[E; Results(1)];
    J_d=Results(2);
    Jd=[Jd;J_d];
    
 
    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,cs,n);
    E=[E; Results(1)];
    J_d=Results(2);
    Jd=[Jd;J_d];
     
    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,cs,n);
    E=[E; Results(1)];
    J_d=Results(2);
    Jd=[Jd;J_d];
    Loss=Jopt-Jd;
    
    J=[J Jd];
    Exitflag=[Exitflag E];
    L=[L Loss];
end
J
Exitflag
Loss
   
