% 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;
p.NT    = 22;
p.NF    = 13;
p.alpha = 2;

Vmaxs=(1500/60);
p.Vmax  = Vmaxs;        % Vmax = [1200 1300 1400 1500 1600 1700 1800]
p.xB=0.0105;
p.F=958;
p.L=778;
p.D=498;
p.LD=1.6;
p.LF=0.8;
p.LV=0.61;
p.xD=0.82;
p.FF0=2;
p.DF0=2;

% CSTR parameters
% p.F0=[400/60 430/60 460/60 490/60 520/60 550/60];
F0s = (460/60);            % F0 = [368 400 430 460 490 520 550 552]
zF0 = 0.9; 

p.F0 = F0s;            % F0 = [368 400 430 460 490 520 550 552]
p.zF0 = zF0;                % zF0 = [0.80 0.85 0.90 0.95 1.0]
p.k1 = 0.341/60; 
% Flags
p.OPTI=0;
p.case_I=1;

if p.case_I==1
    % Constraints
    lb=zeros(p.NT+8,1);
    ub=[ones(p.NT,1); ones(8,1)*Inf];
    
    % xB <= 0.0105
    ub(1)=0.0105;           % xB = [0.009 0.010 0.0105 0.011 0.012]
    % Mr <= 2800;
%     lb(p.NT+7)=2800;
    ub(p.NT+7)=2800;        % Mr = [2400 2600 2800 3000 3200]
    % F0 = fixed; 
    lb(p.NT+8)=p.F0;
    ub(p.NT+8)=p.F0;
 
else
    % Constraints
    lb=zeros(p.NT+8,1);
    ub=[ones(p.NT,1); ones(8,1)*Inf];
    
    % xB <= 0.0105
    ub(1)=0.0105;           % xB = [0.009 0.010 0.0105 0.011 0.012]
    % Mr <= 2800;
    ub(p.NT+7)=2800;        % Mr = [2400 2600 2800 3000 3200]
    % V <= Vmax; 
    ub(p.NT+2)=p.Vmax;
end

% Initial values
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))

Vn=x(p.NT+2)*60
% F0n=-x(p.NT+8)*60
% break 

%%
%Optimal values of V
for i=1:7
    p.F0=(370/60):(30/60):(550/60); 
    lb(p.NT+8)=p.F0(i);
    ub(p.NT+8)=p.F0(i);
    [x,fval,exitflag]=fmincon(@fun,x0,[],[],[],[],lb,ub,@nlcon,options);
    Vopt(i)=x(p.NT+2)*60;
end
 
%% 
%Loss with null space method (Simulink)
%Disturbance in F0
% for i=1:7
%     p.F0=(370/60):(30/60):(550/60); 
%     lb(p.NT+8)=p.F0(i);
%     ub(p.NT+8)=p.F0(i);
%     [x,fval,exitflag]=fmincon(@fun,x0,[],[],[],[],lb,ub,@nlcon,options);
%     Vns(i)=x(p.NT+2)*60;
% end


%% Exact local method in Simulink
%Disturbance in F0
% for i=1:7
%    p.F0=(370/60):(30/60):(550/60); 
%    lb(p.NT+8)=p.F0(i);
%    ub(p.NT+8)=p.F0(i);
%    [x,fval,exitflag]=fmincon(@fun,x0,[],[],[],[],lb,ub,@nlcon,options);
%    Vel(i)=x(p.NT+2)*60;
% 
% end 

load Vopt.mat
load Vns.mat
load Vel.mat
F0=370:30:550;

%Plotting loss with null space method
lossNS= ((Vns-Vopt)./Vopt)*100;
figure(1)
xx=F0;
yy=spline(F0,lossNS,xx);
plot(xx,yy,'b','LineWidth',2)
xlim([370 550])
ylim([0 400])
xlabel('F_0 [kmol/h]')
ylabel('Loss [%]')
title('Null Space Method')

% %Plotting loss with exact local method
% lossEL= ((Vel-Vopt)./Vopt)*100;
% figure(2)
% xx=F0;
% yy=spline(F0,lossEL,xx);
% plot(xx,yy,'b','LineWidth',2)
% xlim([370 550])
% ylim([0 300])
% xlabel('F_0 [kmol/h]')
% ylabel('Loss [%]')
% title('Exact Local Method')

break 

%%
%Calculation of cost with the Brute force method

if p.case_I==1
   
% Calculation of the optimal cost for different F0 values.
VoptF0=zeros(1,17);
for i=1:17
    p.F0=(F0s*0.8):(11.5/60):(F0s*1.2);
    lb(p.NT+8)=p.F0(i);
    ub(p.NT+8)=p.F0(i);
    [x,fval,exitflag]=fmincon(@fun,x0,[],[],[],[],lb,ub,@nlcon,options);
    VoptF0(i)=x(p.NT+2)*60;
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

VF0_DF0=zeros(1,17);
for i=1:17
    p.F0=(F0s*0.8):(11.5/60):(F0s*1.2);
    lb(p.NT+8)=p.F0(i);
    ub(p.NT+8)=p.F0(i);
    [x,fval,exitflag]=fmincon(@fun,x0,[],[],[],[],lb,ub,@nlcon,options);
    VF0_DF0(i)=x(p.NT+2)*60;
end

% break 

%% Calculation of the optimal cost for different zF0 values.
p.F0 = F0s;                         % F0=460
lb(p.NT+8)=p.F0;
ub(p.NT+8)=p.F0;

VoptzF0=zeros(1,11);
for i=1:11
    p.zF0=0.80:0.02:1.0;
    p.zF0=p.zF0(i);
    [x,fval,exitflag]=fmincon(@fun,x0,[],[],[],[],lb,ub,@nlcon,options);
    VoptzF0(i)=x(p.NT+2)*60;
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

VzF0_DF0=zeros(1,11);
for i=1:11
    p.zF0=0.80:0.02:1.0;
    p.zF0=p.zF0(i);
    [x,fval,exitflag]=fmincon(@fun,x0,[],[],[],[],lb,ub,@nlcon,options);
    VzF0_DF0(i)=x(p.NT+2)*60;
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

% Impl.error in F
load yopt_nom_case1.mat           % Matrix with optimal values of the CVs

VieF=zeros(1,9);
for i=1:9
    p.F=A1(1)*0.8/60: A1(1)*0.05/60: A1(1)*1.2/60;
    lb(p.NT+5)=p.F(i);
    ub(p.NT+5)=p.F(i);
    [x,fval,exitflag]=fmincon(@fun,x0,[],[],[],[],lb,ub,@nlcon,options);
    VieF(i)=x(p.NT+2)*60;
end 
lb(p.NT+5)=0;
ub(p.NT+5)=Inf;
 
% Impl.error in L
VieL=zeros(1,9);
for i=1:9
    p.L=A1(2)*0.8/60: A1(2)*0.05/60: A1(2)*1.2/60;
    lb(p.NT+1)=p.L(i);
    ub(p.NT+1)=p.L(i);
    [x,fval,exitflag]=fmincon(@fun,x0,[],[],[],[],lb,ub,@nlcon,options);
    VieL(i)=x(p.NT+2)*60;
end 
lb(p.NT+1)=0;
ub(p.NT+1)=Inf;

% Impl.error in D
VieD=zeros(1,9);
for i=1:9
    p.D=A1(3)*0.8/60: A1(3)*0.05/60: A1(3)*1.2/60;
    lb(p.NT+3)=p.D(i);
    ub(p.NT+3)=p.D(i);
    [x,fval,exitflag]=fmincon(@fun,x0,[],[],[],[],lb,ub,@nlcon,options);
    VieD(i)=x(p.NT+2)*60;
end 
lb(p.NT+3)=0;
ub(p.NT+3)=Inf;

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

% Impl.error in xD
ViexD=zeros(1,9);
for i=1:9
    p.xD=A1(7)*0.8: A1(7)*0.05: A1(7)*1.2;
    lb(p.NT)=p.xD(i);
    ub(p.NT)=p.xD(i);
    [x,fval,exitflag]=fmincon(@fun,x0,[],[],[],[],lb,ub,@nlcon,options);
    ViexD(i)=x(p.NT+2)*60;
end 
lb(p.NT)=0;
ub(p.NT)=1;

% break 

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

VoptMr=zeros(1,11);
for i=1:11
    p.Mr=2300:50:2800;
    lb(p.NT+7)=p.Mr(i);
    ub(p.NT+7)=p.Mr(i);
    [x,fval,exitflag]=fmincon(@fun,x0,[],[],[],[],lb,ub,@nlcon,options);
    VoptMr(i)=x(p.NT+2)*60;
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

VMr_DF0=zeros(1,11);
for i=1:11
    p.Mr=2300:50:2800;
    lb(p.NT+7)=p.Mr(i);
    ub(p.NT+7)=p.Mr(i);
    [x,fval,exitflag]=fmincon(@fun,x0,[],[],[],[],lb,ub,@nlcon,options);
    VMr_DF0(i)=x(p.NT+2)*60;
end 

%%
% Optimal cost for implementation errors in xB
p.F0=F0s;
p.zF0=zF0;
ub(p.NT+7)=2800;

VoptxB=zeros(1,11);
for i=1:11
    p.xB=0.009:0.00015:0.0105;
    lb(1)=p.xB(i);
    ub(1)=p.xB(i);
    [x,fval,exitflag]=fmincon(@fun,x0,[],[],[],[],lb,ub,@nlcon,options);
    VoptxB(i)=x(p.NT+2)*60;
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

VxB_DF0=zeros(1,11);
for i=1:11
    p.xB=0.009:0.00015:0.0105;
    lb(1)=p.xB(i);
    ub(1)=p.xB(i);
    [x,fval,exitflag]=fmincon(@fun,x0,[],[],[],[],lb,ub,@nlcon,options);
    VxB_DF0(i)=x(p.NT+2)*60;
end

else
    
  
% Calculation of the optimal cost for different Vmax values, case II.

% F0oVm=zeros(1,13);
% for i=1:13
%     p.Vmax=1200/60:50/60:1800/60;
%     lb(p.NT+2)=p.Vmax(i);
%     ub(p.NT+2)=p.Vmax(i);
%     [x,fval,exitflag]=fmincon(@fun,x0,[],[],[],[],lb,ub,@nlcon,options);
%     F0oVm(i)=-x(p.NT+8)*60;
% 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

% F0Vm_DF0=zeros(1,13);
% for i=1:13
%     p.Vmax=1200/60:50/60:1800/60;
%     lb(p.NT+2)=p.Vmax(i);
%     ub(p.NT+2)=p.Vmax(i);
%     [x,fval,exitflag]=fmincon(@fun,x0,[],[],[],[],lb,ub,@nlcon,options);
%     F0Vm_DF0(i)=-x(p.NT+8)*60;
% end
% break 

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

% F0ozF0=zeros(1,11);
% for i=1:11
%     p.zF0=0.80:0.02:1.0;
%     p.zF0=p.zF0(i);
%     [x,fval,exitflag]=fmincon(@fun,x0,[],[],[],[],lb,ub,@nlcon,options);
%     F0ozF0(i)=-x(p.NT+8)*60;
% 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_DF0=zeros(1,11);
% for i=1:11
%     p.zF0=0.80:0.02:1.0;
%     p.zF0=p.zF0(i);
%     [x,fval,exitflag]=fmincon(@fun,x0,[],[],[],[],lb,ub,@nlcon,options);
%     F0zF0_DF0(i)=-x(p.NT+8)*60;
% end

% break 
%% Calculation of the cost for implementation errors in c (CVs).
p.zF0 = zF0;                        % zF0=0.9
p.Vmax=Vmaxs;                       % Vmax=1500
load yopt_nom_case2.mat           % Matrix with optimal values of the CVs for case 2

% Impl.error in F
% F0ieF=zeros(1,9);
% for i=1:9
%     p.F=A4(1)*0.8/60: A4(1)*0.05/60: A4(1)*1.2/60;
%     lb(p.NT+5)=p.F(i);
%     ub(p.NT+5)=p.F(i);
%     [x,fval,exitflag]=fmincon(@fun,x0,[],[],[],[],lb,ub,@nlcon,options);
%     F0ieF(i)=-x(p.NT+8)*60;
% end 
% lb(p.NT+5)=0;
% ub(p.NT+5)=Inf;
%  
% % Impl.error in L
% F0ieL=zeros(1,9);
% for i=1:9
%     p.L=A4(2)*0.8/60: A4(2)*0.05/60: A4(2)*1.2/60;
%     lb(p.NT+1)=p.L(i);
%     ub(p.NT+1)=p.L(i);
%     [x,fval,exitflag]=fmincon(@fun,x0,[],[],[],[],lb,ub,@nlcon,options);
%     F0ieL(i)=-x(p.NT+8)*60;
% end 
% lb(p.NT+1)=0;
% ub(p.NT+1)=Inf;
% 
% % Impl.error in D
% F0ieD=zeros(1,9);
% for i=1:9
%     p.D=A4(3)*0.8/60: A4(3)*0.05/60: A4(3)*1.2/60;
%     lb(p.NT+3)=p.D(i);
%     ub(p.NT+3)=p.D(i);
%     [x,fval,exitflag]=fmincon(@fun,x0,[],[],[],[],lb,ub,@nlcon,options);
%     F0ieD(i)=-x(p.NT+8)*60;
% end 
% lb(p.NT+3)=0;
% ub(p.NT+3)=Inf;
% 
% % Impl.error in xD
% F0iexD=zeros(1,9);
% for i=1:9
%     p.xD=A4(7)*0.8: A4(7)*0.05: A4(7)*1.2;
%     lb(p.NT)=p.xD(i);
%     ub(p.NT)=p.xD(i);
%     [x,fval,exitflag]=fmincon(@fun,x0,[],[],[],[],lb,ub,@nlcon,options);
%     F0iexD(i)=-x(p.NT+8)*60;
% end 
% lb(p.NT)=0;
% ub(p.NT)=1;
% 
% % 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;

% F0oMr=zeros(1,11);
% for i=1:11
%     p.Mr=2300:50:2800;
%     lb(p.NT+7)=p.Mr(i);
%     ub(p.NT+7)=p.Mr(i);
%     [x,fval,exitflag]=fmincon(@fun,x0,[],[],[],[],lb,ub,@nlcon,options);
%     F0oMr(i)=-x(p.NT+8)*60;
% 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

F0Mr_F=zeros(1,11);
for i=1:11
    p.Mr=2300:50:2800;
    lb(p.NT+7)=p.Mr(i);
    ub(p.NT+7)=p.Mr(i);
    [x,fval,exitflag]=fmincon(@fun,x0,[],[],[],[],lb,ub,@nlcon,options);
    F0Mr_F(i)=-x(p.NT+8)*60;
end 
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);
% for i=1:11
%     p.xB=0.009:0.00015:0.0105;
%     lb(1)=p.xB(i);
%     ub(1)=p.xB(i);
%     [x,fval,exitflag]=fmincon(@fun,x0,[],[],[],[],lb,ub,@nlcon,options);
%     F0oxB(i)=-x(p.NT+8)*60;
% 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);
% for i=1:11
%     p.xB=0.009:0.00015:0.0105;
%     lb(1)=p.xB(i);
%     ub(1)=p.xB(i);
%     [x,fval,exitflag]=fmincon(@fun,x0,[],[],[],[],lb,ub,@nlcon,options);
%     F0xB_DF0(i)=-x(p.NT+8)*60;
% end
% break

end 

%%
%%Null space method
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];    %Measurements for the nominal case 1 or 2

if p.case_I==1
%% dd1
pert=0.01;
dd1=p.F0*pert;
p.F0 = F0s+dd1;            % F0 = [368 400 430 460 490 520 550 552]

lb(p.NT+8)=p.F0;
ub(p.NT+8)=p.F0;

[x,fval,exitflag]=fmincon(@fun,x0,[],[],[],[],lb,ub,@nlcon,options);
ydd1=[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];

%% dd2
p.F0 = F0s; % F0 = [368 400 430 460 490 520 550 552]
dd2=zF0*pert;

p.zF0 = zF0 +dd2;                % zF0 = [0.80 0.85 0.90 0.95 1.0]

lb(p.NT+8)=p.F0;
ub(p.NT+8)=p.F0;

[x,fval,exitflag]=fmincon(@fun,x0,[],[],[],[],lb,ub,@nlcon,options);
ydd2=[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];

else
%% dd3
pert=0.01;
dd3=p.Vmax*pert;
p.Vmax = Vmaxs+dd3;  

ub(p.NT+2)=p.Vmax;

[x,fval,exitflag]=fmincon(@fun,x0,[],[],[],[],lb,ub,@nlcon,options);
ydd3=[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];
%%
%% dd4
p.Vmax=Vmaxs; 
dd4=zF0*pert;
p.zF0 = zF0 +dd4;

ub(p.NT+2)=p.Vmax;

[x,fval,exitflag]=fmincon(@fun,x0,[],[],[],[],lb,ub,@nlcon,options);
ydd4=[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];
end

if p.case_I==1
F1=[(ydd1-y0)/dd1 (ydd2-y0)/dd2];     % F matrix for case I
H1=null(F1');
H1=H1';
H=H1(1,:);
% % 
% cs1=H*y0
% 
else
% %%
% %Case II: Maximum F0 and F, disturbance d3=Vmax and d4=zF0
% % Two disturbances
F2=[(ydd3-y0)/dd3 (ydd4-y0)/dd4];     % F matrix for case II
H2=null(F2');
H2=H2';
H2=H2(4,:);
% 
% cs2=H2*y0
%cs2=0.2907
end

%% Exact local method 
if p.case_I==1
% Case1
p.F0 = F0s;
p.zF0 = zF0;

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];   % Candidate CVs for the nominal case 1 or 2

Wd=2;
Wn=diag([0.01;
    0.01;
    0.01;
    0.01;
    0.01;
    0.01;
    0.01;
    0.01;
    0.01;
    0.01;
    0.01;
    0.01]);         % Measurement noise matrix

%Finding Gy
% du=u0*0.01;
du=y0*0.01;

%Finding dy/dL
lb(p.NT+1)=(y0(7)+du(7))/60;
ub(p.NT+1)=(y0(7)+du(7))/60;


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

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

%Finding dy/dV
lb(p.NT+1)=0;
ub(p.NT+1)=Inf;
lb(p.NT+2)=(y0(8)+du(8))/60;
ub(p.NT+2)=(y0(8)+du(8))/60;


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

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

%Finding dy/dD
lb(p.NT+1)=0;
ub(p.NT+1)=Inf;
lb(p.NT+2)=0;
ub(p.NT+2)=Inf;
lb(p.NT+3)=(y0(9)+du(9))/60;
ub(p.NT+3)=(y0(9)+du(9))/60;


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

yD=[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];

dydD=(yD-y0)/du(9);

%Finding dy/dB
lb(p.NT+1)=0;
ub(p.NT+1)=Inf;
lb(p.NT+2)=0;
ub(p.NT+2)=Inf;
lb(p.NT+3)=0;
ub(p.NT+3)=Inf;
lb(p.NT+4)=y0(10)+du(10);
ub(p.NT+4)=y0(10)+du(10);

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

yB=[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];

dydB=(yB-y0)/du(10);

%Finding dy/dF
lb(p.NT+1)=0;
ub(p.NT+1)=Inf;
lb(p.NT+2)=0;
ub(p.NT+2)=Inf;
lb(p.NT+3)=0;
ub(p.NT+3)=Inf;
lb(p.NT+4)=0;
ub(p.NT+4)=Inf;
lb(p.NT+5)=y0(11)+du(11);
ub(p.NT+5)=y0(11)+du(11);

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

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

%Finding dy/dF0
lb(p.NT+1)=0;
ub(p.NT+1)=Inf;
lb(p.NT+2)=0;
ub(p.NT+2)=Inf;
lb(p.NT+3)=0;
ub(p.NT+3)=Inf;
lb(p.NT+4)=0;
ub(p.NT+4)=Inf;
lb(p.NT+5)=0;
ub(p.NT+5)=Inf;
lb(p.NT+8)=y0(12)+du(12);
ub(p.NT+8)=y0(12)+du(12);

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

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

Gy=[dydL dydV dydD dydB dydF dydF0];

Y=F1(1:12)*Wd*Wn;
H4=(inv(Y*Y')*Gy)';
Hel=H4(2,:);
% cs4=Hel*y0
%cs4=22.8872

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

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];   % Candidate CVs for the nominal case 2

Wd=2;
Wn=diag([0.01;
    0.01;
    0.01;
    0.01;
    0.01;
    0.01;
    0.01;
    0.01;
    0.01;
    0.01;
    0.01;
    0.01]);         % Measurement noise matrix

%Finding Gy
% du=u0*0.01;
du=y0*0.01;

%Finding dy/dL
lb(p.NT+1)=(y0(7)+du(7))/60;
ub(p.NT+1)=(y0(7)+du(7))/60;


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

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

%Finding dy/dV
lb(p.NT+1)=0;
ub(p.NT+1)=Inf;
lb(p.NT+2)=(y0(8)+du(8))/60;
ub(p.NT+2)=(y0(8)+du(8))/60;


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

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

%Finding dy/dD
lb(p.NT+1)=0;
ub(p.NT+1)=Inf;
lb(p.NT+2)=0;
ub(p.NT+2)=Inf;
lb(p.NT+3)=(y0(9)+du(9))/60;
ub(p.NT+3)=(y0(9)+du(9))/60;


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

yD=[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];

dydD=(yD-y0)/du(9);

%Finding dy/dB
lb(p.NT+1)=0;
ub(p.NT+1)=Inf;
lb(p.NT+2)=0;
ub(p.NT+2)=Inf;
lb(p.NT+3)=0;
ub(p.NT+3)=Inf;
lb(p.NT+4)=y0(10)+du(10);
ub(p.NT+4)=y0(10)+du(10);

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

yB=[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];

dydB=(yB-y0)/du(10);

%Finding dy/dF
lb(p.NT+1)=0;
ub(p.NT+1)=Inf;
lb(p.NT+2)=0;
ub(p.NT+2)=Inf;
lb(p.NT+3)=0;
ub(p.NT+3)=Inf;
lb(p.NT+4)=0;
ub(p.NT+4)=Inf;
lb(p.NT+5)=y0(11)+du(11);
ub(p.NT+5)=y0(11)+du(11);

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

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

%Finding dy/dF0
lb(p.NT+1)=0;
ub(p.NT+1)=Inf;
lb(p.NT+2)=0;
ub(p.NT+2)=Inf;
lb(p.NT+3)=0;
ub(p.NT+3)=Inf;
lb(p.NT+4)=0;
ub(p.NT+4)=Inf;
lb(p.NT+5)=0;
ub(p.NT+5)=Inf;
lb(p.NT+8)=y0(12)+du(12);
ub(p.NT+8)=y0(12)+du(12);

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

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

Gy=[dydL dydV dydD dydB dydF dydF0];

Y=F2(1:12)*Wd*Wn;
H5=(inv(Y*Y')*Gy)';
Hel2=H5(2,:);
% cs5=Hel2*y0
%cs5=1081.3
end
