function alboptc3pf
%% Albumin - Ethanol Separation :
%% 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]

% 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 = 1;
wD       = logspace(-6,6,500);

% calculation of maximum wD
wTwD  = -(4*param.b(5))/(param.b(4)^2*param.c20^2 - ...
          4*param.b(5)*param.b(6)*param.c20^2 + ...
          4*param.b(1)*param.b(5));
wDmax = 1/wTwD;

% find all wD < wDmax 
result = find(wD < wDmax);
wDsing = wD(1:result(end)-1);
wD01   = wD(result(end):end);

% control :
alpha0   = 0;
alpha1   = 1;
alphainf = 1000e1;

% Tolerances :
AbsTol = 1e-8;
RelTol = 1e-6;

% memory allocation :
wTwD      = zeros(length(wDsing),1);
tff       = zeros(length(wDsing),1);
Vw        = zeros(length(wDsing),1);
S         = zeros(length(wDsing),1);
time1tf   = zeros(length(wDsing),1);
time2tf   = zeros(length(wDsing),1);
time3tf   = zeros(length(wDsing),1);
state11tf = zeros(length(wDsing),1);
state12tf = zeros(length(wDsing),1);
state13tf = zeros(length(wDsing),1);
state21tf = zeros(length(wDsing),1);
state22tf = zeros(length(wDsing),1);
state23tf = zeros(length(wDsing),1);
state31tf = zeros(length(wDsing),1);
state32tf = zeros(length(wDsing),1);
state33tf = zeros(length(wDsing),1);

%% Singular arc :
for i = 1 : length(wDsing)
  S(i) = -(param.b(5)*param.wT*param.c10^2 + param.b(4)*param.wT*param.c10*...
          param.c20 + param.b(6)*param.wT*param.c20^2 - wD(i) - param.b(1)*param.wT);
end

%% Simulation :

for i = 1 : length(wDsing)
        x01m      = [param.c10;param.c20;0];
        options1m = odeset('Events',@event1m,'AbsTol',AbsTol,'RelTol',RelTol); 
        tspan1m   = [0 param.N];
        [time1m,state1m,te1m] = ode45(@modelstart,tspan1m,x01m,options1m,param,wD(i),S(i));
 
        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,wD(i),S(i));
 
        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,wD(i),S(i));
        
        tff(i)        = te3m;
        Vw(i)         = state3m(end,3);
        time1tf(i)    = time1m(end);
        time2tf(i)    = time2m(end);
        time3tf(i)    = time3m(end);
        state11tf(i)  = state1m(end,1);
        state12tf(i)  = state1m(end,2);
        state13tf(i)  = state1m(end,3);
        state21tf(i)  = state2m(end,1);
        state22tf(i)  = state2m(end,2);
        state23tf(i)  = state2m(end,3);
        state31tf(i)  = state3m(end,1);
        state32tf(i)  = state3m(end,2);
        state33tf(i)  = state3m(end,3);
end

x01m01      = [param.c10;param.c20;0];
options1m01 = odeset('Events',@event1m01,'AbsTol',AbsTol,'RelTol',RelTol);
tspan1m01   = [0 param.N];
[time1m01,state1m01,te1m01] = ode45(@model,tspan1m01,x01m01,options1m01,alpha0,param,wD(i),S(i));           

options2m01 = odeset('Events',@event2m01,'AbsTol',AbsTol,'RelTol',RelTol);
x02m01      = state1m01(end,:);
tspan2m01   = [time1m01(end) param.N];
[time2m01,state2m01,te2m01] = ode45(@model,tspan2m01,x02m01,options2m01,alpha1,param,wD(i),S(i));
        
%% Pareto Front :

% economically optimal : 
tfm = 2.398;
Vwm = 0.0352;

% minimum time : 
tft = 2.3451;
Vwt = 0.0675;

% minimum diluant :
tfd = 2.4244;
Vwd = 0.0315;

% control [10, as, 0, 1] --> [as, 0, 1]
tfs = 2.3666;
Vws = 0.0449;

figure
line        = 1.5;
marker_size = 8;
plot(tft,Vwt,'bs','markersize',marker_size)                                % minimum time
hold on
plot(tfd,Vwd,'bo','markersize',marker_size)                                % minimum diluant
plot(tfm,Vwm,'bd','markersize',marker_size)                                % multi-objective optimal
plot(tfs,Vws,'bh','markersize',marker_size)                                % switching strategy [10, as, 0, 1] --> [as, 0, 1]
plot(tff(1),Vw(end),'b*','markersize',marker_size)
plot(tff,Vw,'linewidth',line)
plot([tff(1);tff(1)],[Vw(1);Vw(end)],'b--','linewidth',line)                                 % J = tf
plot([tff(1);tff(end)],[Vw(end);Vw(end)],'b--','linewidth',line)                 % J = Vw
legend('minimum time','minimum diluant','economically optimal',sprintf('optimal control structure\nswitching point'),'utopia point')        % optimal control structure switching point
xlabel('tf')
ylabel('Vw')
axis([2.341 2.429 0.029 0.07])
%fs=6; fz=11;
%exportfig(gcf,'AlbEth_case9_PF.eps','width',fs,'fontmode','fixed','fontsize',fz,'Color','rgb');
set(gcf, 'PaperUnits', 'centimeters') 
set(gcf, 'PaperPosition', [3 10 11 8.5])
print -depsc 'AlbEth_case9_PF.eps'

function fun = model(t,x,u,param,wD,S)
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 = modelstart(t,x,param,wD,S)
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 S > 0
      alpha = 0;
elseif S < 0
      alpha = 1; 
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 fun = modelsing(t,x,param,wD,S)
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] = event1m(t,y,param,wD,S)
  c1opt = (((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)*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));
  value = y(1) - c1opt;   
  isterminal = 1;
  direction  = 0;
   function [value,isterminal,direction] = event1m01(t,y,u,param,wD,S)
value      = y(1) - param.c1f;
isterminal = 1;
direction  = 0;function [value,isterminal,direction] = event2m(t,y,param,wD,S)
value      = y(1) - param.c1f;
isterminal = 1;
direction  = 1;
function [value,isterminal,direction] = event2m01(t,y,u,param,wD,S)
value      = y(2) - param.c2f;
isterminal = 1;
direction  = 0;function [value,isterminal,direction] = event3m(t,y,u,param,wD,S)
value      = y(2) - param.c2f;
isterminal = 1;
direction  = -1;