%close all
clc
clear all
%Try to optimize distillation column "A" for different disturbances!

% Define nominal inputs!

% Nominal inputs
LT=2.70629;                          % Reflux
VB=3.20629;                          % Boilup
D=0.5;                               % Distillate
B=0.5;                               % Bottoms
F=1.0;                               % Feedrate
zF=0.5;                              % Feed composition
qF=1.0;                              % Feed liquid fraction
Uinit = [LT VB D B F zF qF]';
global NT d
NT = 41 ;

%Now we split this set in three: LT and VB are our unused inputs which can
%be used for optimization. D and B are used to control the holdups in the
%column. F, zF and qF are considered disturbances (but we leave qF
%unchanged for now).

u0 = Uinit(1:2) ; % L and V
pV = 0.01;
F = 1.4 ;
d = [F zF qF pV]' ;

xb = 0.007:0.0001:0.011 ;
J = zeros(size(xb));
L = J ;
V = J ;
Vmax = u0(2)*1.25*ones(size(xb));

options = optimset('Display','off') ;
for i=1:length(xb)
     uopt = fsolve(@col_eqconstr,u0,options,xb(i)) ;
     J(i) = col_obj(uopt) ;
     L(i) = uopt(1) ; V(i) = uopt(2) ;
     disp(['Calculation number ' num2str(i) ' finished'])
end

figure(1),plot((1-xb),L,'-r',(1-xb),V,'b',(1-xb),Vmax,'gr')
figure(2),plot((1-xb(1:14)),J(1:14),'--',(1-xb(14:31)),J(14:31),'-',(1-xb(31:end)),J(31:end),'--');

pV = 0.0135 ; 
d = [F zF qF pV]' ;
for i=1:length(xb)
     uopt = fsolve(@col_eqconstr,u0,options,xb(i)) ;
     J(i) = col_obj(uopt) ;
     L(i) = uopt(1) ; V(i) = uopt(2) ;
     disp(['Calculation number ' num2str(i) ' finished'])
end
figure(3),plot((1-xb(1:14)),J(1:14),'--',(1-xb(14:31)),J(14:31),'-',(1-xb(31:end)),J(31:end),'--');

pV = 0.018 ; 
d = [F zF qF pV]' ;
for i=1:length(xb)
     uopt = fsolve(@col_eqconstr,u0,options,xb(i)) ;
     J(i) = col_obj(uopt) ;
     L(i) = uopt(1) ; V(i) = uopt(2) ;
     disp(['Calculation number ' num2str(i) ' finished'])
end
figure(4),plot((1-xb(1:14)),J(1:14),'--',(1-xb(14:31)),J(14:31),'-',(1-xb(31:end)),J(31:end),'--');


