gareal_tourney.m

一个启发式搜索算法 给初学者用的 随便参考着用 提供了多种策略

% J. C. Spall, October 1999
% Written in support of text, Introduction to Stochastic Search and Optimization, 2003
% This program runs a GA with real-number coding.  Elitism is used 
% and the mutation operator is simply the addition of a Gaussian
% random vector to the non-elite elements.
% 
% The user is expected to set a variable 'expect_fn' representing the 
% expected number of function evaluations allowed.  The main "while"
% loop in the code tests for when the expected number of function
% evaluations exceeds 'expect_fn' and terminates the iteration on the 
% first generation that would require a no. of function evaluations
% exceeding this value.  The while loop expects that the mutation 
% operator is invoked (P_m represents mutation stand. deviation here; 
% P_m n.e. 0), so that all non-elite elements are
% guaranteed to change at every iteration (differs from the bit form
% where it is possible to repeat a loss function call within and across
% iterations).
%  
% This code handles noisy loss functions via the fact that the variables
% 'loss' and 'lossfinal' can call different functions.  The global 
% variable 'sigma' is used to generate the scaling for noise in the 
% loss measurements.
%
global p N thetamin thetamax
cases=20;         		%no. of Monte Carlo cases  
N=80;          			%population size
p=2;            			%dimension of theta
P_c=.8;         			%crossover rate
P_m=0.05;      			%mutation sigma
elite=2;                %no. of population elements for "elitism" 
                        %(should pick s.t. N-elite = even no.)
expect_fn=4000;         	%expected no. of function evaluations                        
thetamin=zeros(p,1);  		%lower bounds on theta elements
thetamax=10*ones(p,1);   	%upper bounds: ditto
thetamin_0=zeros(p,1);   	%lower bounds on theta for initialization 
thetamax_0=10*ones(p,1);   %upper bounds: ditto 
lossfinaleval ='loss_example9_5';  	%choice of loss function for final perf.
                              	%evaluation (no noise)
loss='loss_example9_5';    		%loss function for use in algorithm evaluations
lossfinalsq=0;          %variable for cum.(over 'cases')squared loss values
lossfinal=0;            %variable for cum. loss values
rand('seed',31415927);
randn('seed',3111113)
global z sigma;       	%declaration of random var. used for normal noise
                      	%generation in loss fn. calls given seed above;
                      	%also sigma in noise (noise is dependent on theta)
sigma=0;            		%multiplier in loss noise (lower bd. to s.d.)          
fitness=zeros(N,1);                %dim.-setting statement
theta=zeros(p,1);                  %dim.-setting statement
%
thetapop=zeros(N,p);  %dummy statement setting dimensions of matrix
							 %containing the theta values  
thetapop_noelite=zeros(N-elite,p);%dummy statement setting
                                  %dimension of matrix used to temp.
                                  %store offspring  
parent1=zeros(1,p); %dummy dimension statement
parent2=zeros(1,p); %ditto 
% User warning if P_m = 0 (to ensure that proper number of function evaluations
% are taken; code assumes all non-elite elements change every iteration)
if P_m==0
   disp('WARNING: P_m should be > 0');
else
end

% Begin outer loop of 'cases' Monte Carlo runs
%
navg=0;
for c=1:cases
   c
   %
   % Initial Random Population
   %
   % Statements below allow intialization of the search using a random
   % placement in natural theta space.  Points are placed uniformly dist.
   % on the hypercube defined by thetamax_0,thetamin_0.  
   %
   for i=1:N
      thetapop(i,:)=((thetamax_0-thetamin_0).*rand(p,1)+thetamin_0)';
   end
   % Evaluation of fitness for initial population
   for i=1:N
      fitness(i)=-feval(loss,thetapop(i,:)');%N by 1 vector of fitness values
   end
   [bestfitness,bestindex]=max(fitness);  bestindex
   %
   % Main loop of GA
   %
   index=zeros(elite,1);
   %cumfit=zeros(N-elite,1);
   n=1;               
   while n<=(expect_fn-N)/(N-elite)
      % Invoke elitism for 'elite' best chromosomes
      temp_fit=fitness;
      for j=1:elite
         [junk,index(j)]=max(temp_fit);
         temp_fit(index(j))=min(temp_fit);
      end
      for j=1:2:N-elite-1
         % Tournament selection below.  Uses tournament size = 2.  
            cand=ceil(N*rand(2,1));  %picks two candidate parents
                                     %from which one will be chosen   
         if fitness(cand(1))>fitness(cand(2))
            parent1_idx=cand(1);
         else
            parent1_idx=cand(2);
         end   
         cand=ceil(N*rand(2,1));  %ditto-for picking other parent   
         if fitness(cand(1))>fitness(cand(2))
            parent2_idx=cand(1);
         else
            parent2_idx=cand(2);
         end
         % Crossover with above two parents
         if rand < P_c
            cross_point=ceil(rand*(p-1)); %an integer 1,...,p-1
            parent1=thetapop(parent1_idx,:);
            %above picks off the chrom. for parent1_idx
            parent2=thetapop(parent2_idx,:);
            %above picks off the chrom. for parent2_idx
            child1=[parent1(1:cross_point),...  
               parent2(cross_point+1:p)];   %merges parents1&2
            child2=[parent2(1:cross_point),...
               parent1(cross_point+1:p)];   %merges parents1&2
            thetapop_noelite(j,:)=child1;%assigns child1 to new chrom. matrix
            thetapop_noelite(j+1,:)=child2;%ditto for child2 
         else
            parent1=thetapop(parent1_idx,:);
               %above picks off the chrom. for parent1_idx
            parent2=thetapop(parent2_idx,:);
               %above picks off the chrom. for parent2_idx
            thetapop_noelite(j,:)=parent1;%assigns parent1 to chrom. matrix
            thetapop_noelite(j+1,:)=parent2;%ditto for parent2 
         end
      end
      %
      % Mutation operator applied to thetapop_noelite vector (elite  
      % chromosomes not subject to mutation)by addition of Gauss.
      % random vector with diag. cov. matrix (standard devs. = P_m)
      %
      thetapop_noelite=thetapop_noelite+P_m*randn(N-elite,p);
      %
      % Invoke constraints to above population
      %
      for j=1:N-elite
         thetapop_noelite(j,:)=min(thetapop_noelite(j,:),thetamax');
         thetapop_noelite(j,:)=max(thetapop_noelite(j,:),thetamin');
      end
      %
      % Put elite chromosomes (theta values) back in population
      %
      for j=1:elite
        thetapop(j,:)=thetapop(index(j),:);
      end
      thetapop(elite+1:N,:)=thetapop_noelite;
      for i=1:N
         fitness(i)=-feval(loss,thetapop(i,:)');
      end
      n=n+1;
   end
   [bestfitness,bestindex]=max(fitness);
   lossvalue=feval(lossfinaleval,thetapop(bestindex,:)');
   lossfinalsq=lossfinalsq+lossvalue^2;
   lossfinal=lossfinal+lossvalue;
   navg=navg+n-1;				%counter for # of iters. taken
end
navg/cases
% standard dev. of loss values
if cases ~= 1 
	disp('standard deviation of mean loss value') 
   sd=((cases/(cases-1))^.5)*(lossfinalsq/cases-(lossfinal/cases)^2)^.5;
   sd=sd/(cases^.5)
else
end
% normalized loss values
disp('mean loss value over "cases" runs') 
mean=lossfinal/cases