function albe
%% Albumin - Ethanol Separation : Economically optimal 
%% Choose case : 
case39 = input('Choose the case number (3,9): ');

if (case39 == 3)
    %% Case 3:
    % initial and final conditions :
    x0 = [15; 194.3];  param.c10 = x0(1);  param.c20 = x0(2);   % [kg/m3]
    xf = [80; 0.1];    param.c1f = xf(1);  param.c2f = xf(2);   % [kg/m3]
elseif(case39 == 9)
    %% Case 9:
    % initial and final conditions :
    x0 = [15; 194.3];  param.c10 = x0(1);  param.c20 = x0(2);   % [kg/m3]
    xf = [240; 0.1];   param.c1f = xf(1);  param.c2f = xf(2);   % [kg/m3]
else
    disp('Wrong choice !!!')
    return;
end
% initial volume : 
param.V0 = 0.0666;                                              % [m3]

% parameters : 
param.b = [2.877; 1.698e-1; 1.874e-2; 5.708e-4; -2.394e-4; 9.334e-5];
param.K = 1/(param.c10*param.V0);
param.N = 10000000000;

% weight coefficients : 
param.wT = 2;
param.wD = 10;

% control :
alpha0   = 0;
alpha1   = 1;
alphainf = 1000e1;
alphaVVD = 1/(1+((log(param.c1f) - log(param.c10))/(log(param.c20) - log(param.c2f))));

% Tolerances :
AbsTol = 1e-13;
RelTol = 1e-12;

%% minimum time : 

x01t      = [param.c10;param.c20;0];
options1t = odeset('Events',@event1t,'AbsTol',AbsTol,'RelTol',RelTol);
tspan1t   = [0 param.N];
[time1t,state1t,te1t] = ode45(@model,tspan1t,x01t,options1t,alpha1,param);

options2t = odeset('Events',@event2t,'AbsTol',AbsTol,'RelTol',RelTol);
x02t      = state1t(end,:);
tspan2t   = [time1t(end) param.N];
[time2t,state2t,te2t] = ode45(@modelsing,tspan2t,x02t,options2t,param);

options3t = odeset('Events',@event3t,'AbsTol',AbsTol,'RelTol',RelTol);
x03t      = state2t(end,:);
tspan3t   = [time2t(end) param.N];
[time3t,state3t,te3t] = ode45(@model,tspan3t,x03t,options3t,alpha1,param);

%% multi-objective : 

x01m      = [param.c10;param.c20;0];
options1m = odeset('Events',@event1m,'AbsTol',AbsTol,'RelTol',RelTol);
tspan1m   = [0 param.N];
[time1m,state1m,te1m] = ode45(@model,tspan1m,x01m,options1m,alpha0,param);           

options2m = odeset('Events',@event2m,'AbsTol',AbsTol,'RelTol',RelTol);
x02m      = state1m(end,:);
tspan2m   = [time1m(end) param.N];
[time2m,state2m,te2m] = ode45(@modelsing,tspan2m,x02m,options2m,param);

options3m = odeset('Events',@event3m,'AbsTol',AbsTol,'RelTol',RelTol);
x03m      = state2m(end,:);
tspan3m   = [time2m(end) param.N];
[time3m,state3m,te3m] = ode45(@model,tspan3m,x03m,options3m,alpha1,param);

%% minimum diluant : 

x01d      = [param.c10;param.c20;0];
options1d = odeset('Events',@event1d,'AbsTol',AbsTol,'RelTol',RelTol);
tspan1d   = [0 param.N];
[time1d,state1d,te1d] = ode45(@model,tspan1d,x01d,options1d,alpha0,param);           

options2d = odeset('Events',@event2d,'AbsTol',AbsTol,'RelTol',RelTol);
x02d      = state1d(end,:);
tspan2d   = [time1d(end) param.N];
[time2d,state2d,te2d] = ode45(@model,tspan2d,x02d,options2d,alpha1,param);
 
% %% VVD :
% 
% x01VVD      = [param.c10;param.c20;0];
% optionsVVD = odeset('Events',@eventVVD,'AbsTol',AbsTol,'RelTol',RelTol);
% tspanVVD   = [0 param.N];
% [timeVVD,stateVVD,teVVD] = ode45(@model,tspanVVD,x01VVD,optionsVVD,alphaVVD,param);
% 
% %% C-CVD :
% 
% options1CVD = odeset('Events',@event1CVD,'AbsTol',AbsTol,'RelTol',RelTol);
% c01CVD      = [param.c10;param.c20;0];
% tspan1CVD   = [0 10000];
% [time1CVD,state1CVD,te1CVD] = ode45(@model,tspan1CVD,c01CVD,options1CVD,alpha0,param);
% 
% options2CVD = odeset('Events',@event2CVD,'AbsTol',AbsTol,'RelTol',RelTol);
% c02CVD      = state1CVD(end,:);
% tspan2CVD   = [time1CVD(end) 10000];
% [time2CVD,state2CVD,te2CVD] = ode45(@model,tspan2CVD,c02CVD,options2CVD,alpha1,param);

% %% plot flux : 
% 
% time  = [time1t;time2t;time3t]/time3t(end);
% c1    = [state1t(:,1);state2t(:,1);state3t(:,1)]; 
% c2    = [state1t(:,2);state2t(:,2);state3t(:,2)];
% 
% q = 1./(param.b(1) + param.b(2).*c1 + param.b(3).*c2 + ...
%     param.b(4).*c1.*c2 + param.b(5).*c1.^2 + param.b(6).*c2.^2);
% 
% outflow = c2.*q;
% 
% plot(time,outflow)
% xlabel('time')
% ylabel('c_2*q(c_1,c_2)')
% title('Case B')
%% state diagram : 
line        = 1.5;
marker_sizes = 6;
marker_size  = 9;

b  = [2.877; 1.698e-1; 1.874e-2; 5.708e-4; -2.394e-4; 9.334e-5];
c1 = [10:5:250];
c2t = -(param.b(4).*c1 - (param.b(4)^2*c1.^2 - 4*param.b(5)*param.b(6)*...
       c1.^2 + 4*param.b(1)*param.b(6)).^(1/2))/(2*param.b(6));

c2m = ((param.wT*(param.wT*param.b(4)^2*c1.^2 - 4*param.b(5)*param.b(6)...
      *param.wT*c1.^2 + 4*param.b(6)*param.wD + 4*param.b(1)*param.b(6)...
      *param.wT)).^(1/2) - param.b(4)*c1*param.wT)/(2*param.b(6)*param.wT);

figure
plot([state1d(:,1);state2d(:,1)],...
     [state1d(:,2);state2d(:,2)],'g--','linewidth',line)                     % minimum diluant
hold on
plot([state1m(:,1);state2m(:,1);state3m(:,1)],...
     [state1m(:,2);state2m(:,2);state3m(:,2)],'r','linewidth',line)        % economically optimal
plot([state1t(:,1);state2t(:,1);state3t(:,1)],...
     [state1t(:,2);state2t(:,2);state3t(:,2)],'-.','linewidth',line)            % minimum time




plot(c1,c2t,'bo','MarkerEdgeColor','b','MarkerFaceColor','b',...
                                       'markersize',marker_sizes)
                                   
plot(c1,c2m,'rd','MarkerEdgeColor','r','MarkerFaceColor','r',... 
'markersize',marker_sizes)
plot(param.c10,param.c20,'go',param.c1f,param.c2f,'rx','markersize',marker_size)
legend('minimum diluant','economically optimal','minimum time',...
       'S_T(c_1,c_2) = 0','S_E(c_1,c_2) = 0','Location','SouthWest')
xlabel('c1')
ylabel('c2')
if (case39 == 3)
    axis([9.2 100.6 -5 220])
%    fs=6; fz=11;
%    exportfig(gcf,'AlbEth_case3_state.eps','width',fs,'fontmode','fixed',...
%    'fontsize',fz,'Color','rgb');
set(gcf, 'PaperUnits', 'centimeters') 
set(gcf, 'PaperPosition', [3 10 16 10])
print -depsc 'AlbEth_case3_state.eps'
elseif(case39 == 9)
    axis([8.8 252.5 -5 220])
%    fs=6; fz=11;
%    exportfig(gcf,'AlbEth_case9_state.eps','width',fs,'fontmode','fixed',...
%    'fontsize',fz,'Color','rgb');
set(gcf, 'PaperUnits', 'centimeters') 
set(gcf, 'PaperPosition', [3 10 16 10])
print -depsc 'AlbEth_case9_state.eps'
   
end


%% control diagram :     
figure
if (case39 == 9)
% minimum diluant : 
a10d     = alpha0*ones(length(time1d),1);
a21d     = alpha1*ones(length(time2d),1);

plot([time1d;time2d],[a10d;a21d],'g--','linewidth',line)
axis([0 2.7 -0.1 1.1])
hold on

% economically optimal :
a10m     = alpha0*ones(length(time1m),1);
a31m     = alpha1*ones(length(time3m),1);

alphasgm = (0.5.*param.b(4).*state2m(:,1).*state2m(:,2) + param.b(5).*...
            state2m(:,1).^2)./(param.b(4).*state2m(:,1).*state2m(:,2) + ...
            param.b(5).*state2m(:,1).^2 + param.b(6).*state2m(:,2).^2);
pom  = find((alphasgm >= 0) & (alphasgm <= 1));
pozm = pom(end);
a20m = alpha0*ones(length(time2m(pozm:end)),1);

plot([time1m;time2m(1:pozm);time2m(pozm:end);time3m],[a10m; ...
      alphasgm(1:pozm); a20m; a31m],'r','linewidth',line)

% minimum time : 
a11t     = alpha1*ones(length(time1t),1);
a31t     = alpha1*ones(length(time3t),1);
alphasgt = (0.5.*param.b(4).*state2t(:,1).*state2t(:,2) + param.b(5).*...
            state2t(:,1).^2)./(param.b(4).*state2t(:,1).*state2t(:,2) + ...
            param.b(5).*state2t(:,1).^2 + param.b(6).*state2t(:,2).^2);

pot = find((alphasgt >= 0) & (alphasgt <= 1));
pozt = pot(end);
a20t = alpha0*ones(length(time2t(pozt:end)),1);        

plot([time1t;time2t(1:pozt);time2t(pozt:end);time3t],[a11t;...
      alphasgt(1:pozt); a20t; a31t],'-.','linewidth',line)
legend('minimum diluant','economically optimal','minimum time',...
       'Location','SouthEast')
   xlabel('time')
ylabel('alpha')
%fs=6; fz=11;
%exportfig(gcf,'AlbEth_case9_control.eps','width',fs,'fontmode','fixed',...
%'fontsize',fz,'Color','rgb');
set(gcf, 'PaperUnits', 'centimeters') 
set(gcf, 'PaperPosition', [3 10 16 10])
print -depsc 'AlbEth_case9_control.eps'

end
if (case39 == 3)
% minimum diluant : 
a10d     = alpha0*ones(length(time1d),1);
a21d     = alpha1*ones(length(time2d),1);

plot([time1d;time2d],[a10d;a21d],'g--','linewidth',line)
axis([0 te2d -0.1 1.1])
hold on
% economically optimal :
a10m     = alpha0*ones(length(time1m),1);
a31m     = alpha1*ones(length(time3m),1);
alphasgm = (0.5.*param.b(4).*state2m(:,1).*state2m(:,2) + param.b(5).*...
            state2m(:,1).^2)./(param.b(4).*state2m(:,1).*state2m(:,2) + ...
            param.b(5).*state2m(:,1).^2 + param.b(6).*state2m(:,2).^2);
        
plot([time1m;time2m;time3m],[a10m; alphasgm; a31m],'r','linewidth',line)
% minimum time :        
a11t      = alpha1*ones(length(time1t),1);
a31t      = alpha1*ones(length(time3t),1);
alphasgt = (0.5.*param.b(4).*state2t(:,1).*state2t(:,2) + param.b(5).*...
            state2t(:,1).^2)./(param.b(4).*state2t(:,1).*state2t(:,2) + ...
            param.b(5).*state2t(:,1).^2 + param.b(6).*state2t(:,2).^2);

plot([time1t;time2t;time3t],[a11t; alphasgt; a31t],'-.','linewidth',line)
legend('minimum diluant','economically optimal','minimum time',...
       'Location','SouthEast')
   xlabel('time')
ylabel('alpha') 
%fs=6; fz=11;
%exportfig(gcf,'AlbEth_case3_control.eps','width',fs,'fontmode','fixed',...
%'fontsize',fz,'Color','rgb');
set(gcf, 'PaperUnits', 'centimeters') 
set(gcf, 'PaperPosition', [3 10 16 10])
print -depsc 'AlbEth_case3_control.eps'

end

function fun = model(t,x,u,param)
c1 = x(1); c2 = x(2); alpha = u;

q = 1/(param.b(1) + param.b(2)*c1 + param.b(3)*c2 + ...
    param.b(4)*c1*c2 + param.b(5)*c1^2 + param.b(6)*c2^2);

dc1dt = param.K*c1^2*q*(1 - alpha);
dc2dt = -c1*c2*param.K*q*alpha;
daqdt = alpha*q;

fun = [dc1dt;dc2dt;daqdt];

function fun = modelsing(t,x,param)
c1 = x(1); c2 = x(2);

alphasig = (0.5*param.b(4)*c1*c2 + param.b(5)*c1^2)/(param.b(4)*c1*c2 + ...
            param.b(5)*c1^2 + param.b(6)*c2^2);

if ((alphasig <= 0) || (alphasig >= 1))    
   alpha = 0;
else
   alpha = alphasig;
end

q = 1/(param.b(1) + param.b(2)*c1 + param.b(3)*c2 + ...
    param.b(4)*c1*c2 + param.b(5)*c1^2 + param.b(6)*c2^2);

dc1dt = param.K*c1^2*q*(1 - alpha);
dc2dt = -c1*c2*param.K*q*alpha;
daqdt = alpha*q;

fun = [dc1dt;dc2dt;daqdt];

function [value,isterminal,direction] = eventVVD(t,y,u,param)
value      = y(1) - param.c1f;
isterminal = 1;
direction  = 0;

function [value,isterminal,direction] = event1CVD(t,y,u,param)
value      = y(1) - param.c1f;
isterminal = 1;
direction  = 1;

function [value,isterminal,direction] = event1d(t,y,u,param)
value      = y(1) - param.c1f;
isterminal = 1;
direction  = 1;

function [value,isterminal,direction] = event1m(t,y,u,param)
value      = y(1) - (((param.wT*(param.wT*param.b(4)^2*y(2)^2 - 4*param.b(5)*param.b(6)*param.wT*y(2)^2 + 4*param.b(5)*param.wD + 4*param.b(1)*param.b(5)*param.wT))^(1/2) - param.b(4)*y(2)*param.wT)/(2*param.b(5)*param.wT));
isterminal = 1;
direction  = 0;

function [value,isterminal,direction] = event1t(t,y,u,param)
value      = y(1) - (-(param.b(4)*y(2) - (param.b(4)^2*y(2)^2 - 4*param.b(5)*param.b(6)*y(2)^2 + 4*param.b(1)*param.b(5))^(1/2))/(2*param.b(5)));
isterminal = 1;
direction  = -1;

function [value,isterminal,direction] = event2CVD(t,y,u,param)
value      = y(2) - param.c2f;
isterminal = 1;
direction  = -1;

function [value,isterminal,direction] = event2d(t,y,u,param)
value      = y(2) - param.c2f;
isterminal = 1;
direction  = -1;

function [value,isterminal,direction] = event2m(t,y,param)
value      = y(1) - param.c1f;
isterminal = 1;
direction  = 1;

function [value,isterminal,direction] = event2t(t,y,param)
value      = y(1) - param.c1f;
isterminal = 1;
direction  = 1;

function [value,isterminal,direction] = event3m(t,y,u,param)
value      = y(2) - param.c2f;
isterminal = 1;
direction  = -1;

function [value,isterminal,direction] = event3t(t,y,u,param)
value      = y(2) - param.c2f;
isterminal = 1;
direction  = -1;