%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Genetic Algorithm (GA)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Kevin Passino
% Version: 7/21/98, latest change 10/22/03 with correction from
% Phan Tran Ho Truc, Mr., B. Eng., lecturer in HCM City University of Technology
% to lines 379-384 of old program (lines 380-381 of this program)
%
% Notes: This program has evolved (hopefully to achieve a higher fitness!)
% over time using programming ideas from several persons including
% LaMoyne Porter, Will Lennon, Jonathan Cook, and Jim Gremling.
%
% This program simulates the optimzation of a simple function with
% a base-10 genetic algorithm. First, we outline some background
% information on GAs and define some terminology. Following this,
% we explain how to use the program for your own purposes and as an
% example we show how to solve an optimization problem.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Background Information on GAs:
%
% A genetic algorithm is an optimization method that
% manipulates a string of numbers in a manner similar to how chromosomes
% are changed in biological evolution. An initial population made up of
% strings of numbers is chosen at random or is specified by the
% user. Each string of numbers is called a "chromosome" or an
% "individual," and each number slot is called a "gene." A
% set of chromosomes forms a population. Each chromosome
% represents a given number of traits which are the actual
% parameters that are being varied to optimize the "fitness function".
% The fitness function is a performance index that we seek to maximize.
%
% The operation of the GA proceeds in steps. Beginning with the initial
% population, "selection" is used to choose which chromosomes
% should survive to form a "mating pool." Chromosomes are chosen based on
% how fit they are (as computed by the fitness function) relative
% to the other members of the population. More fit individuals end up
% with more copies of themselves in the mating pool so that they
% will more significantly effect the formation of the next
% generation. Next, several operations are taken on the mating
% pool. First, "crossover" (which represents mating, the exchange
% of genetic material) occurs between parents.
% To perform crossover, a random spot is picked
% in the chromosome, and the genes after this spot are switched
% with the corresponding genes of the other parent. Following this,
% "mutation" occurs. This is where some genes are randomly changed
% to other values. After the crossover and mutation operations
% occur, the resulting strings form the next generation
% and the process is repeated. A termination criterion
% is used to specify when the GA should end (e.g., the maximum
% number of generations or until the fitness stops increasing).
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Some GA Definitions for this Program:
%
% A GENE is a digit position that can take on certain values
% Ex: In a base-10 algorithm a gene can hold any number between 0 and 9.
%
% A CHROMOSOME is a string of genes.
% Ex: In base-10 1234567890 could be a chromosome of length 10.
%
% A TRAIT is a decimal number which is decoded from a chromosome.
% Normally, a chromosome is a concatenation of several TRAITS.
%
% An INDIVIDUAL is the object that the GA is attempting to optimize.
% An individual is described by its chromosome.
% The individual's traits determine its fitness.
%
% A POPULATION is a set of individuals (set of chromosomes).
%
% FITNESS FUNCTION is the objective function to be optimized which provides
% the mechanism for evaluating each string (we maximize the fitness
% function).
%
% SELECTION: a fitter string receives a higher number of offspring
% and thus has a higher chance of surviving in the subsequent generation.
%
% CROSSOVER is an operation which swaps genes between two chromosome.
%
% MUTATION is flipping a digit in the chromosome after crossover.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% How to use the program for your own purposes:
%
% The main thing that needs to be changed in ga.m is the fitness function
% and a few parameters. It is now set to minimize a function z=f(x,y)
% that is a sum of scaled translated Gaussian distributions with x and y
% between 0 and 30 (specified by the elements of HIGHTRAIT and LOWTRAIT).
% The function to be optimized must be changed in the section specified
% in the program below.
%
% In addition, the following parameters should be modified according
% to your needs:
%
% 1. NUM_TRAITS - the number of traits represented by each chromosome
% (e.g. if we were trying to optimize the above function
% we could have NUM_TRAITS=2).
% Example: Chromosome: 123456789012
% Possible two traits: 123456, 789012
%
% 2. HIGHTRAIT,LOWTRAIT - arrays with NUM_TRAITS elements with each element
% specifying the upper (lower) limit on each trait (often
% known a priori for an optimization problem).
% For example, we may only want to find the optimum values
% of a function over a certain region. For instance,
% we could choose HIGHTRAIT=[30 30] and LOWTRAIT=[0 0]
% for the above optimization problem.
%
% 3. SIG_FIGS - an array with NUM_TRAITS elements with each element
% specifying the number of significant figures in each
% trait (there must be at least one significant figure).
% Example: SIG_FIGS=[6 6] with above possible traits.
%
% 4. DECIMAL - an array with NUM_TRAITS elements with each element
% specifying the number of digits to the left of the decimal
% point for each trait. Choose DECIMAL to be the
% largest number integer part in HIGHTRAIT or LOWTRAIT
% to avoid saturation.
% (Example: HIGHTRAIT=[12.9 9.9], LOWTRAIT=[5.5 -5.5],
% DECIMAL=[2 1] so for our optimization problem
% let DECIMAL=[2 2])
%
% With this we can explain how to "decode" a chromosome:
%
% The first gene of each trait determines the sign of the
% trait. For our base-10 algorithm if this gene is 0-4, the trait is
% negative and if this gene is 5-9, the trait is positive.
% The remaining genes determine the size of the trait. The
% second gene is the most significant digit, while the last gene
% is the least significant digit.
%
% To determine the relative magnitude of a trait, this software uses
% the DECIMAL[] constant. This number determines how many digits to the
% left of the decimal to place the first digit of the trait (=second gene
% on the trait).
% Ex: Suppose trait #i = 123456
%
% If DECIMAL[i] = Then Trait[i]=
%
% -2 -0.0023456
% -1 -0.023456
% 0 -0.23456
% 1 -2.3456
% 2 -23.456
% etc...
%
% Note also that you must set LOWTRAIT, HIGHTRAIT, DECIMAL, and SIG_FIGS
% so that if we saturate a trait at the value of LOWTRAIT or HIGHTRAIT
% then the genes on the trait must be able to represent LOWTRAIT or
% HIGHTRAIT. For example you would not choose HIGHTRAIT(1)=100
% and DECIMAL(1)=-2 with SIG_FIGS(1)=1. Why? Hence, you must choose
% the elements of LOWTRAIT and HIGHTRAIT so that they have the same
% decimal position and number of significant figures as their
% corresponding traits.
%
%
% 5. MUTAT_PROB - the probability that a mutation will occur in any given
% gene (i.e., that a gene will be switched to another
% element of the alphabet that is being used). Note that
% for base 10 we randomly choose a number between 0-9,
% with the exception of the number that was being mutated.
%
% 6. CROSS_PROB - the probability that a crossover will occur between two
% parents (standard swapping of genes is used).
%
% 7. SELF_ENTERED - if this flag is set to 1, an initial
% population can be entered by the user. If it is 0,
% a random initial population is created. We enter the
% initial population in a special convenient way for
% the user.
%
% 8. POP_SIZE - the number of individuals in a population (fixed).
% This should be an even number - since pairs mate.
%
% 9. ELITISM - if this flag is set to one, the most fit individual in each
% generation will be carried over to the next generation
% without being modified in any way by the genetic operators.
%
% 10. DELTA, EPSILON - these are for the termination criterion.
% The program is terminated if the fitness
% changes less than EPSILON over DELTA
% generations.
%
% 11. MAX_GENERATION - the number of generations that the program will
% produce before stopping. This indicates the length
% of the program will run. If it takes too long
% to run, this number can be lowered.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Next, we provide some ideas on how to run to the program to illustrate
% its operation. Proceed as follows:
%
% (1) Experiment with the intial population. Enter your own initial
% population (to do this you modify the code and let SELF_ENTERED = 1).
% Currently, the code sets the intial population to be all zeros if
% you let SELF_ENTERED = 1. Another choice you may want to consider
% is to place the initial members of the population on some type of
% uniform grid to try to make sure that at least one of the intial
% points is in the proximity of the final point. Currently, if you let
% SELF_ENTERED = 0, the program automatically generates a random initial
% population with a uniform distribution on the allowable range of each
% trait.
% (2) Next, you can experiment with the settings for the crossover and
% mutation probabilities. First, try MUTAT_PROB=0.05; CROSS_PROB=0.8.
% If you lower the crossover probability you will generally find that
% there will be less local search near the best individuals (i.e.,
% it will tend to less often interpolate between good individuals to
% try to find a solution). If you raise the mutation probability it
% will tend to cause more random excursions into possibly unexplored
% regions. Too high of a mutation probability can cause degradation
% to random search, where too low of a mutation probability can cause
% the GA to lock in on a local minima and never find the global one (at
% least in a reasonable amount of time).
% (3) Next, you can experiment with the population size. Larger population
% sizes require more computation time but they can speed convergence.
% For instance, for a larger population size you will have more guesses
% of where the global maximum is for the intial population and hence
% it is more likely that there will be at least one good intial guess.
% There are similar effects as the algorithm runs.
% (4) Next, turn on elitism. You will find that elitism removes some of the
% effects of having a higher mutation rate since it makes sure that
% the best candidate is not disturbed.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear % Initialize momery
rand('state',0) % Reset the random number generator so each time you re-run
% the program with no changes you will get the same results.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% User: Input your own parameters below (the ones given are for the
% optimization problem given above):
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
NUM_TRAITS=2; % Number of traits in each individual
HIGHTRAIT=[30 30]; % Upper limit of a trait
LOWTRAIT=[0 0]; % Lower limit of a trait
SIG_FIGS=[6 6]'; % Number of genes in each trait
DECIMAL=[2 2]; % Order of magnitude the trait
MUTAT_PROB=0.05; % Probability of mutation (typically <.1)
CROSS_PROB=0.8; % Probability of crossover (typically near 1)
SELF_ENTERED=0; % "0": a random initial population.
% "1": a specified initial population
% If you choose "1", enter it in the program.
POP_SIZE=20; % Number of individuals in the population
ELITISM=0; % Elitism ON/OFF, 1/0
DELTA=100; % Number of generations to be counted
% for the termination criteria.
EPSILON = 0.01; % Range that the fitness must change
% in the termination criteria.
MAX_GENERATION=1000; % Number of times the loop will run before
% giving up on the EPSILON-DELTA termination cond.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Next, we want to make the initial population of size = POP_SIZE called pop.
% Each column represents a chromosome and each element in that column
% representing a gene. There will be POP_SIZE chromosomes with each of
% them having CHROM_LENGTH genes. Now, rather than directly create this
% matrix we allow the user to enter the more easily specified values of the
% traits (pop is coded to be operated on by a genetic algorithm while trait is
% a matrix with each of actual values of the parameters). Taking this approach
% you can easily enter an intial population of parameter guesses on the domain
% of the function to be maximized and the first iteration in the GA will
% convert the values to a form (stored in pop) that the genetic operators can
% work on.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
popcount=1; % Initialize the generation count,
% set it to one, the first population
if SELF_ENTERED == 0 % Make a random initial population for
% base-10 operation by specifying the
% matrix of initial traits
for pop_member = 1:POP_SIZE
for current_trait = 1:NUM_TRAITS,
trait(current_trait,pop_member,popcount)=...
(rand-(1/2))*(HIGHTRAIT(current_trait)-LOWTRAIT(current_trait))+...
(1/2)*(HIGHTRAIT(current_trait)+LOWTRAIT(current_trait));
% This starts the population off with numbers chosen randomly
% on the allowed range of variation, for this example.
end
end
else
for pop_member = 1:POP_SIZE
for current_trait = 1:NUM_TRAITS,
trait(current_trait,pop_member,popcount)=0;
% To start with a guess where all the population members are zero
end
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Next, we need to set some values based on the values input by
% the user. In particular, we determine the length of the
% chromosome and start point of each trait. This information
% will be used later in the main algorithm.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
CHROM_LENGTH=sum(SIG_FIGS)+NUM_TRAITS; % Length of the chromosome is
% the number of sig. figs. plus
% the number of sign positions
TRAIT_START(1)=1; % Initialize: the first trait
% starts at the first digit
% (this is the sign digit)
for current_trait=1:NUM_TRAITS, % Determine the start point of the
% other traits - it is the start of
% the last trait plus the no. of sig.
% figs. plus one for sign
TRAIT_START(current_trait+1)=...
TRAIT_START(current_trait)+SIG_FIGS(current_trait)+1;
% Yes, we compute the TRAIT_START for one extra trait - this
% is used for convenience in the code below.
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% In this next loop the fitness is calculated, the children are
% created and it repeats until the EPSILON-DELTA termination condition
% is satisfied or MAX_GENERATION is reached. This is the main
% loop of the GA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
while popcount <= MAX_GENERATION
% First, fix bad traits (i.e., ones that are out of the range
% specified by HIGHTRAIT and LOWTRAIT) by saturation at the extremes
for pop_member = 1:POP_SIZE
for current_trait = 1:NUM_TRAITS,
if trait(current_trait,pop_member,popcount)>HIGHTRAIT(current_trait)
% The trait has went higher than the upper
% bound so let the trait equal to the
% HIGHTRAIT bound.
trait(current_trait,pop_member,popcount)=HIGHTRAIT(current_trait);
% Now consider the other case:
elseif trait(current_trait,pop_member,popcount) rand % If true then crossover occurs
site = rand*CHROM_LENGTH; % Choose site for crossover
site = site-rem(site,1)+1; % and make it a valid integer
% number for a site
% The next two lines form the child by the swapping of genetic
% material between the parents
child(1:site,parent_number1)=parent_chrom(1:site,parent_number1);
child(site+1:CHROM_LENGTH,parent_number1)=...
parent_chrom(site+1:CHROM_LENGTH,parent_number2);
else % No crossover occurs
% Copy non-crossovered chromosomes into next generation
% In this case we simply take one parent and make them
% the child.
child(:,parent_number1)=parent_chrom(:,parent_number1);
end
end % End the "if ELITISM..." statement
end % End "for parent_number1=..." loop
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Mutate children.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Here, we mutate to a different allele
% with a probability MUTAT_PROB
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for pop_member= 1:POP_SIZE,
if ELITISM ==1 & pop_member==bestmember % If elitism on, and
% have the elite member
child(:,pop_member)=child(:,pop_member); % Do not mutate
% the elite member
else
for site = 1:CHROM_LENGTH,
if MUTAT_PROB > rand % If true then mutate
rand_gene=rand*10; % Creat a random gene
% If it is the same as the one already there then
% generate another random allele in the alphabet
while child(site,pop_member) == rand_gene-rem(rand_gene,1),
rand_gene=rand*10;
end;
% If it is not the same one, then mutate
child(site,pop_member)=rand_gene-rem(rand_gene,1);
% If takes a value of 10 (which it cannot
% mutate to) then try again (this is a very low probability
% event (most random number generators generate numbers
% on the *closed* interval [0,1] and this is why this line
% is included).
if rand_gene == 10
site=site-1;
end
end % End "if MUTAT_PROB > rand ...
end % End for site... loop
end % End "if ELITISM..."
end % End for pop_member loop
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Create the next generation (this completes the main part of the GA)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
pop=child; % Create next generation (children
% become parents)
popcount=popcount+1; % Increment to the next generation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Next, we have to convert the population (pop) to the base-10
% representation (called trait) so that we can check if the traits
% all still lie in the proper ranges specified by HIGHTRAIT and LOWTRAIT
% at the beginning of the next time around the loop.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for pop_member = 1:POP_SIZE
for current_trait = 1:NUM_TRAITS,
trait(current_trait,pop_member,popcount)=0; % Initialize variables
place_pointer=1;
% Change each of the coded traits on the chromosomes into base-10 traits:
% For each gene on the current_trait past the sign digit but before the
% next trait find its real number amount and hence after finishing
% the next loop trait(current_trait,pop_member,popcount) will be the base-10
% number representing the trait
for gene=TRAIT_START(current_trait)+1:TRAIT_START(current_trait+1)-1,
place=DECIMAL(current_trait)-place_pointer;
trait(current_trait,pop_member,popcount)=...
trait(current_trait,pop_member,popcount)+...
(pop(gene,pop_member))*10^place;
place_pointer=place_pointer+1;
end
% Determine sign of the traits and fix
% trait(current_trait,pop_member,popcount) so that it has the right sign:
if pop(TRAIT_START(current_trait),pop_member) < 5
trait(current_trait,pop_member,popcount)=...
-trait(current_trait,pop_member,popcount);
end
end % Ends "for current_trait=..." loop
end % Ends "for pop_member=..." loop
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Terminate the program when the best fitness has not changed
% more than EPSILON over the last DELTA generations. It would also
% make sense to use avefitness rather than bestfitness in this test.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if popcount > DELTA+1 & ...
max(abs(bestfitness(popcount-DELTA:popcount-1)-...
bestfitness(popcount-DELTA-1:popcount-2)))<=EPSILON
break;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end % End "for pop_count=..." loop - the main loop.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Next, provide some plots of the results of the simulation.
% These are for the example function that we seek to maximize; however
% for other applications it is not difficult to modify this code
% to show plots that are informative for the application at hand.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
t=1:popcount-1; % For use in plotting
x=0:31/100:30; % For our function the range of values we are considering
y=x;
% Compute the function that we are trying to find the maximum of.
for jj=1:length(x)
for ii=1:length(y)
z(ii,jj)=...
+5*exp(-0.1*((x(jj)-15)^2+(y(ii)-20)^2))...
-2*exp(-0.08*((x(jj)-20)^2+(y(ii)-15)^2))...
+3*exp(-0.08*((x(jj)-25)^2+(y(ii)-10)^2))...
+2*exp(-0.1*((x(jj)-10)^2+(y(ii)-10)^2))...
-2*exp(-0.5*((x(jj)-5)^2+(y(ii)-10)^2))...
-4*exp(-0.1*((x(jj)-15)^2+(y(ii)-5)^2))...
-2*exp(-0.5*((x(jj)-8)^2+(y(ii)-25)^2))...
-2*exp(-0.5*((x(jj)-21)^2+(y(ii)-25)^2))...
+2*exp(-0.5*((x(jj)-25)^2+(y(ii)-16)^2))...
+2*exp(-0.5*((x(jj)-5)^2+(y(ii)-14)^2));
end
end
% First, show the actual function to be maximized and its contour map
figure(1)
clf
surf(x,y,z);
colormap(jet)
% Use next line for generating plots to put in black and white documents.
colormap(white);
xlabel('x');
ylabel('y');
zlabel('z');
title('Fitness function');
figure(2)
clf
contour(x,y,z,25)
colormap(jet)
% Use next line for generating plots to put in black and white documents.
colormap(gray);
xlabel('x');
ylabel('y');
zlabel('z');
title('Fitness function (contour map)');
hold on
% Next, add on the evolution of the other population
% members over time.
for pop_member = 1:POP_SIZE
xx=squeeze(trait(1,pop_member,:));
yy=squeeze(trait(2,pop_member,:));
traitplot=plot(xx,yy,'k.');
set(traitplot,'MarkerSize',6);
end
% Then, add on the evolution of the best individual over
% time that the GA finds.
btraitplot=plot(bestindividual(1,:),bestindividual(2,:),'ko');
set(btraitplot,'MarkerSize',4)
hold off
% Next, plot some data on the operation of the GA
figure (3)
clf
subplot(211)
plot(t,bestfitness,'k--',t,avefitness,'k-',t,worstfitness,'k-.')
xlabel('Generation')
ylabel('Best, average, and worst fitness')
title('Best and worst fitness vs. generation')
subplot(212)
plot(t,bestindividual(1,:),'k-',t,bestindividual(2,:),'k--')
xlabel('Generation')
ylabel('Best individuals')
title('Best individuals vs. generation')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% End of program
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%