function parfront
x0=[3.3, 5.5]; c10=x0(1); c20=x0(2);
xf=[9.04, 0.6]; c1f=xf(1); c2f=xf(2);
V0=100;

b0=63.42; b1=-12.439; b2=-7.836;
% optimal strategy alpha={0,1,\infty}

WT=logspace(-5,5,50);
% WT = 0.02537;
for i=1:length(WT) 
    [tf(i),Vw(i)]=find_tf_vw(WT(i),c10,c20,c1f,c2f,b0,b1,b2,V0);
end

line        = 1.5;
marker_size = 9;
epsilon     = 1.0025;

figure
plot(4.531,Vw(end),'bs','markersize',marker_size)       % minimum time
hold on
plot(8.91,32.8,'bo','markersize',marker_size)           % minimum diluant
plot(4.756,36.6,'bd','markersize',marker_size)          % minimum multi-objective
plot(4.5309, 32.7555*epsilon,'*','markersize',marker_size)
plot(tf,Vw,'linewidth',line)
plot([4.5309,4.5309],[32.7555,Vw(end)],'b--','linewidth',line)                                % J = tf
plot([4.5309*epsilon,tf(1)*epsilon],[32.7555*epsilon,32.7555*epsilon],'b--','linewidth',line) % J = Vw
axis([4.3 9.05 32.4 48.9])
xlabel('J1')
ylabel('J2')
legend('minimum time','minimum diluant','multi-objective optimal','utopia point')
% set(xy,'Fontsize',30);
% % print(gcf, '-depsc2', 'PF_SOLP2')
% opts = struct('FontMode','fixed','FontSize',12);
% exportfig(gcf, 'PF_SOLP2.eps', opts);
% fs=6; fz=11;
% exportfig(gcf,'PF_SOLP.eps','width',fs,'fontmode','fixed','fontsize',fz,'Color','rgb');
set(gcf, 'PaperUnits', 'centimeters') 
set(gcf, 'PaperPosition', [3 10 11 8.25])
print -depsc 'parfront_lactprot.eps'


 c1=logspace(log10(10),log10(70),30);
 figure(2)
 
 wT=1; wD=eps; WT=wT/wD;
 Sc2t=exp(-(b0 + (wT - (-wT*(4*b1*wD - wT + 4*b2*wD))^(1/2))/(2*wD) + b1*log(c1))/b2);
 
 wT=0; wD=1; WT=wT/wD;
 Sc2d=exp(-(b0 + (wT - (-wT*(4*b1*wD - wT + 4*b2*wD))^(1/2))/(2*wD) + b1*log(c1))/b2);
 
 wT=0.0132; wD=0.001; WT=wT/wD;
 Sc2m=exp(-(b0 + (wT - (-wT*(4*b1*wD - wT + 4*b2*wD))^(1/2))/(2*wD) + b1*log(c1))/b2);
 
 wT=1; wD=eps; WT=wT/wD;
 [tft,Vwt,dt1t,dt2t,c11t,c12t,c22t]=find_tf_vw(WT,c10,c20,c1f,c2f,b0,b1,b2,V0);
 
 wT=0; wD=1; WT=wT/wD;
 [tfd,Vwd,dt1d,dt2d,c11d,c12d,c22d]=find_tf_vw(WT,c10,c20,c1f,c2f,b0,b1,b2,V0);
 
 wT=0.0132; wD=0.001; WT=wT/wD;
 [tfm,Vwm,dt1m,dt2m,c11m,c12m,c22m]=find_tf_vw(WT,c10,c20,c1f,c2f,b0,b1,b2,V0);
 
 %% first phase 
 plot([c10;c11d],[c20;c20],'g--','linewidth',line)
 hold on
 plot([c10;c11m],[c20;c20],'r','linewidth',line)
 plot([c10;c11t],[c20;c20],'b-.','linewidth',line)
 legend('minimum diluant','multi-objective optimal','minimum time','Location','SouthEast')
 
 plot(c10,c20,'go','MarkerSize',marker_size)
 plot(c1f,c2f,'rx','MarkerSize',marker_size)
 
 
 %% third phase 
 plot([c1f;c12d],[c2f;c22d],'g--','linewidth',line)
 plot([c1f;c12m],[c2f;c22m],'r','linewidth',line)
 plot([c1f;c12t],[c2f;c22t],'b-.','linewidth',line)
 
 %% singular phase :
 c1t=linspace(c11t,c12t,10); wT=1; wD=eps;
 c2t=exp(-(b0 + (wT - (-wT*(4*b1*wD - wT + 4*b2*wD))^(1/2))/(2*wD) + b1*log(c1t))/b2);
 
 c1d=linspace(c11d,c12d,10); wT=0; wD=1;
 c2d=exp(-(b0 + (wT - (-wT*(4*b1*wD - wT + 4*b2*wD))^(1/2))/(2*wD) + b1*log(c1d))/b2);
 
 c1m=linspace(c11m,c12m,10); wT = 0.0132; wD = 0.001;
 c2m=exp(-(b0 + (wT - (-wT*(4*b1*wD - wT + 4*b2*wD))^(1/2))/(2*wD) + b1*log(c1m))/b2);
 
 plot(c1d,c2d,'g--','linewidth',line)
 plot(c1m,c2m,'r','linewidth',line)
 plot(c1t,c2t,'b-.','linewidth',line)
 
 plot(c1,Sc2t,'o',c1,Sc2d,'o',c1,Sc2m,'o');
 xlabel('c1'); ylabel('c2');
 axis([0 71 0 6])
% fs=6; fz=12;
% exportfig(gcf,'mulobj_mtim_mdil_lactprot.eps','width',fs,'fontmode','fixed','fontsize',fz,'Color','rgb');
set(gcf, 'PaperUnits', 'centimeters') 
set(gcf, 'PaperPosition', [3 10 16 10])
print -depsc 'mulobj_mtim_mdil_lactprot.eps'

dend = 5.3;
 droz = 8.91-dend;
 % graf riadenia : 
 figure(3)
 plot([0;dend],[0;0],'g--','linewidth',line)
 hold on
 plot([0;dt1m],[0;0],'r-','linewidth',line)
 plot([0;dt1t],[0;0],'b-.','linewidth',line)
 plot([dt1t;dt1t],[0;0.61],'b-.',[dt1t;dt1t+dt2t],[0.61;0.61],'b-.','linewidth',line)
 plot([dt1m;dt1m],[0;0.61],'r-',[dt1m;dt1m+dt2m],[0.61;0.61],'r-','linewidth',line)
 plot([dend;dend],[0;0.61],'g--',[dend;dend+droz],[0.61;0.61],'g--','linewidth',line)
 plot([dt1t+dt2t;dt1t+dt2t],[0.61;1.1],'b-.',dt1t+dt2t,1.1,'^','MarkerEdgeColor','b','MarkerFaceColor','b','linewidth',line,'MarkerSize',marker_size)
 plot([dt1m+dt2m;dt1m+dt2m],[0.61;1.1],'r',dt1m+dt2m,1.1,'^','MarkerEdgeColor','r','MarkerFaceColor','r','linewidth',line,'MarkerSize',marker_size)
 plot([dend+droz;dend+droz],[0.61;1.1],'g--',dend+droz,1.1,'^','MarkerEdgeColor','g','MarkerFaceColor','g','linewidth',line,'MarkerSize',marker_size)
 xlabel('time')
 ylabel('alpha')
 legend('minimum diluant','multi-objective optimal','minimum time','Location','SouthEast')
 axis([0 12.2 -0.05 1.13])
 set(gca, 'box','on');
% fs=6; fz=12;
% exportfig(gcf,'alpha_SOLP.eps','width',fs,'fontmode','fixed','fontsize',fz,'Color','rgb');
set(gcf, 'PaperUnits', 'centimeters') 
set(gcf, 'PaperPosition', [3 10 16 10])
print -depsc 'alpha_SOLP.eps'

function [tf,Vw,delt1,delt2,c11,c12,c22]=find_tf_vw(WT,c10,c20,c1f,c2f,b0,b1,b2,V0)

% singular surface S=wD*q^2+wT*(q+c1*q1+c2*q2);
% wT=1; wD=eps;
% wT=1e-7; wD=1;
wT=WT; wD=1;
c2=c20;
c11=exp(-(b0+(wT-(-wT*(4*b1*wD-wT+4*b2*wD))^(1/2))/(2*wD)+b2*log(c2))/b1);


% 1st phase (c2=const., c10 -> c11)
delt1=((c2^(b2/b1)*exp(b0/b1)*mfun('Li',exp(-b0/b1)/(c11*c20^(b2/b1))))/b1-...
    (c2^(b2/b1)*exp(b0/b1)*mfun('Li',exp(-b0/b1)/(c10*c20^(b2/b1))))/b1)*c10*V0;
Vw1=0;
V1=c10*V0/c11;

% 2nd phase (c11->c12, c20->c22, q=const., alpha_vvd)
q=(-wT/wD+sqrt(wT^2/wD^2-4*wT/wD*(b1+b2)))/2/wD;
% q=-wT+sqrt(wT^2-4*wT*wD*(b1+b2))/2/wD;
% q=-b1-b2
alpha=b1/(b1+b2);
c12=(c11^alpha*c20*c1f*c2f^alpha)/(c20^alpha*c1f^alpha*c2f);
c22=c12*c2f/c1f;
delt2=(-1/c12+1/c11)/q*c10*V0/(1-alpha);
Vw2=alpha*q*delt2;
V2=(alpha-1)*q*delt2+V1; % V2=c10*V0/c12

%3rd phase (c12->c1f, c22->c2f)
delt3=0;
Vw3=(c12/c1f-1)*V2;

% final results
tf=delt1+delt2+delt3;
Vw=Vw1+Vw2+Vw3;    
