%***************************************************************
%           CALCULATING THE LOSS IN THE OBJECT FUNCTION
%
%   - calculates the loss, L at different disturbances, while
%     keeping the controlled variable, c constant
%   - the following programs are used:
%           broyden1    - solves f(x)=0 with Broydens method
%           multieffect - combines the HP and LP column
%           colamod1    - the HP column model
%           colamod2    - the LP column model
%   - inputs:  Results_optimisation
%   - outputs: Results_loss
%
% H.K. Engelien 28/06/2001
%***************************************************************

NT = 16;                            % No. stages in HP column
Nt = 26;                            % No. stages in LP column

wr = 0.6488;                 % Relative cost (mole/MJ)

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 ];

r=size(cs0);

% Holdups in top and btm at optimum
MB_HP=Xfinal_opt(NT+1+1);
MD_HP=Xfinal_opt(2*NT+1);
MB_LP=Xfinal_opt(3*NT+Nt+1+1);
MD_LP=Xfinal_opt(3*NT+2*Nt+1);
m = [MB_HP MD_HP MB_LP MD_LP];

J = [];
Loss = [];
Exitflag=[];
Jopt = [Y(:,3)];       % Optimal values of J for each disturbance

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

% Reboiler duty in HP column
Hvap_btm_HP = 35440;         % From HYSYS at 5.2bar & 0.53 MeOH 


for i= 1:r(2)
    Jd = [];
    E=[];
    cs= cs0(i);                 % Constant controlled variable
    n=i
    
    % Initial values
    u = [D_HPopt;B_HPopt;LT_HPopt;VB_HPopt;D_LPopt;B_LPopt;VD_LPopt;LT_LPopt]; %VD_HPopt;
    ux0= [Xfinal_opt(2:3*NT+3*Nt+1);u];       % Initial values for broyden routine
    B0 = [];                    % Jacobi matrix

    % Nominal operating     
    d=[F_HP; zF_HP; qF_HP];
    [ux,B] = broyden1('multieffect_mod',ux0,B0,d,n,cs,m);
    QB_HP = ux(3*NT+3*Nt+4)*Hvap_btm_HP;    % (kJ/h)
    J_d = (ux(3*NT+3*Nt+1) + ux(3*NT+3*Nt+5))*(1000/3600)-wr*QB_HP/(1000*3600);    % Objective (mol/s)
    Jd=[Jd;J_d];
   
    %Disturbance d1 : F + 20%
    d1=F_HP+F_HP*0.2;          
    d=[d1; zF_HP; qF_HP];
    options=optimset; options.Display='iter';options.MaxFunEvals=150000;options.MaxIter=3000;
    [ux,B,exitflag] = fsolve('multieffect_mod',ux0,options,d,n,cs,m);
    %[ux,B] = broyden1('multieffect_mod',ux,B0,d,n,cs,m);
    
    QB_HP = ux(3*NT+3*Nt+4)*Hvap_btm_HP;    % (kJ/h)
    J_d = (ux(3*NT+3*Nt+1) + ux(3*NT+3*Nt+5))*(1000/3600)-wr*QB_HP/(1000*3600);    % Objective (mol/s)
    Jd=[Jd;J_d];
    E=[exitflag];
    
    %Disturbance d2 : F - 20%
    d2=F_HP-0.2*F_HP;          
    d=[d2; zF_HP; qF_HP];
    options=optimset; options.Display='iter';options.MaxFunEvals=150000;options.MaxIter=3000;
    [ux,B] = fsolve('multieffect_mod',ux0,options,d,n,cs,m);
    %[ux,B] = broyden1('multieffect_mod',ux,B0,d,n,cs,m);
    
    QB_HP = ux(3*NT+3*Nt+4)*Hvap_btm_HP;    % (kJ/h)
    J_d = (ux(3*NT+3*Nt+1) + ux(3*NT+3*Nt+5))*(1000/3600)-wr*QB_HP/(1000*3600);    % Objective (mol/s)
    Jd=[Jd;J_d];
    E=[E;exitflag];
    
    J = [J Jd];             % Objective function
    Exitflag=[Exitflag E];
    
    L = Jopt - Jd;
    Loss = [Loss L];
end
    % Evaluate loss, L at all disturbances and c's
    J
    Exitflag
    Loss
   

  