% 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=0;             % 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 

%% Plotting loss with the Null Space Method and Exact Local Method
% In order for the plots to be made correctly, the equality constraints
% with the septoints found above need to be added to nlcon.m for each of
% the plots.
if p.case_I==1
    %% Exact Local Method, Case I
    Vel=zeros(1,17);                        % Preallocation of V-values
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 given value
    ub(p.NT+8)=p.F0(i);                     % Keeping F0 at given value
    [x,fval,exitflag]=fmincon(@fun,x0,[],[],[],[],lb,ub,@nlcon,options);
    Vel(i)=x(p.NT+2)*60;                    % Cost with the exact local method
end
load VoptF0                                 % Optimal cost for disturbance in F0, saved earlier.
loss_el=((Vel-VoptF0)./VoptF0)*100;         % Calculating the loss
F0=368:11.5:552;                            % F0 values for the x-axis
% Making the plot of the loss versus F0
figure (2)
xx=F0;
yy=spline(F0,loss_el,xx);
plot(xx,yy,'b','LineWidth',2.5)
xlim([370 550])
ylim([0 1])
set(gca,'XTick',400:50:550);
set(gca,'YTick',0:0.2:1);
xlabel('F_0 [kmol/h]','FontSize',12,'FontWeight','bold')
ylabel('Loss [%]','FontSize',12,'FontWeight','bold')
title('Exact Local Method - Case 1','FontSize',12,'FontWeight','bold')

else
    %% Exact Local Case II
    F0el=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);
    F0el(i)=-x(p.NT+8)*60;                  % Cost with the exact local method
end
load F0oVm                                  % Optimal cost for disturbance in Vmax, saved earlier.
loss_el=((F0el-F0oVm)./abs(F0oVm))*100;     % Calculating the loss
Vm=1200:50:1800;                            % Vmax-values for the x-axis
% Making the plot of the loss versus Vmax
figure (4)
xx=Vm;
yy=spline(Vm,loss_el,xx);
plot(xx,yy,'b','LineWidth',2.5)
xlim([1200 1800])
ylim([0 0.2])
set(gca,'XTick',1200:200:1800);
set(gca,'YTick',0:0.05:0.2);
xlabel('V_m_a_x [kmol/h]','FontSize',12,'FontWeight','bold')
ylabel('Loss [%]','FontSize',12,'FontWeight','bold')
title('Exact Local Method - Case 2','FontSize',12,'FontWeight','bold')   
end
