% 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 

%% Exact local method 
if p.case_I==1
% Case1
p.F0 = F0s;
lb(p.NT+8)=p.F0;
ub(p.NT+8)=p.F0;
p.zF0 = zF0;

% Running the optimization
[x,fval,exitflag]=fmincon(@fun,x0,[],[],[],[],lb,ub,@nlcon,options);

% Measurements for the nominal case I, the measurements are xB,
% x6, x10, x14, x18, xD, L, V, D, B, F, F0, in that order.
y0=[x(1);
    x(6);
    x(10);
    x(14);
    x(18);
    x(22);    
    x(23)*60;
    x(24)*60;
    x(25)*60;
    x(26)*60;
    x(27)*60;
    x(30)*60];

F0=460/60;                                  % Optimal value of F0
% dd=[0.01*0.9*ones(6,1); 0.01*F0*ones(6,1)];
dd=0.01*F0*ones(12,1);                      % Magnitude of disturbance in F0
Wd=diag(dd);                                % Diagonal scaling matrix for disturbances

pert=0.02;                                  % Magnitude of measurement error, 2%
Wn=diag([0.01; 0.01; 0.01; 0.01; 0.01; 0.01; pert*y0(7:12)]);   % Diagonal scaling matrix for measurement noise

% Finding Gy, the steady-state gain matrix
du=y0*0.01;                                 % Input change, 1%

%Finding dy/dL
lb(p.NT+1)=(y0(7)+du(7))/60;                % Increasing L with 1% from the optimal value
ub(p.NT+1)=(y0(7)+du(7))/60;                % Increasing L with 1% from the optimal value

% Running the optimization
[x,fval,exitflag]=fmincon(@fun,x0,[],[],[],[],lb,ub,@nlcon,options);
% Calculating the new measurements with step in L
yL=[x(1);
    x(6);
    x(10);
    x(14);
    x(18);
    x(22);    
    x(23)*60;
    x(24)*60;
    x(25)*60;
    x(26)*60;
    x(27)*60;
    x(30)*60];

dydL=(yL-y0)/du(7);                         % Calculating dy/dL

%Finding dy/dV
lb(p.NT+1)=0;                               % Setting L back to its optimal value
ub(p.NT+1)=Inf;                             % Setting L back to its optimal value
lb(p.NT+2)=(y0(8)+du(8))/60;                % Increasing V with 1% from the optimal value
ub(p.NT+2)=(y0(8)+du(8))/60;                % Increasing V with 1% from the optimal value

% Running the optimization
[x,fval,exitflag]=fmincon(@fun,x0,[],[],[],[],lb,ub,@nlcon,options);
% Calculating the new measurements with step in V
yV=[x(1);
    x(6);
    x(10);
    x(14);
    x(18);
    x(22);    
    x(23)*60;
    x(24)*60;
    x(25)*60;
    x(26)*60;
    x(27)*60;
    x(30)*60];

dydV=(yV-y0)/du(8);                         % Calculating dy/dV

%Finding dy/dF
lb(p.NT+2)=0;                               % Setting V back to its optimal value
ub(p.NT+2)=Inf;                             % Setting V back to its optimal value
lb(p.NT+5)=(y0(11)+du(11))/60;              % Increasing F with 1% from the optimal value
ub(p.NT+5)=(y0(11)+du(11))/60;              % Increasing F with 1% from the optimal value

% Running the optimization
[x,fval,exitflag]=fmincon(@fun,x0,[],[],[],[],lb,ub,@nlcon,options);
% Calculating the new measurements with step in F.
yF=[x(1);
    x(6);
    x(10);
    x(14);
    x(18);
    x(22);    
    x(23)*60;
    x(24)*60;
    x(25)*60;
    x(26)*60;
    x(27)*60;
    x(30)*60];

dydF=(yF-y0)/du(11);                        % Calculating dy/dV

lb(p.NT+5)=0;                               % Setting F back to its optimal value
ub(p.NT+5)=Inf;                             % Setting F back to its optimal value

Gy=[dydL dydV dydF];                        % Steady-state gain matrix dy/du

Y=(F1(:,1))'*Wd*Wn;                         % Finding the matrix Y=F'*Wd*Wn
H4=(inv(Y*Y')*Gy)';                         % Finding the matrix H
Hel=H4(1,:);                                % Selecting one of the rows
cs4=Hel*y0;                                 % Set-point of the measurement combination
%cs4=-0.0011;                               % The calculated set-point

else 
    
%% Case2
p.zF0 = zF0;
p.Vmax = Vmaxs;

% Running optimization
[x,fval,exitflag]=fmincon(@fun,x0,[],[],[],[],lb,ub,@nlcon,options);

% Measurements for the nominal case II, the measurements are xB,
% x6, x10, x14, x18, xD, L, V, D, B, F, F0, in that order.
y0=[x(1);
    x(6);
    x(10);
    x(14);
    x(18);
    x(22);    
    x(23)*60;
    x(24)*60;
    x(25)*60;
    x(26)*60;
    x(27)*60;
    x(30)*60]; 

Vm=1500/60;                                 % Optimal value of Vmax
dd=0.01*Vm*ones(12,1);                      % Magnitude of disturbance in Vmax
Wd=diag(dd);                                % Diagonal scaling matrix for disturbances                    

pert=0.02;
Wn=diag([0.01; 0.01; 0.01; 0.01; 0.01; 0.01; pert*y0(7:12)]);   % Diagonal scaling matrix for measurement noise

%Finding Gy, the steady-state gain matrix
% du=u0*0.01;
du=y0*0.01;                                 % Input change, 1%

%Finding dy/dL
lb(p.NT+1)=(y0(7)+du(7))/60;                % Increasing L with 1% from the optimal value
ub(p.NT+1)=(y0(7)+du(7))/60;                % Increasing L with 1% from the optimal value

% Running the optimization
[x,fval,exitflag]=fmincon(@fun,x0,[],[],[],[],lb,ub,@nlcon,options);

% Calculating the new measurements with step in L.
yL=[x(1);
    x(6);
    x(10);
    x(14);
    x(18);
    x(22);    
    x(23)*60;
    x(24)*60;
    x(25)*60;
    x(26)*60;
    x(27)*60;
    x(30)*60];

dydL=(yL-y0)/du(7);                         % Calculating dy/dL

%Finding dy/dF
lb(p.NT+1)=0;                               % Setting L back to its optimal value
ub(p.NT+1)=Inf;                             % Setting L back to its optimal value
lb(p.NT+5)=(y0(11)+du(11))/60;              % Increasing F with 1% from the optimal value
ub(p.NT+5)=(y0(11)+du(11))/60;              % Increasing F with 1% from the optimal value

% Running the optimization
[x,fval,exitflag]=fmincon(@fun,x0,[],[],[],[],lb,ub,@nlcon,options);

% Calculating the new measurements with step in F
yF=[x(1);
    x(6);
    x(10);
    x(14);
    x(18);
    x(22);    
    x(23)*60;
    x(24)*60;
    x(25)*60;
    x(26)*60;
    x(27)*60;
    x(30)*60];

dydF=(yF-y0)/du(11);                        % Calculating dy/dF

%Finding dy/dF0
lb(p.NT+5)=0;                               % Setting F back to its optimal value
ub(p.NT+5)=Inf;                             % Setting F back to its optimal value
lb(p.NT+8)=(y0(12)+du(12))/60;              % Increasing F0 with 1% from the optimal value
ub(p.NT+8)=(y0(12)+du(12))/60;              % Increasing F0 with 1% from the optimal value

% Running the optimization
[x,fval,exitflag]=fmincon(@fun,x0,[],[],[],[],lb,ub,@nlcon,options);

% Calculating the new measurements
yF0=[x(1);
    x(6);
    x(10);
    x(14);
    x(18);
    x(22);    
    x(23)*60;
    x(24)*60;
    x(25)*60;
    x(26)*60;
    x(27)*60;
    x(30)*60];

dydF0=(yF0-y0)/du(12);                      % Calculating dy/dF0

lb(p.NT+8)=0;                               % Setting F0 back to its optimal value
ub(p.NT+8)=Inf;                             % Setting F0 back to its optimal value

Gy=[dydL dydF dydF0];                       % Steady-state gain matrix dy/du
Y=F2(1:12)*Wd*Wn;                           % Finding the matrix Y=F'*Wd*Wn
H5=(inv(Y*Y')*Gy)';                         % Finding the matrix H
Hel2=H5(1,:);                               % Selecting one of the rows
cs5=Hel2*y0                                 % Set-point of the measurement combination
%cs5=-0.0030                                % The calculated set-point
end

