genetic algorithm does not give actual calculated value of provided fitness function

 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
------------------------------

[Seth DeLand]

Hi Tony,

I've ran your code a few times now, and I'm not able to reproduce the inconsistency. The value for fval seems to always be equal to modele_ga_time(sol). If you're still getting different values for fval and modele_ga_time(sol), I would check to make sure there aren't any issues with the global variables. Because the objective function is dependent on the globals, changing a value of a global variable will change the value of the objective function.

One other thing to watch out for: the very last time that the Genetic Algorithm evaluates your objective function, it is highly likely that this is not at the best point that it has found. If there are any global variables that are set during this final evaluation, and those values aren't set to the correct value the next time you call modele_ga_time (to verify the result), that would be an issue you need to address.

In general, I try to avoid using global variables to pass extra parameters to GA. There's a good section in the documentation that provides some alternatives that are less risky: http://www.mathworks.com/help/toolbox/optim/ug/brhkghv-7.html

Answer 1

[mehdi]

Hi Tony, I have same problem. I'm using GA in Global Optimization Toolbox to solve a MINLP. The optimal value for fitness function by GA is different by the one that I found using optimal decision variables in the model..sometimes they are very close and sometimes the gap is significant. However, when there is no integer decision variable in the optimization, I have no problem and the optimal value from GA is the same as the one for the model using optimal decision variables.

Did you find any answer for your problem, if yes, would you please let me know about it? Regards, Mehdi