% Script for optimization of reactor, separator and recycle process.
% The numerical description of the process is taken from Larsson et al (2003)

clc
clear all
global p;

% Column parameters
p.qF    = 1;            % Liquid fraction in column feed [-]
p.NT    = 22;           % Number of stages in distillation column [-]
p.NF    = 13;           % Feed stage in distillation column [-]
p.alpha = 2;            % Relative volatility between component A and B [-]

Vmaxs=(1500/60);        % The optimal value of vapour boilup in case II  [kmol/h]
p.Vmax  = Vmaxs;        % Making the vapour boilup a global parameter  [kmol/h]
p.xB=0.0105;            % Bottom composition (original value) [mol A/mol]
p.F=958;                % Column feed (original value) [kmol/h]
p.L=778;                % Reflux (original value) [kmol/h]
p.D=498;                % Recycle (original value) [kmol/h]
p.LD=1.6;               % Ratio L/D (original value) [-]
p.LF=0.8;               % Ratio L/F (original value) [-]
p.LV=0.61;              % Ratio L/V (original value) [-]
p.xD=0.82;              % Distillate composition (original value) [mol A/mol]
p.FF0=2;                % Ratio F/F0 (original value) [-]
p.DF0=2;                % Ratio D/F0 (original value) [-]

% CSTR parameters
F0s = (460/60);         % The optimal value of reactor feed (case I)  [kmol/h]
zF0 = 0.9;              % The optimal value of reactor feed composition [mol A/mol]

p.F0 = F0s;             % Making the reactor feed a global parameter  [kmol/h]
p.zF0 = zF0;            % Making the reactor feed a global parameter [-]
p.k1 = 0.341/60;        % Reaction rate constant [min^-1]
p.Mr=2800;              % Liquid holdup in reactor [kmol/h]
% Flags
p.OPTI=0;
p.case_I=1;             % Switching between case I and case II

if p.case_I==1          % Case I
    % Constraints
    lb=zeros(p.NT+8,1);                 % Lower bounds for the 30 variables
    ub=[ones(p.NT,1); ones(8,1)*Inf];   % Upper bounds for the 30 variables
    
    % Active constraints for case I
    % xB <= 0.0105
    ub(1)=0.0105;           
    % Mr <= 2800;
    ub(p.NT+7)=2800;        
    % F0 = fixed; 
    lb(p.NT+8)=p.F0;
    ub(p.NT+8)=p.F0;
 
else                    % Case II
    % Constraints
    lb=zeros(p.NT+8,1);                 % Lower bounds for the 30 variables
    ub=[ones(p.NT,1); ones(8,1)*Inf];   % Upper bounds for the 30 variables
    
    % Active constraints for case II
    % xB <= 0.0105
    ub(1)=0.0105;         
    % Mr <= 2800;
    ub(p.NT+7)=2800;        
    % V <= Vmax; 
    ub(p.NT+2)=p.Vmax;
end

% Initial values of the 30 variables
x0= [ones(1,p.NT)*0.5 10 15 5 5 1.1 0.5 1000 400/60]';

% fmincon options
options = optimset('TolFun',10e-6,'TolCon',10e-6,'MaxFunEvals',1e4,...
'Display','none','Algorithm','sqp','Diagnostics','off'...
);
% fmincon
[x,fval,exitflag]=fmincon(@fun,x0,[],[],[],[],lb,ub,@nlcon,options);
x0=x;

% Results
if p.case_I==1
    casename='case I: min operation cost(energy)\n';
else
    casename='case II: max production rate\n';
end
results_fmincon=sprintf(strcat(...
    casename,...
    'feed rate,             F0[kmol/h]      = %1$0.4d\n',...
    'reactor effluent,      F[kmol/h]       = %2$0.4d\n',...
    'vapor boilup,          V[kmol/h]       = %3$0.4d\n',...
    'reflux,                L[kmol/h]       = %4$0.4d\n',...
    'recycle (distilate),   D[kmol/h]       = %5$0.4d\n',...
    'recycle composition,   xD[molA/mol]    = %6$0.4d\n',...
    'bottom composition,    xB[molA/mol]    = %7$0.4d\n',...
    'reactor composition,   zF[molA/mol]    = %8$0.4d\n',...
    'reactor holdup,        Mr[kmol/h]      = %9$0.4d\n'...
    ),x(p.NT+8)*60,x(p.NT+5)*60,x(p.NT+2)*60,x(p.NT+1)*60, x(p.NT+3)*60,...
    x(p.NT), x(1), x(p.NT+6),x(p.NT+7))
% Nominal results
Vnom=x(p.NT+2)*60           % Nominal vapour boilup
F0n=-x(p.NT+8)*60           % Nominal negative feed 
% break 

%% Calculation of cost with the Brute force method
if p.case_I==1
%% Calculation of the optimal cost for different F0 values, case I.
VoptF0=zeros(1,17);                         % Preallocation of the optimal values of V
for i=1:17
    p.F0=(F0s*0.8):(11.5/60):(F0s*1.2);     % Vector of F0-values
    lb(p.NT+8)=p.F0(i);                     % Keeping F0 at a given value
    ub(p.NT+8)=p.F0(i);                     % Keeping F0 at a given value
    [x,fval,exitflag]=fmincon(@fun,x0,[],[],[],[],lb,ub,@nlcon,options);
    VoptF0(i)=x(p.NT+2)*60;                 % Vector of optimal values of V
end
 
% Calculation of the cost when F, L, D, L/D, L/F, L/V, xD, F/F0 and D/F0
% is kept constant at their optimal values with different F0 values

load yopt_nom_case1.mat                     % Matrix with optimal values of the CVs
% lb(p.NT+5)=A1(1)/60;                      % Keeping F constant
% ub(p.NT+5)=A1(1)/60;                      % Keeping F constant
% lb(p.NT+1)=A1(2)/60;                      % Keeping L constant
% ub(p.NT+1)=A1(2)/60;                      % Keeping L constant
% lb(p.NT+3)=A1(3)/60;                      % Keeping D constant
% ub(p.NT+3)=A1(3)/60;                      % Keeping D constant
% lb(p.NT)=A1(7);                           % Keeping xD constant
% ub(p.NT)=A1(7);                           % Keeping xD constant

% The remaining variables L/D, L/F, L/V, F/F0 and D/F0 was kept constant by
% adding a non-linear equality constraint in the script named nlcon.m
% because the lower and upper bounds can not be used for keeping those variables constant. 

VF0_F=zeros(1,17);                          % Vector of V-values with the different CVs kept constant
for i=1:17
    p.F0=(F0s*0.8):(11.5/60):(F0s*1.2);     % Vector of F0-values
    lb(p.NT+8)=p.F0(i);                     % Keeping F0 at a given value
    ub(p.NT+8)=p.F0(i);                     % Keeping F0 at a given value
    [x,fval,exitflag]=fmincon(@fun,x0,[],[],[],[],lb,ub,@nlcon,options);
    VF0_F(i)=x(p.NT+2)*60;                  % Vector of V-values with CV constant
end
% break 

%% Calculation of the optimal cost for different zF0 values.
p.F0 = F0s;                                 % Setting feed rate back to F0=460
lb(p.NT+8)=p.F0;                            % Keeping F0 at the given value
ub(p.NT+8)=p.F0;                            % Keeping F0 at the given value

VoptzF0=zeros(1,11);                        % Preallocation of the optimal values of V
for i=1:11
    p.zF0=0.80:0.02:1.0;                    % Vector of zF0-values
    p.zF0=p.zF0(i);                         
    [x,fval,exitflag]=fmincon(@fun,x0,[],[],[],[],lb,ub,@nlcon,options);
    VoptzF0(i)=x(p.NT+2)*60;                % Vector of optimal values of V
end


% Calculation of the cost when F, L, D, L/D, L/F, L/V, xD, F/F0 and D/F0
% is kept constant at their optimal values with different zF0 values

load yopt_nom_case1.mat                     % Matrix with optimal values of the CVs
% lb(p.NT+5)=A1(1)/60;                      % Keeping F constant
% ub(p.NT+5)=A1(1)/60;                      % Keeping F constant 
% lb(p.NT+1)=A1(2)/60;                      % Keeping L constant
% ub(p.NT+1)=A1(2)/60;                      % Keeping L constant
% lb(p.NT+3)=A1(3)/60;                      % Keeping D constant
% ub(p.NT+3)=A1(3)/60;                      % Keeping D constant
% lb(p.NT)=A1(7);                           % Keeping xD constant
% ub(p.NT)=A1(7);                           % Keeping xD constant

% The remaining variables L/D, L/F, L/V, F/F0 and D/F0 was kept constant by
% adding a non-linear equality constraint in the script named nlcon.m
% because the lower and upper bounds can not be used for keeping those variables constant. 

VzF0_DF0=zeros(1,11);                       % Vector of V-values with the different CVs kept constant
for i=1:11
    p.zF0=0.80:0.02:1.0;                    % Vector of zF0-values
    p.zF0=p.zF0(i);
    [x,fval,exitflag]=fmincon(@fun,x0,[],[],[],[],lb,ub,@nlcon,options);
    VzF0_DF0(i)=x(p.NT+2)*60;               % Vector Vs with CVs kept constant
end
% break 

%% Calculation of the cost for implementation errors in c (CVs).
p.F0 = F0s;                                 % F0=460
lb(p.NT+8)=p.F0;
ub(p.NT+8)=p.F0;
p.zF0 = zF0;                                % zF0=0.9
load yopt_nom_case1.mat                     % Matrix with optimal values of the CVs

% Implementation error in F
VieF=zeros(1,9);                            % Preallocation of Vs for implementation error in F
for i=1:9
    p.F=A1(1)*0.8/60: A1(1)*0.05/60: A1(1)*1.2/60;  % Vector of F-values with +-20% impl. error
    lb(p.NT+5)=p.F(i);                              % Keeping F at given value
    ub(p.NT+5)=p.F(i);                              % Keeping F at given value
    [x,fval,exitflag]=fmincon(@fun,x0,[],[],[],[],lb,ub,@nlcon,options);
    VieF(i)=x(p.NT+2)*60;                           % Vector of V-values for +-20% impl.error in F
end 
lb(p.NT+5)=0;
ub(p.NT+5)=Inf;
 
% Implementation error in L
VieL=zeros(1,9);                                    % Preallocation of Vs for implementation error in L
for i=1:9
    p.L=A1(2)*0.8/60: A1(2)*0.05/60: A1(2)*1.2/60;  % Vector of L-values with +-20% impl. error
    lb(p.NT+1)=p.L(i);                              % Keeping L at given value
    ub(p.NT+1)=p.L(i);                              % Keeping L at given value
    [x,fval,exitflag]=fmincon(@fun,x0,[],[],[],[],lb,ub,@nlcon,options);
    VieL(i)=x(p.NT+2)*60;                           % Vector of V-values for +-20% impl.error in L
end 
lb(p.NT+1)=0;
ub(p.NT+1)=Inf;

% Impl.error in D
VieD=zeros(1,9);                                    % Preallocation of Vs for implementation error in D
for i=1:9
    p.D=A1(3)*0.8/60: A1(3)*0.05/60: A1(3)*1.2/60;  % Vector of D-values with +-20% impl. error
    lb(p.NT+3)=p.D(i);                              % Keeping D at given value
    ub(p.NT+3)=p.D(i);                              % Keeping D at given value
    [x,fval,exitflag]=fmincon(@fun,x0,[],[],[],[],lb,ub,@nlcon,options);
    VieD(i)=x(p.NT+2)*60;                           % Vector of V-values for +-20% impl.error in D
end 
lb(p.NT+3)=0;
ub(p.NT+3)=Inf;

% Impl.error in xD
ViexD=zeros(1,9);                                   % Preallocation of Vs for implementation error in xD
for i=1:9
    p.xD=A1(7)*0.8: A1(7)*0.05: A1(7)*1.2;          % Vector of xD-values with +-20% impl. error
    lb(p.NT)=p.xD(i);                               % Keeping xD at given value
    ub(p.NT)=p.xD(i);                               % Keeping xD at given value
    [x,fval,exitflag]=fmincon(@fun,x0,[],[],[],[],lb,ub,@nlcon,options);
    ViexD(i)=x(p.NT+2)*60;                          % Vector of V-values for +-20% impl.error in xD
end 
lb(p.NT)=0;
ub(p.NT)=1;

% The operational cost for +-20% implementation errors in the remaining
% variables (L/D, L/F, L/V, F/F0 and D/F0) was found by adding equality
% constraints in the script nlcon.m, and were then typed in below for
% simplicity. Again, this is because the lower and upper bounds could not
% be used to keep ratios constant. 

% Impl.error in L/D 
Vie_LD=[1278.9 1277.6 1276.6 1276 1275.7 1276 1276.7 1278.1 1280];

% Impl.error in L/F
Vie_LF=[1281.1 1279.0 1277.3 1276.2 1275.7 1276.2 1277.9 1281.1 1286.2];

% Impl.error in L/V
Vie_LV=[1286.7 1283 1279.6 1276.9 1275.7 1277.7 1286.6 1311 1369.1];

% Impl.error in F/F0
Vie_FF0=[7117.3 1460.3 1299.8 1279.1 1275.7 1277.4 1280.9 1285.1 1289.4];

% Impl.error in D/F0
Vie_DF0=[1303.3 1286.7 1279.5 1276.5 1275.7 1276.3 1277.5 1279.2 1281.2];
% break 

%% Optimal cost for implementation errors in Mr
p.F0=F0s;
p.zF0=zF0;

VoptMr=zeros(1,11);                                 % Preallocation of the optimal values of V
for i=1:11
    p.Mr=2300:50:2800;                              % Vector of Mr-values (infeasible with Mr>2800 kmol/h)
    lb(p.NT+7)=p.Mr(i);                             % Keeping Mr at given value
    ub(p.NT+7)=p.Mr(i);                             % Keeping Mr at given value
    [x,fval,exitflag]=fmincon(@fun,x0,[],[],[],[],lb,ub,@nlcon,options);
    VoptMr(i)=x(p.NT+2)*60;                         % Vector of optimal V-values with impl.error in Mr
end 

% Calculation of the cost when F, L, D, L/D, L/F, L/V, xD, F/F0 and D/F0
% is kept constant at their optimal values with different Mr values

load yopt_nom_case1.mat                             % Matrix with optimal values of the CVs for case I
% lb(p.NT+5)=A1(1)/60;                              % Keeping F constant
% ub(p.NT+5)=A1(1)/60;                              % Keeping F constant 
% lb(p.NT+1)=A1(2)/60;                              % Keeping L constant
% ub(p.NT+1)=A1(2)/60;                              % Keeping L constant
% lb(p.NT+3)=A1(3)/60;                              % Keeping D constant
% ub(p.NT+3)=A1(3)/60;                              % Keeping D constant
% lb(p.NT)=A1(7);                                   % Keeping xD constant
% ub(p.NT)=A1(7);                                   % Keeping xD constant

% The remaining variables L/D, L/F, L/V, F/F0 and D/F0 was kept constant by
% adding a non-linear equality constraint in the script named nlcon.m
% because the lower and upper bounds can not be used for keeping those variables constant. 

VMr_F=zeros(1,11);                                  % Preallocation of V-values with CVs kept constant
for i=1:11
    p.Mr=2300:50:2800;                              % Vector of Mr-values
    lb(p.NT+7)=p.Mr(i);                             % Keeping Mr at given value
    ub(p.NT+7)=p.Mr(i);                             % Keeping Mr at given value
    [x,fval,exitflag]=fmincon(@fun,x0,[],[],[],[],lb,ub,@nlcon,options);
    VMr_F(i)=x(p.NT+2)*60;                          % Vector of V-values with impl. error in CVs
end 
lb(p.NT+7)=0;
ub(p.NT+7)=2800;

%%
% Optimal cost for implementation errors in xB
p.F0=F0s;
p.zF0=zF0;

VoptxB=zeros(1,11);                                 % Preallocation of optimal values of V with impl. error in xB
for i=1:11
    p.xB=0.009:0.00015:0.0105;                      % Vector of xB-values (infeasible with xB>0.0105 mol A/mol)
    lb(1)=p.xB(i);                                  % Keeping xB at given value
    ub(1)=p.xB(i);                                  % Keeping xB at given value
    [x,fval,exitflag]=fmincon(@fun,x0,[],[],[],[],lb,ub,@nlcon,options);
    VoptxB(i)=x(p.NT+2)*60;                         % Vector of optimal V-values with impl. error in xB
end


% Calculation of the cost when F, L, D, L/D, L/F, L/V, xD, F/F0 and D/F0
% is kept constant at their optimal values with different Mr values

load yopt_nom_case1.mat                             % Matrix with optimal values of the CVs
% lb(p.NT+5)=A1(1)/60;                              % Keeping F constant
% ub(p.NT+5)=A1(1)/60;                              % Keeping F constant 
% lb(p.NT+1)=A1(2)/60;                              % Keeping L constant
% ub(p.NT+1)=A1(2)/60;                              % Keeping L constant
% lb(p.NT+3)=A1(3)/60;                              % Keeping D constant
% ub(p.NT+3)=A1(3)/60;                              % Keeping D constant
% lb(p.NT)=A1(7);                                   % Keeping xD constant
% ub(p.NT)=A1(7);                                   % Keeping xD constant

% The remaining variables L/D, L/F, L/V, F/F0 and D/F0 was kept constant by
% adding a non-linear equality constraint in the script named nlcon.m
% because the lower and upper bounds can not be used for keeping those variables constant. 

VxB_F=zeros(1,11);                                  % Preallocation of V-values with CVs kept constant
for i=1:11
    p.xB=0.009:0.00015:0.0105;                      % Vector of xB-values
    lb(1)=p.xB(i);                                  % Keeping xB at given value
    ub(1)=p.xB(i);                                  % Keeping xB at given value      
    [x,fval,exitflag]=fmincon(@fun,x0,[],[],[],[],lb,ub,@nlcon,options);
    VxB_F(i)=x(p.NT+2)*60;                          % Vector of V-values with CVs kept constant
end
lb(1)=0;
ub(1)=1;

else
    
%% Case II
% Calculation of the optimal cost for different Vmax values.
F0oVm=zeros(1,13);                              % Preallocation of the optimal values of F0
for i=1:13
    p.Vmax=1200/60:50/60:1800/60;               % Vector of Vmax-values
    lb(p.NT+2)=p.Vmax(i);                       % Keeping Vmax at given value
    ub(p.NT+2)=p.Vmax(i);                       % Keeping Vmax at given value
    [x,fval,exitflag]=fmincon(@fun,x0,[],[],[],[],lb,ub,@nlcon,options);
    F0oVm(i)=-x(p.NT+8)*60;                     % Vector of optimval F0-values
end
 
% Calculation of the cost when F, L, D, L/D, L/F, L/V, xD, F/F0 and D/F0
% is kept constant at their optimal values with different Vmax values

load yopt_nom_case2.mat                         % Matrix with optimal values of the CVs
% lb(p.NT+5)=A4(1)/60;                          % Keeping F constant
% ub(p.NT+5)=A4(1)/60;                          % Keeping F constant 
% lb(p.NT+1)=A4(2)/60;                          % Keeping L constant
% ub(p.NT+1)=A4(2)/60;                          % Keeping L constant
% lb(p.NT+3)=A4(3)/60;                          % Keeping D constant
% ub(p.NT+3)=A4(3)/60;                          % Keeping D constant
% lb(p.NT)=A4(7);                               % Keeping xD constant
% ub(p.NT)=A4(7);                               % Keeping xD constant

% The remaining variables L/D, L/F, L/V, F/F0 and D/F0 was kept constant by
% adding a non-linear equality constraint in the script named nlcon.m
% because the lower and upper bounds can not be used for keeping those variables constant. 

F0Vm_CV=zeros(1,13);                            % Preallocation of F0-values
for i=1:13
    p.Vmax=1200/60:50/60:1800/60;               % Vector of Vmax-values
    lb(p.NT+2)=p.Vmax(i);                       % Keeping Vmax at given value
    ub(p.NT+2)=p.Vmax(i);                       % Keeping Vmax at given value
    [x,fval,exitflag]=fmincon(@fun,x0,[],[],[],[],lb,ub,@nlcon,options);
    F0Vm_CV(i)=-x(p.NT+8)*60;                   % Vector of F0-values with CVs constant
end
% break 

%% Calculation of the optimal cost for different zF0 values.
p.Vmax = Vmaxs;                                 % Vmax=1500
lb(p.NT+2)=p.Vmax;                              % Keeping Vmax at given value
ub(p.NT+2)=p.Vmax;                              % Keeping Vmax at given value


F0ozF0=zeros(1,11);                             % Preallocation of optimal values of F0
for i=1:11
    p.zF0=0.80:0.02:1.0;                        % Vector of zF0-values
    p.zF0=p.zF0(i);
    [x,fval,exitflag]=fmincon(@fun,x0,[],[],[],[],lb,ub,@nlcon,options);
    F0ozF0(i)=-x(p.NT+8)*60;                    % Vector of optimal F0-values
end

% Calculation of the cost when F, L, D, L/D, L/F, L/V, xD, F/F0 and D/F0
% is kept constant at their optimal values with different zF0 values

load yopt_nom_case2.mat                         % Matrix with optimal values of the CVs
% lb(p.NT+5)=A4(1)/60;                          % Keeping F constant
% ub(p.NT+5)=A4(1)/60;                          % Keeping F constant 
% lb(p.NT+1)=A4(2)/60;                          % Keeping L constant
% ub(p.NT+1)=A4(2)/60;                          % Keeping L constant
% lb(p.NT+3)=A4(3)/60;                          % Keeping D constant
% ub(p.NT+3)=A4(3)/60;                          % Keeping D constant
% lb(p.NT)=A4(7);                               % Keeping xD constant
% ub(p.NT)=A4(7);                               % Keeping xD constant

F0zF0_CV=zeros(1,11);                               % Preallocation of F0-values
for i=1:11
    p.zF0=0.80:0.02:1.0;                            % Vector with zF0-values
    p.zF0=p.zF0(i);
    [x,fval,exitflag]=fmincon(@fun,x0,[],[],[],[],lb,ub,@nlcon,options);
    F0zF0_CV(i)=-x(p.NT+8)*60;                      % Vector with F0-values with CVs constant
end
% break 
%% Calculation of the cost for implementation errors in c (CVs).
p.zF0 = zF0;                                        % zF0=0.9
p.Vmax=Vmaxs;                                       % Vmax=1500
ub(p.NT+2)=p.Vmax;
load yopt_nom_case2.mat                             % Matrix with optimal values of the CVs for case 2

% Impl.error in F
F0ieF=zeros(1,9);                                   % Preallocation of F0-vector
for i=1:9
    p.F=A4(1)*0.8/60: A4(1)*0.05/60: A4(1)*1.2/60;  % Vector of F-values with +-20% implerror
    lb(p.NT+5)=p.F(i);                              % Keeping F at given value
    ub(p.NT+5)=p.F(i);                              % Keeping F at given value
    [x,fval,exitflag]=fmincon(@fun,x0,[],[],[],[],lb,ub,@nlcon,options);
    F0ieF(i)=-x(p.NT+8)*60;                         % Vector of F0-values for impl.error in F
end 
lb(p.NT+5)=0;
ub(p.NT+5)=Inf;
 
% Impl.error in L
F0ieL=zeros(1,9);                                   % Preallocation of F0-vector
for i=1:9
    p.L=A4(2)*0.8/60: A4(2)*0.05/60: A4(2)*1.2/60;  % Vector of L-values with +-20% implerror
    lb(p.NT+1)=p.L(i);                              % Keeping L at given value
    ub(p.NT+1)=p.L(i);                              % Keeping L at given value
    [x,fval,exitflag]=fmincon(@fun,x0,[],[],[],[],lb,ub,@nlcon,options);
    F0ieL(i)=-x(p.NT+8)*60;                         % Vector of F0-values for impl.error in L
end 
lb(p.NT+1)=0;
ub(p.NT+1)=Inf;

% Impl.error in D
F0ieD=zeros(1,9);                                   % Preallocation of F0-vector
for i=1:9
    p.D=A4(3)*0.8/60: A4(3)*0.05/60: A4(3)*1.2/60;  % Vector of D-values with +-20% implerror
    lb(p.NT+3)=p.D(i);                              % Keeping D at given value
    ub(p.NT+3)=p.D(i);                              % Keeping D at given value
    [x,fval,exitflag]=fmincon(@fun,x0,[],[],[],[],lb,ub,@nlcon,options);
    F0ieD(i)=-x(p.NT+8)*60;                         % Vector of F0-values for impl.error in D
end 
lb(p.NT+3)=0;
ub(p.NT+3)=Inf;

% Impl.error in xD
F0iexD=zeros(1,9);                                  % Preallocation of F0-vector
for i=1:9
    p.xD=A4(7)*0.8: A4(7)*0.05: A4(7)*1.2;          % Vector of xD-values with +-20% implerror
    lb(p.NT)=p.xD(i);                               % Keeping xD at given value
    ub(p.NT)=p.xD(i);                               % Keeping xD at given value
    [x,fval,exitflag]=fmincon(@fun,x0,[],[],[],[],lb,ub,@nlcon,options);
    F0iexD(i)=-x(p.NT+8)*60;                        % Vector of F0-values for impl.error in xD
end 
lb(p.NT)=0;
ub(p.NT)=1;

% The operational cost for +-20% implementation errors in the remaining
% variables (L/D, L/F, L/V, F/F0 and D/F0) was found by adding equality
% constraints in the script nlcon.m, and were then typed in below for
% simplicity. Again, this is because the lower and upper bounds could not
% be used to keep ratios constant. 

% Impl.error in L/D
F0ie_LD=[-497.2457 -497.5115 -497.7107 -497.8362 -497.8799 -497.8332 -497.6872 -497.4323 -497.0590];

% Impl.error in L/F
F0ie_LF=[-496.89 -497.2817 -497.5937 -497.8027 -497.8799 -497.7899 -497.4908 -496.9346 -496.0690];

% Impl.error in L/V
F0ie_LV=[-495.8496 -496.5349 -497.1659 -497.6631 -497.8799 -497.5409 -496.1306 -492.7389 -486.0753];

% Impl.error in F/F0
F0ie_FF0=[-466.3790 -485.9074 -494.4255 -497.2917 -497.8799 -497.5584 -496.8714 -496.0507 -495.2004];

% Impl.error in D/F0
F0ie_DF0=[-493.2664 -495.8115 -497.1334 -497.7255 -497.8799 -497.7692 -497.4975 -497.1291 -496.7046];
% break

%% Optimal cost for implementation errors in Mr
p.zF0=zF0;
p.Vmax=Vmaxs;
ub(p.NT+2)=p.Vmax;

F0oMr=zeros(1,11);                              % Preallocation of the optimal values of F0
for i=1:11
    p.Mr=2300:50:2800;                          % Vector of Mr-values (infeasible with MR>2800 kmol/h)
    lb(p.NT+7)=p.Mr(i);                         % Keeping Mr at given value
    ub(p.NT+7)=p.Mr(i);                         % Keeping Mr at given value
    [x,fval,exitflag]=fmincon(@fun,x0,[],[],[],[],lb,ub,@nlcon,options);
    F0oMr(i)=-x(p.NT+8)*60;                     % Vector of optimal F0-values with impl. error in Mr
end

% Calculation of the cost when F, L, D, L/D, L/F, L/V, xD, F/F0 and D/F0
% is kept constant at their optimal values with different Mr values

load yopt_nom_case2.mat                         % Matrix with optimal values of the CVs
% lb(p.NT+5)=A4(1)/60;                          % Keeping F constant
% ub(p.NT+5)=A4(1)/60;                          % Keeping F constant 
% lb(p.NT+1)=A4(2)/60;                          % Keeping L constant
% ub(p.NT+1)=A4(2)/60;                          % Keeping L constant
% lb(p.NT+3)=A4(3)/60;                          % Keeping D constant
% ub(p.NT+3)=A4(3)/60;                          % Keeping D constant
% lb(p.NT)=A4(7);                               % Keeping xD constant
% ub(p.NT)=A4(7);                               % Keeping xD constant

% The remaining variables L/D, L/F, L/V, F/F0 and D/F0 was kept constant by
% adding a non-linear equality constraint in the script named nlcon.m
% because the lower and upper bounds can not be used for keeping those variables constant. 

F0Mr_F=zeros(1,11);                             % Preallocation of F0-values
for i=1:11  
    p.Mr=2300:50:2800;                          % Vector of Mr-values
    lb(p.NT+7)=p.Mr(i);                         % Keeping Mr at given value
    ub(p.NT+7)=p.Mr(i);                         % Keeping Mr at given value
    [x,fval,exitflag]=fmincon(@fun,x0,[],[],[],[],lb,ub,@nlcon,options);
    F0Mr_F(i)=-x(p.NT+8)*60;                    % Vector of F0-values with CVs kept constant
end 
lb(p.NT+7)=0;
ub(p.NT+7)=2800;
% break 
%% Optimal cost for implementation errors in xB
ub(p.NT+7)=2800;        % Mr=2800
p.zF0=zF0;              % zF0=0.9

F0oxB=zeros(1,11);                              % Preallocation of optimal F0-values
for i=1:11  
    p.xB=0.009:0.00015:0.0105;                  % Vector of xB-values (infeasible with xB>0.0105 mol A/mol)
    lb(1)=p.xB(i);                              % Keeping xB at given value
    ub(1)=p.xB(i);                              % Keeping xB at given value
    [x,fval,exitflag]=fmincon(@fun,x0,[],[],[],[],lb,ub,@nlcon,options);
    F0oxB(i)=-x(p.NT+8)*60;                     % Vector of optimal F0-values with impl. error in xB
end

% Calculation of the cost when F, L, D, L/D, L/F, L/V, xD, F/F0 and D/F0
% is kept constant at their optimal values with different Mr values

load yopt_nom_case2.mat                         % Matrix with optimal values of the CVs
% lb(p.NT+5)=A4(1)/60;                          % Keeping F constant
% ub(p.NT+5)=A4(1)/60;                          % Keeping F constant 
% lb(p.NT+1)=A4(2)/60;                          % Keeping L constant
% ub(p.NT+1)=A4(2)/60;                          % Keeping L constant
% lb(p.NT+3)=A4(3)/60;                          % Keeping D constant
% ub(p.NT+3)=A4(3)/60;                          % Keeping D constant
% lb(p.NT)=A4(7);                               % Keeping xD constant
% ub(p.NT)=A4(7);                               % Keeping xD constant

F0xB_DF0=zeros(1,11);                           % Preallocation of F0-values
for i=1:11
    p.xB=0.009:0.00015:0.0105;                  % Vector of xB-values
    lb(1)=p.xB(i);                              % Keeping xB at given value
    ub(1)=p.xB(i);                              % Keeping xB at given value
    [x,fval,exitflag]=fmincon(@fun,x0,[],[],[],[],lb,ub,@nlcon,options);
    F0xB_DF0(i)=-x(p.NT+8)*60;                  % Vector of F0-values with CVs kept constant
end
lb(1)=0;
ub(1)=1;
% break
end 