% 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