genetic.m

基本遗传算法程序包

function [xopt,stats,options,bestf,fgen,lgen] = genetic(fun, ...
          x0,options,vlb,vub,bits,P1,P2,P3,P4,P5,P6,P7P,P8,P9,P10)
%GENETIC tries to maximize a function using a simple genetic algorithm.
%	X=GENETIC('FUN',X0,OPTIONS,VLB,VUB) uses a simple (haploid) 
%       genetic algorithm to find a maximum of the fitness function 
%       FUN (usually an M-file: FUN.M).  The user may define all or
%       part of an initial population X0 (or supply an empty argument
%       in which case an initial population will be chosen randomly
%       between the lower and upper bounds VLB and VUB.  Use OPTIONS
%       to specify optional parameters such as population size and
%       maximum number of generations produced.  Type HELP GOPTIONS
%       for more information.
%
%       The default algorithm uses a fixed population size, OPTIONS(11),
%       and no generational overlap.  Three genetic operations: 
%       reproduction, crossover, and mutation are performed during 
%       procreation.  The probability that an individual of the 
%       population will reproduce is proportional to its fitness.
%       Individuals chosen for reproduction are mated at random.
%       Mating produces two offspring (re: constant population size.)
%       Crossover in mating occurs with probability Pc=OPTIONS(12)
%       and the crossover index is randomly selected.  Each feature
%       of the offspring can mutate independently with probability
%       Pm=OPTIONS(13).  Default options are OPTIONS(11:13)=[30 1 0].
%       The default maximum generations OPTIONS(14) is 100.
%
%	X=GENETIC('FUN',X0,OPTIONS,VLB,VUB,BITS) allows the user to 
%	define the number of BITS used to code non-binary parameters
%	as binary strings.  Note: length(BITS) must equal length(VLB).
%
%	X=GENETIC('FUN',X0,OPTIONS,VLB,VUB,BITS,P1,P2,...) allows up 
%	to ten arguments, P1, P2, ... to be passed directly to FUN.
%	F=FUN(X,P1,P2,...).
%
%	[X,STATS,OPTIONS,BESTF,FGEN,LGEN]=GENETIC(<ARGS>)
%          STATS   - [max min mean std] for each generation
%          OPTIONS - options used
%          BESTF   - Fitness of indivadual X (i.e.: best fitness)
%          FGEN    - first generation population
%          LGEN    - last generation population
%
%       The algorithm implemented here is taken from Genetic
%       Algorithms in Search, Optimization, and Machine Learning,
%       David E. Goldberg, Addison-Wiley Publishing Company, Inc.,
%       1989.
%
%	Copyright (c) 1993 by the MathWorks, Inc.
%	Andrew Potvin 1-10-93.

% OPT_STOP used by user to halt optimization prematurely
% OPT_STEP will be true during the evaluation of the last individual's
% cost function.  Often used to determine when to update graphics.
global OPT_STOP OPT_STEP
OPT_STOP = 0;

% Argument and error checking
if nargin<4,
   error('No population bounds given.')
elseif (size(vlb,1)~=1) | (size(vub,1)~=1),
   % Remark: this will change if algorithm accomodates matrix variables
   error('VLB and VUB must be row vectors')
elseif (size(vlb,2)~=size(vub,2)),
   error('VLB and VUB must have the same number of columns.')
elseif (size(vub,2)~=size(x0,2)) & (size(x0,1)>0),
   error('X0 must all have the same number of columns as VLB and VUB.')
elseif any(vlb>vub),
   error('Some lower bounds greater than upper bounds')
else
   x0_row = size(x0,1);
   for i=1:x0_row,
      if any(x0(x0_row,:)<vlb) | any(x0(x0_row,:)>vub),
         error('Some initial population not within bounds.')
      end % if initial pop not within bounds
   end % for initial pop
end % if nargin<4   

if nargin<6,
   bits = [];
elseif (size(bits,1)~=1) | (size(bits,2)~=size(vlb,2)),
   % Remark: this will change if algorithm accomodates matrix variables
   error('BITS must have one row and length(VLB) columns')
elseif any(bits~=round(bits)) | any(bits<1),
   error('BITS must be a vector of integers >0')
end % if nargin<6

% Form string to call for function evaluation
if ~( any(fun<48) | any(fun>122) | any((fun>90) & (fun<97)) | ...
    any((fun>57) & (fun<65)) ), 
   % Only alphanumeric implies must be a function
   evalstr = [fun '(x'];
   for i=1:nargin-6,
      evalstr = [evalstr,',P',int2str(i)];
   end
   evalstr = [evalstr, ')'];
else
   evalstr = fun;
end

% Determine all options
% Remark: add another options index for type of termination criterion
if size(options,1)>1,
   error('OPTIONS must be a row vector')
else
   options = foptions(options);
   if options(11)==0,
      % Default size_pop
      options(11) = 30;
   end
   if options(12)==0,
      % Default Pc
      options(12) = 1;
   end
   if options(14)==0,
      % Default max_gen
      options(14) = 100;
   end
end
PRINTING = options(1);
terminate = options(2);
size_pop = options(11);
Pc = options(12);
Pm = options(13);
max_gen  = options(14);
% Ensure valid options: e.q. Pc,Pm,size_pop,max_gen>0, Pc,Pm<1
if any([Pc Pm size_pop max_gen]<0) | any([Pc Pm]>1),
   error('Some Pc,Pm,size_pop,max_gen<0 or Pc,Pm>1')
end

% Encode fitness (cost) function if necessary
ENCODED = any(any(([vlb; vub; x0]~=0) & ([vlb; vub; x0]~=1))) |  ...
          ~isempty(bits);
if ENCODED,
   [fgen,lchrom] = encode(x0,vlb,vub,bits);
else
   fgen = x0;
   lchrom = size(vlb,2);
end

% Display warning if odd number in initial population
if rem(size_pop,2)==1,
   disp('Warning: pop_size should be even.  Adding 1 to population.')
   size_pop = size_pop +1;
end

% Form random initial population if not enough supplied by user
if size(fgen,1)<size_pop,
   fgen = [fgen; (rand(size_pop-size(fgen,1),lchrom)<0.5)];
end
xopt = vlb;
bestf = -Inf;
new_gen = fgen;

% Header display
if PRINTING>=1,
   if ENCODED,
      disp('Cost function encoding as binary successful.')
      disp('')
      fgen = decode(fgen,vlb,vub,bits);
   end
   disp('                   Fitness statistics')
   disp('Generation Maximum      Minimum         Mean    Std. dev.')
end

% Set up main loop
STOP_FLAG = 0;
for generation = 1:max_gen+1,
   options(10) = generation-1;
   old_gen = new_gen;

   % Get fitness of each string in population
   % Decode first if necessary
   if ENCODED,
      x_pop = decode(old_gen,vlb,vub,bits);
   else
      x_pop = old_gen;
   end
   for i=1:size_pop,
      x = x_pop(i,:);
      if i==size_pop,
         OPT_STEP = 1;
      else
         OPT_STEP = 0;
      end
      fitness(i) = eval(evalstr);
   end
   [max_fit,INDEX] = max(fitness);
   stats = [stats; max_fit min(fitness) mean(fitness) std(fitness)];
   if max_fit>bestf,
      bestf = max_fit;
      xopt = x_pop(INDEX(1),:);
   else
      % Remark: may want to regenerate to guarantee cost decrease
      % Remark: be careful not to get stuck in infinite loop
   end

   % Display if necessary
   % Remark: consider alternate printing options
   if PRINTING>=1,
      disp([sprintf('%5.0f %12.6g %12.6g ',generation-1, ...
                     stats(generation,1),stats(generation,2)), ...
            sprintf('%12.6g %12.6g  ',stats(generation,3), ...
                    stats(generation,4))]);
   end

   % Check for termination
   % Remark: add more termination options
   if terminate>0,
      if generation>5,
         if ( (stats(generation,1)-stats(generation-5,1))/ ...
               stats(generation,1)<terminate ) & ...
            (stats(generation,1)>stats(generation-5,1)),
            STOP_FLAG = 1;
         end
      elseif generation>1,
         if ( (stats(generation,3)-stats(generation-1,3))/ ...
               stats(generation,3)<terminate ) & ...
            (stats(generation,3)>stats(generation-1,3)),
            STOP_FLAG = 1;
         end
      end
   end
   if STOP_FLAG | OPT_STOP,
      fprintf('\n')
      if STOP_FLAG,
         disp('Genetic algorithm converged.')   
      else
         disp('Genetic algorithm terminated by user.')
      end
      return
   end

   % Reproduce
   new_gen = reproduc(old_gen,fitness);

   % Mate
   new_gen = mate(new_gen);

   % Crossover
   new_gen = xover(new_gen,Pc);

   % Mutate
   new_gen = mutate(new_gen,Pm);
   if ENCODED,
      lgen = decode(new_gen,vlb,vub,bits);
   else
      lgen = new_gen;
   end

end % for max_gen

% Maximum number of generations reached without termination
if PRINTING>=1,
   fprintf('\n')
   disp('Maximum number of generations reached without termination')
   disp('criterion met.  Either increase maximum generations')
   disp('or ease termination citerion.')
end

% end genetic