function parfront
% c1 = macro-solute (pectin)
% c2 = micro-solute (glucose+fructose+sacharose)

% Type of problem: (Parameter "N")
% 1 --> minimum time
% 2 --> minimum diluant
% 3 --> multi-objective optimal

x0 = [0.15  , 120];   c10 = x0(1); c20 = x0(2);        % c0 [kg/m3]
xf = [11.913, 0.1];   c1f = xf(1); c2f = xf(2);        % cf [kg/m3]
V0 = 100;                                              % m3                          

% parameters : 
clim    = 32.383;                                      % kg/m3
k       = 8.34e-6*3600;                                % m/s
A       = 105.6;                                       % m2  
epsilon = 1.02;                                        % eps = 2%

% Weight coefficients : 
nn      = 1000;
wT      = 1;                                           % -
wD      = logspace(-6,6,nn);                          % -

% Calculating optimal time and water consumption for different values of weight coefficients
for i=1:length(wD) 
  copt(i) = clim*exp((wT - (wT^2 + 4*k*A*wT*wD(i))^(1/2))/(2*A*k*wD(i)));
  [tf(i),Vw(i),c11,c12,c22,dt1,dt2,dt3]=find_tf_vw(wT,wD(i),c10,c20,c1f,c2f,V0,clim,k,A,3);
  str1 = fprintf('Step: %d/%d, wTwd: %d \n', i,nn,wT/wD(i));
end

% calculating time and water consumption (minimum time problem): 
wTt = 1;
wDt = 0;
[tft,Vwt,c11t,c12t,c22t,dt1t,dt2t,dt3t]=find_tf_vw(wTt,wDt,c10,c20,c1f,c2f,V0,clim,k,A,1);

% calculating time and water consumption (minimum diluant problem): 
wTd = 0;
wDd = 1;
[tfd,Vwd,c11d,c12d,c22d,dt1d,dt2d,dt3d]=find_tf_vw(wTd,wDd,c10,c20,c1f,c2f,V0,clim,k,A,2);

% Calculating multi-objective optimal point :
wTm = 3.6;      % euro/h 
wDm = 10;       % euro/m3 
[tfm,Vwm,c11m,c12m,c22m,dt1m,dt2m,dt3m]=find_tf_vw(wTm,wDm,c10,c20,c1f,c2f,V0,clim,k,A,3);

% Calculating almost minimum time :
wTm1 = 1;      % euro/h 
wDm1 = 0.1;       % euro/m3 
[tfm1,Vwm1,c11m1,c12m1,c22m1,dt1m1,dt2m1,dt3m1]=find_tf_vw(wTm1,wDm1,c10,c20,c1f,c2f,V0,clim,k,A,3);

% Calculating real minimum diluant :
wTd1 = 0.01;      % euro/h 
wDd1 = 1;       % euro/m3 
[tfd1,Vwd1,c11d1,c12d1,c22d1,dt1d1,dt2d1,dt3d1]=find_tf_vw(wTd1,wDd1,c10,c20,c1f,c2f,V0,clim,k,A,3);

marker_size = 7;
% pareto-front
figure
plot(tft, Vwt,'bs','markersize',marker_size)                            % minimum time
hold on
plot(tfd1,Vwd1,'bo','markersize',marker_size)                             % minimum diluant
plot(tfm,Vwm,'bd','markersize',marker_size)                                % multi-objective optimal
plot(tfm1, Vwm1,'bx','markersize',marker_size)                            % almost min time 
plot(tft,Vwd,'b*','markersize',marker_size)                  % utopia point
plot(tf,Vw)
plot([tft;tft],[Vw(1);Vw(end)],'b--')                                  % J = tf
plot([tf(1);tf(end)],[Vwd;Vwd],'b--')              % J = Vw
%%% plot(tf(1),Vw(end)*epsilon,'b*')
legend('minimum time','minimum diluant','economically optimal','almost min time','utopia point')
axis([7.3 42.23 3.0 9])
xlabel('tf')
ylabel('Vw')
set(gcf, 'PaperUnits', 'centimeters') 
set(gcf, 'PaperPosition', [3 10 11 8.25])
print -depsc 'PF_applejuice_1.eps'
clc
disp('==========================================================================')
str2 = fprintf('Minimum time problem:    tf = %d [h], Vw = %d [m3]\n',tft,Vwt);
disp('==========================================================================')
str3 = fprintf('Minimum diluant problem: tf = %d [h], Vw = %d [m3]\n',tfd,Vwd);
disp('==========================================================================')
str4 = fprintf('Multi-objective optimal: tf = %d [h], Vw = %d [m3]\n',tfm,Vwm);
disp('==========================================================================')
str5 = fprintf('wt/wd=10: tf = %d [h], Vw = %d [m3]\n',tfm1,Vwm1);
disp('==========================================================================')
str6 = fprintf('wt/wd=0.01: tf = %d [h], Vw = %d [m3]\n',tfd1,Vwd1);

function [tf,Vw,c11,c12,c22,delt1,delt2,delt3]=find_tf_vw(wT,wD,c10,c20,c1f,c2f,V0,clim,k,A,N)

if (N == 1)
  copt = clim/exp(1);
elseif (N == 2)
  copt = clim;
else
  copt = clim*exp((wT - (wT^2 + 4*k*A*wT*wD)^(1/2))/(2*A*k*wD));
end

% ---------------------------------------------------------------------------------

     % 1st phase (c2 = const, c10 -> c11) pure ultrafiltration (alpha = 0)
     % until c1 = copt
     
     c11   = copt;
     delt1 = ((c10*V0)/(k*A*clim))*(mfun('Li',clim/c10) - mfun('Li',clim/c11));
     if (N==2)
       delt1 = Inf;
     end
     Vw1   = 0;

% ---------------------------------------------------------------------------------

     % 2nd phase (c1 = const, c12 -> c11) diafiltration
     c12 = copt;

     % CVD step until : c1/c2 = c1f/c2f  (if c1f < copt(c11->c12))
     if (c1f < c12)  
       c22   = c12*c2f/c1f;
       if (N ~=2) 
	 delt2 = ((c10*V0)/(k*A*c12))*((log((c20*c1f)/(c12*c2f)))/(log(clim/c12))); 
	 Vw2   = (k*A*log(clim/c12))*delt2;
       else
	 delt2 = Inf;
	 Vw2 = (c10*V0/clim)*(log((c20*c1f)/(clim*c2f)));
       end
     end
     
     % CVD step until c2 = c2f (if c1f > copt)
     if (c1f > c12)
       c22   = c2f;
       if (N ~=2) 
	 delt2 = ((c10*V0)/(k*A*c12))*(log(c20/c2f)/log(clim/c12)); 
	 Vw2   = (k*A*log(clim/c12))*delt2;
       else
	 delt2 = Inf;
	 Vw2 = (c10*V0/clim)*log(c20/c2f); 
       end
     end

% ---------------------------------------------------------------------------------

     % 3rd phase
     c13 = copt;

     % pure dilution (alpha = \infty) when c1f < copt
     if (c1f < c13)
        delt3 = 0;
        Vw3   = c10*V0*(1/c1f - 1/c13);
     end
     
     % concentration step (ultrafiltration --> alpha = 0) when c1f > copt
     if (c1f > c13)
        delt3 = ((c10*V0)/(k*A*clim))*(mfun('Li',clim/c13) - mfun('Li',clim/c1f));
        Vw3   = 0;
     end

% ---------------------------------------------------------------------------------

     % final results : 
     tf = delt1 + delt2 + delt3;
     Vw = Vw1 + Vw2 + Vw3;
%     [wT, wD, delt1 delt2 delt3 Vw1 Vw2 Vw3 ]
