Hi,
I'm using the genetic algorithm to minimize a function. However the value of the (fitness) function when the program is done is different than the fitness function value returned by the function "ga" (fval). For instance, ga(@modele_ga_time,nbparam,[],[],[],[],lb,ub,[],options) returns a fval of 0.043 where modele_ga_time is defined as function f = modele_ga_time(x), and the actual calculated value for "f" is 3.4742!! Any idea of why this discrepancy. Below are the files used (main: "inversion_ga_time.m", the fitness function: "modele_ga_time.m", the dataset used: "test.dat" and the function integrated in modele_ga_time.m "fonction.m" Thanks very much in advance!!
-------------------------------- -------------------------------- inversion_ga_time.m -------------------------------- -------------------------------- clear all
global nbcoeff model f nbsample nbratios ratios_calc ratios_cib l238 l234 l230 k238 k234 k230 U8i U4i Th0i F238 F234 F230 k234_8 F234_8 f_k238 k230_4 f_k238 f230_4
nbsol = 1; model = 6; choice = 1; % line chosen where to started picking samples in loaded data file nbsample = 6; nbratios = 2; % tolerror = 0.01; % error allowed between calculated and observed ratio % modele % 5: loss/gain of U isotopes only % 6: loss/gain of all nuclides % 7: loss of U isotopes only, no gain (no f values)
% load Bisley.dat; load test.dat; % target values for i = 1:nbsample, for j =1:nbratios, ratios_cib(i,j) = test(i+choice-1,j); % ratios_cib(i,j) = Bisley(i+choice-1,j); end end
% Initial ratios % r48i = 0.993; r04i = 0.814; r02i = 1; r48i = 1; r04i = 1; % r48i = Bisley(choice+1,1); r04i = Bisley(choice+1,2); r02i = Bisley(choice+1,4);
% Initial U concentration U_init = 0.2; % in ppm
% parameters: log(k238) log(k234/8) log(f/k238)log(f234/8) log(k230/4) log(f230/4) log(time(yr)) % lower bound values for parameters low = [-7 0 -2 0 -7 -7 2]; % upper bound values for parameters up = [-3 log10(2) +2 log10(2) +2 +2 6]; % lower bound values for initial pop lowi = [-5 0 +0 0 -6 -3 4]; % upper bound values for initial pop upi = [-4 log10(1.2) +1 log10(1.2) -4 +0 5];
if model == 5, % parameters: log(k238) log(k234/8) log(f/k238) log(f234/8) log(time(yr)) lb = [low(1:4) low(7)*ones(1,nbsample)]; % lower bound values for parameters ub = [up(1:4) up(7)*ones(1,nbsample)]; % upper bound values for parameters elseif model == 6, % parameters: log(k238) log(k234/8) log(f/k238) log(f234/8) log(k230/4) log(f230/4) log(time(yr)) lb = [low(1:6) low(7)*ones(1,nbsample)]; % lower bound values for parameters ub = [up(1:6) up(7)*ones(1,nbsample)]; % upper bound values for parameters lbi = [lowi(1:6) lowi(7)*ones(1,nbsample)]; % lower bound values for initial pop ubi = [upi(1:6) upi(7)*ones(1,nbsample)]; % upper bound values for initial pop elseif model == 7, % parameters: log(k238) log(k234/8) log(time(yr)) lb = [low(1:2) low(7)*ones(1,nbsample)]; % lower bound values for parameters ub = [up(1:2) up(7)*ones(1,nbsample)]; % upper bound values for parameters end
options = gaoptimset('PopInitRange',[lbi; ubi],'FitnessLimit',0.005,'Generations',500,'TolFun',1e-10,'TolCon','',... 'MigrationDirection','both','CrossoverFraction','','PopulationSize',nbsample*nbratios*2,... 'MutationFcn',@mutationadaptfeasible,'CrossoverFcn',@crossoverheuristic,'PlotFcns',{@gaplotbestf,@gaplotbestindiv,@gaplotstopping,@gaplotscorediversity}); nbcoeff = length(lb) - nbsample;
% decay constants (yr-1) l238 = 0.1551e-9; l234 = 2.826e-6; l230 = 9.158e-6;
% intial abundances (nbr atoms) U8i = U_init*1e-6/238*6.02e23; U4i = r48i*l238/l234*U8i; Th0i = r04i*l234/l230*U4i;
nbparam = length(lb);
j=0; while j<nbsol, [sol fval exitflag output population] = ga(@modele_ga_time,nbparam,[],[],[],[],lb,ub,[],options); j=j+1; disp(j); for k = 1:nbparam, msol(j,k) = sol(k); end mfunc(j) = fval; disp(10^msol(j,nbcoeff+1)); disp(mfunc(j)); func(j) = (f); % Th230_232_i(j) = r02i; mratios(:,:,j) = ratios_calc(:,:); end
---------------------------------------------------------- ---------------------------------------------------------- modele_ga_time.m ------------------------------------------------------- -------------------------------- function f = modele_ga_time(x)
global nbcoeff model f nbsample nbratios ratios_calc ratios_cib l238 l234 l230 k238 k234 k230 U8i U4i Th0i F238 F234 F230 k234_8 F234_8 f_k238 k230_4 f_k238 f230_4
if model == 5, k238 = 10^x(1); k234_8 = 10^x(2); k230_4 = 1e-6; f_k238 = 10^x(3); f234_8 = 10^x(4); f230_4 = 1e-6; elseif model == 6, k238 = 10^x(1); k234_8 = 10^x(2); k230_4 = 10^x(5); f_k238 = 10^x(3); f234_8 = 10^x(4); f230_4 = 10^x(6); elseif model == 7, k238 = 10^x(1); k234_8 = 10^x(2); k230_4 = 1e-6; f_k238 = 1e-10*k238; f234_8 = 1; f230_4 = 1e-6; end
k234 = k234_8*k238; k230 = k230_4*k234; f238 = f_k238*k238; F238 = f238*l238*U8i; f234 = f234_8*f238; F234 = f234*l234*U4i; f230 = f230_4*f234; F230 = f230*l230*Th0i;
for i = 1:nbsample, [t N] = ode23(@fonction,[0 10^x(i+nbcoeff)], [U8i U4i Th0i]); nbtime = length(t); ratios_calc(i,:) = [l234/l238*N(nbtime,2)./N(nbtime,1) l230/l234*N(nbtime,3)./N(nbtime,2)]; end
f = sum(sum((ratios_calc - ratios_cib).^2));
--------------------------------------------------- --------------------------------------------------- fonction.m ---------------------------------------------------- -------------------------------- function yp = fonction(t,N)
global l238 l234 l230 k238 k234 k230 F238 F234 F230
% N(1) = U238; N(2) = U234; N(3) = Th230;
yp = [-k238*N(1) + F238/l238; l238*N(1) - (l234 + k234)*N(2) + F234/l234; l234*N(2) - (l230 + k230)*N(3) + F230/l230];
--------------------------------------------------- --------------------------------------------------- test.dat ---------------------------------------------------- ---------------------------------------------------- % test dataset % k238 = 1e-5; k234_8 = 1.2; k232_8 = 1e-6; k230_4 = 1e-6; f_k238 = 10; F234_8 = 1.5; f232_8 = 2; f230_2 = 1;
1.018 1.024 400.000 1.093 1.115 2400.000 1.202 1.248 7400.000 1.259 1.330 12400.000 1.292 1.392 17400.000 1.313 1.444 22400.000 1.325 1.490 27400.000 1.333 1.533 32400.000 1.337 1.572 37400.000 1.340 1.610 42400.000 1.340 1.646 47400.000 1.340 1.664 50000.000 ------------------------------
