How Genetic Algorithms Work, with a MATLAB Implementation

A Concrete Example: Finding the Maximum of a Two-Variable Function

I’ve been learning about genetic algorithms lately, so I picked up a classic textbook problem to walk through the whole pipeline.

Problem: find the maximum of the following two-variable function

Individual Encoding

A genetic algorithm operates on strings, so x1 and x2 must be encoded first. Here we use binary integers: both x1 and x2 are integers in the range 1–7, each represented by a 3-bit unsigned binary integer. We concatenate x1 and x2 to encode them together as a single individual; for example, the genotype 101110 corresponds to the phenotype x=[5, 6]. The phenotype and the genotype can be converted back and forth easily.

Generating the Initial Population

The essence of a genetic algorithm is to evolve a population, so we first need a batch of initial individuals. In this example the initial population size is set to 4, and each individual’s genotype is generated randomly, for example:

011101, 101011, 011100, 111001

Fitness Calculation

A genetic algorithm uses fitness as the metric for judging how good or bad an individual is, which in turn determines its chances of reproduction (whether it is eliminated or carried forward). The optimization goal in this example is to find a non-negative maximum, so we use the function value directly as the fitness—the larger the function value, the better the individual and the higher its fitness.

Selection Operation

The selection operation passes the fitter individuals into the next generation according to some rule. The higher an individual’s fitness, the more chances it gets to reproduce. In this example we determine which individuals advance to the next generation using a probability proportional to fitness. The specific steps are as follows:

Step one, compute the sum of the fitness of all individuals:

Step two, compute the relative fitness of each individual:

This is the probability of each individual being carried over to the next generation.

Step three, the probability values form a set of intervals, with all the probabilities summing to 1.

Step four, generate a random number between 0 and 1 for each individual (4 in total in this example), and determine which individual is selected based on which probability interval the random number falls into.

Crossover Operation

The crossover operation is the main way a genetic algorithm produces new individuals: with some probability, two individuals exchange parts of their chromosomes. This example uses single-point crossover:

1. First pair up the population randomly (No. 1 with No. 2, No. 3 with No. 4)
2. Randomly set the crossover point
3. Exchange the portions of genes between the paired chromosomes

As you can see, the newly produced individuals 111101 and 111011 both have higher fitness than the two original individuals.

Mutation Operation

The mutation operation changes the gene values at one or several loci of an individual with a small probability; it is another way of producing new individuals. This example uses basic bit mutation:

1. First determine the mutation locus of each individual (the numbers in the table below indicate the locus at which the mutation point is placed)
2. Flip the bit at the mutation point with some probability

The New Generation P(t+1)

After running one round of selection, crossover, and mutation on the population P(t), we obtain the new generation P(t+1):

MATLAB Implementation

The pipeline above is implemented in MATLAB using the University of Sheffield’s Genetic Algorithm Toolbox (Sheffield GA Toolbox). First, define a custom objective function:

function z=test_shang(x,y)
z=x.^2+y.^2;
end

The main body of the genetic algorithm:

NIND=40;        % number of individuals
NVAR=2;         % number of variables (dimension) is 2
MAXGEN=100;     % maximum number of generations
PRECI=20;       % number of binary bits per variable
GGAP=0.9;       % generation gap of 0.9, discarding the 10% of individuals with the lowest fitness

trace=zeros(2,MAXGEN); % algorithm performance tracking; create a 2*MAXGEN matrix of zeros

Chrom=crtbp(NIND,PRECI*NVAR); % initialize the population

%FieldD=[10;-1;2;1;0;1;1];
FieldD=[rep(length(Chrom)/2,[1,2]); rep([-1;2],[1,2]); rep([1;0;1;1],[1,2])];

% Build the region descriptor; from left to right it is
% FieldD=[len lb ub code scale lbin ubin]
% len: the length of each substring contained in Chrom; for 1 variable it equals the number of binary bits per variable
%     rep(length(Chrom); for 2 variables it is rep(length(Chrom)/2
% lb,ub: the lower and upper bounds of each variable, here the interval [-1,2]
% code: encoding mode; code=1 is standard binary encoding
% scale: usually arithmetic scaling, scale=0
% lbin,ubin: whether to include the bounds; 0 excludes the bounds, 1 includes them

gen=0; % start at generation 0

variable=bs2rv(Chrom,FieldD);  % convert binary strings to the phenotype
ObjV=test_shang(variable(:,1),variable(:,2)); % initialize; ObjV is the optimization objective

while gen<MAXGEN
    FitnV=ranking(ObjV);                          % assign fitness values
    SelCh=select('sus',Chrom,FitnV,GGAP);         % selection operation
    SelCh=recombin('xovsp',SelCh,0.7);            % recombination
    SelCh=mut(SelCh);                             % mutation
    variable=bs2rv(SelCh,FieldD);                 % convert binary to the phenotype
    ObjVSel=test_shang(variable(:,1),variable(:,2)); % compute the objective values of the offspring
    [Chrom, ObjV]=reins(Chrom,SelCh,1,1,ObjV,ObjVSel); % reinsertion
    variable=bs2rv(Chrom,FieldD); % convert the initial population to the phenotype
    % Chrom is the initial population; SelCh holds the intermediate values, in binary
    gen=gen+1;
    [Y,I]=max(ObjV); hold on; % find the maximum and minimum within the population
    trace(1,gen)=max(ObjV);
    trace(2,gen)=sum(ObjV)/length(ObjV);
end

variable=bs2rv(Chrom,FieldD);
hold on grid;

figure(1);
plot(variable(:,1),variable(:,2),'o');

figure(2);
plot(trace(1,:));
hold on;
plot(trace(2,:),'-.');grid
legend('change in the solution','change in the population mean')

The trace matrix records the best solution and the population mean for each generation, so we can finally plot the convergence curve.

Sheffield GA Toolbox: Quick Function Reference

The FieldD Structure

FieldD=[len  lb  ub  code  scale  lbin  ubin]

• len: the length of each substring contained in Chrom (PRECI)
• lb, ub: the lower and upper bounds of each variable, respectively
• code: the encoding mode; 1 for standard binary encoding, 2 for Gray code
• scale: the scaling type; 0 for arithmetic scaling, 1 for logarithmic scaling
• lbin, ubin: whether to include the bounds; 0 excludes the bounds, 1 includes them

The multi-dimensional case:

FieldD=[rep(len(Chrom),[1,dimension]); rep([lb;ub],[1,dimension]); rep([code;scale;lbin;ubin],[1,dimension])];

Population Initialization: crtbp, crtrp, crtbase

crtbp: create a binary population

[Chrom,Lind,BaseV]=crtbp(Nind, Lind*dimension);  % for multi-dimensional functions, multiply by the dimension
[Chrom,Lind,BaseV]=crtbp(Nind, BaseV);
BaseV = crtbase([Nind,Lind],[a,b]);

Nind is the number of individuals in the population, and Lind specifies the length of an individual or chromosome.

crtrp: create a real-valued initial population

Chrom=crtrp(Nind, FieldDR);

FieldDR is a 2×Nvar matrix whose first row is the lower bound of the values and whose second row is the upper bound; the length of the Nvar vector is the length of the chromosome.

bs2rv: convert a binary string to a real value

Phen=bs2rv(Chrom, FieldD);

len denotes the chromosome length; lb and ub are the lower and upper bounds of the variable values, respectively; code(i)=1 is standard binary encoding and code(i)=0 is Gray code; scale(i)=0 is arithmetic scaling and scale(i)=1 is logarithmic scaling; lbin and ubin indicate whether the values include the bounds—0 excludes the bounds and 1 includes them.

Fitness Calculation: ranking, scaling

ranking: rank-based fitness assignment

FitV=ranking(ObjV, RFun, SUBPOP);

• RFun(1): the linear ranking scalar in [1, 2], default 2; for nonlinear ranking it is in [1, length(ObjV)-2] and specifies the selective pressure
• RFun(2): the ranking method; 0 for linear ranking (default), 1 for nonlinear ranking
• SUBPOP: the number of subpopulations in ObjV, default 1

scaling: linear fitness calculation

FitnV = scaling(ObjV, Smul);

Note: linear scaling is not suitable when the objective function returns negative fitness values.

High-Level Selection Function: select

SelCh=select(SEL_F, Chrom, FitnV);
SelCh=select(SEL_F, Chrom, FitnV, GGAP);
SelCh=select(SEL_F, Chrom, FitnV, GGAP, SUBPOP);

SEL_F is a string naming a low-level selection function, such as rws or sus; GGAP specifies the generation gap, with a default of 1, and it can also be greater than 1, allowing the offspring to outnumber the parents.

rws: roulette wheel selection

NewChrIx=rws(FitnV, Nsel);

Use roulette wheel selection to choose Nsel individuals from the population; NewChrIx holds the indices of the selected individuals.

sus: stochastic universal sampling

NewChrIx=sus(FitnV, Nsel);

reins: reinsert the offspring subpopulation into the population (required when there is a generation gap)

Chrom=reins(Chrom, SelCh);
Chrom=reins(Chrom, SelCh, SUBPOP);
Chrom=reins(Chrom, SelCh, SUBPOP, InsOpt, ObjVch);
[Chrom, ObjVch]=reins(Chrom, SelCh, SUBPOP, InsOpt, ObjVch, ObjVSel);

• InsOpt(1): the method for offspring to replace parents; 0 for uniform selection, 1 for fitness-based selection, default 0
• InsOpt(2): the proportion of offspring reinserted into each subpopulation, default 1
• ObjVch: the objective values of the individuals in Chrom, required for fitness-based reinsertion
• ObjVSel: the objective values of the individuals in SelCh, required when the number of offspring exceeds the number reinserted into the population

Crossover Functions

Recombination rule: each odd-numbered row is paired with the next even-numbered row; if the matrix OldChrom has an odd number of rows, the last odd row is not crossed over and is appended as the last row of NewChrom.

recombin: high-level individual recombination function

NewChrom=recombin(REC_F, Chrom);
NewChrom=recombin(REC_F, Chrom, RecOpt);
NewChrom=recombin(REC_F, Chrom, RecOpt, SUBPOP);

REC_F is a string naming a low-level recombination function, for example recdis, recint, reclin, xovdp, xovdprs, xovmp, xovsh, xovshrs, xovsp, xovsprs; RecOpt is the crossover probability.

Common low-level crossover functions:

Mutation Functions: mutate, mut, mutbga

mutate: high-level mutation function

NewChrom=mutate(MUT_F, OldChrom, FieldDR, MutOpt, SUBPOP);

MUT_F is a string naming a low-level mutation function (such as mut or mutbga); SUBPOP defaults to 1 and determines the number of subpopulations in OldChrom.

mut: discrete mutation operator

NewChrom=mut(OldChrom, Pm);
NewChrom=mut(OldChrom, Pm, BaseV);

Pm is the mutation probability, default 0.7/Lind.

mutbga: mutation for real-valued populations

NewChrom=mutbga(OldChrom, FieldDR, MutOpt);

• MutOpt(1): the mutation probability, default 1/Nvar
• MutOpt(2): a scalar in [0, 1] that shrinks the recombination range, default 1 (no shrinking)

Other Functions

rep: matrix replication function

MatOut=rep(MatIn, REPN);

REPN is a two-element vector; REPN(1) specifies the number of vertical replications and REPN(2) specifies the number of horizontal replications. rep is a low-level replication function that is usually not called directly; it is invoked by many functions in the GA Toolbox.

migrate: migrate individuals between subpopulations

Chrom = migrate(Chrom, SUBPOP);
Chrom = migrate(Chrom, SUBPOP, MigOpt);
Chrom = migrate(Chrom, SUBPOP, MigOpt, ObjV);
[Chrom,ObjV] = migrate(Chrom, SUBPOP, MigOpt, ObjV);

MigOpt is an optional vector with three parameters:

• MigOpt(1): the migration probability, default 0.2
• MigOpt(2): the migration selection method; 0 for uniform migration (default), 1 for fitness-based migration
• MigOpt(3): the migration topology; 0 for a fully connected (mesh) structure (default), 1 for a neighborhood structure, 2 for a ring structure

ObjV is a column vector of the objective values corresponding to all individuals, required when MigOpt(2)=1.