In [45]:
%matplotlib inline
In [46]:
from ga_auxiliary import *
from ga_solver import *
CIVL4250 Numerical Methods in Engineering
Dr Dorival Pedroso
The code for this notebook is available in https://github.com/cpmech/CIVL4250py (browse 2015 => week07)
A slightly more general Python code for teaching/learning Simple Genetic Algorithms is availalble in https://github.com/cpmech/tlga
Notes:
A Genetic Algorithm (GA) is a computational method to find successive improvements to a set of initial solutions in optimisation problems. In other words, it is an optimisation method.
The essence of GAs is to mimic Darwin's natural evolution described by genetics and based on the survival of the fittest. The algorithm comprises of steps such as fitness assessment, selection, crossover and mutation -- all try to make an analogy to Nature.
Although it sounds complicated, GAs are very simple. To illustrate, suppose the minimum value of $$ y(x) = -x \, \sin{(x)} $$ has to be found in $[0, 4\,\pi]$.
The function $y(x)$ is known as model function and is directly related to the objective function explained later on. Its results will serve to compute the objective values that are the key component to decide whether one $x$ value is "good" or "bad".
This function is plotted below
In [47]:
# model function
def y(x): return -x * sin(x)
# plot model function
xcurve = linspace(0.0, 4*pi, 101)
plot(xcurve, y(xcurve), label='y(x)')
Gll('x', 'y')
where we can see that both local and global (minima) values exist. The first is around $x=2$ and the second is around $x=8$. Because of the probable existence of local minima, conventional methods have difficulties when solving this kind of problem whereas a simple GA can succeed.
In [48]:
# initialise random numbers' generator
seed(1234) # use a fixed seed, so every time we run this code
# we will get the same results
# population
ninds = 10 # number of individuals: population size
xmin, xmax = 0.0, 4.0*pi # limits
X = Random(ninds, xmin, xmax) # generate numbers between 0 and 4*pi
# just for the sake of showing nice numbers in this note:
X = array(X, dtype=int) # truncate values
X = array(X, dtype=float) # revert to float
print 'Population:\n', X
For illustration purposes, these values, the initial population, are plotted below
In [49]:
plot(xcurve, y(xcurve), label='y(x)')
plot(X, y(X), 'ro', label='population')
Gll('x', 'y', 'upper left')
where we can see that they are all over the range from $0$ to $4\,\pi$ and that three values are good candidates, namely $x=7$, $x=9$ and $x=9$ again.
In the GA, each one of these values $x \in X$ is called an individual. All of them ($X$) is known as a population.
Each individual is further subdivided into bases and the resulting set defines a chromosome. Any kind of subdivision could be applied to $x$. Here, we apply a simple random subdivision.
For example, if $x=10$, its chromosome is obtained as follows:
In [50]:
nbases = 5
c10 = SimpleChromo(10, nbases)
print c10, ' sum(c10) =', sum(c10)
where the number 10 is split into 1.363, 1.908, 2.603, etc. and the resulting sum gives $x=10$ back. Thus, the representation model for mapping individuals (actually genes) into their chromosome is just a simple addition formula.
The two following formula implements the representation model for future use:
In [51]:
# compute chromosome from x
def cFcn(x): return SimpleChromo(x, nbases)
# compute x from chromosome
def xFcn(c): return sum(c)
During reproduction process of a GA, the bases of the chromosome from different parents will be mixed with each other in order to create new individuals. During this process, crossover and mutation may happen.
Now, the chromosomes of all individuals in the population $X$ are computed as follows:
In [52]:
C = array([cFcn(x) for x in X])
print C
The objective function is similar to the model function; but also implements the constraint condition $0 \leq x \leq 4\,\pi$. To do this, one technique is to apply a penalty to values outside this range.
Furthermore, it is more convenient for the definition of a general GA solver to implement the objective function as a function of the chromosome.
The adopted objective function is:
In [53]:
def objFcn(c):
x = xFcn(c)
if x < xmin or x > xmax: return 100.0*(1.0+abs(x))
return y(x)
where we use $100(1+|x|)$ as a penalty for values outside $[0, 4\,\pi]$.
All objective values can now be computed by means of
In [54]:
Y = array([objFcn(c) for c in C]) # objective values
PrintPop(C, Y, xFcn, showC=True)
Before the reproduction stage, individuals must be selected (in pairs) and, to do so, a measure of fitness must be defined.
For our minimisation problem, the fitness function $f(x)$ can be simply a linear model that maps $y(x)$ to a value in $[0, 1]$ as follows $$ f(x) = \frac{y_{max} - y(x)}{y_{max} - y_{min}} $$ Thus, if $x_{best}$ is such that $y(x_{best})=y_{min}$, then $f(x_{best})=1$. On the other hand, if $x_{worst}$ is such that $y(x_{worst})=y_{max}$, then $f(x_{worst})=0$. Therefore, $1$ means "good" (fittest) and $0$ means "bad".
Furthermore, for the next computations, the population must be sorted in decreasing order of fitness; i.e. the best ones are listed first.
The fitness of all individuals and their sorting are computed as follows:
In [55]:
# compute fitness from objective values
F = Fitness(Y)
# sort in decreasing order of fitness
C, Y, F = SortPop(C, Y, F)
PrintPop(C, Y, xFcn, F)
We can see that individuals $x=7$ and $x=9$ (there are two of those) are among the best ones (smallest $y$ values). This can be easily verified in the above figure; the corresponding points are around the global minimum.
The next step before selecting individuals is to compute the probability of each individual going to the reproduction stage. This is simply accomplished by dividing the fitness value by the sum of fitness: $$ p_i = \frac{f_i}{\sum_k^{nind} f_k} $$ where $p_i$ is the probability of individual $i$ with $x_i$.
Later, the cumulated probability will also be required and we use $M$ to hold the corresponding values.
In [56]:
# probabilities
P = F / sum(F)
# cumulated probabilities
M = cumsum(P)
PrintPop(C, Y, xFcn, F, P, M)
In a simple GA, it is common to replace the population with the same number of individuals. Therefore, in a population with $N_{ind}$ individuals, $N_{ind}/2$ pairs are selected.
The selection has to be in such a way that the genetic data of the fittest individuals have higher chance of survival; i.e. be passed to the new generation. That's when the probability (of survival) come into place.
One simple and common method to perform a biased selection of an individual is the so-called Roulette Wheel selection. In this method, one can think that the roulette is made in such a way that each bin corresponds to an individual and its size is proportional to the individual's probability. Then the roulette is spun in one direction and a ball is thrown into it along an opposite direction. One individual is then found where the ball lands. The process is repeated to select more individuals.
To visualise the selection algorithm using the roulette wheel method, the following figure is created where the bins are distributed over the horizontal line.
Suppose 5 individuals are to be selected.
In [57]:
PlotProbBins(X, P)
where we see that the bins for individuals $x=2$ and $x=7$ are the largest ones.
With the aid of above figure, suppose 5 individuals are to be selected. First, we generate 5 random numbers in $[0, 1)$
In [58]:
sample = random(5)
print sample
Now, looking at the figure and comparing the cumulated probabilities, we can see that the following individuals are selected: $$ selected = [2.0, 5.0, 12.0, 2.0, 2.0] $$
A function to do this is available in the auxiliary file; for example
In [59]:
S = RouletteSelect(M, 5, sample)
print 'indices of selected individuals:\n', S
With two individuals selected, there is a probability ($p_{c}$) that their genetic data will be mixed during reproduction. This known as the crossover. Mathematically, the crossover is carried out by randomly selecting the index ($pos$) of a base which is then used to partially copy the data into the new generation.
This can be done by means of
In [60]:
# crossover between the first two individuals
pos = randint(1, nbases-1)
A, B = C[0], C[1]
a = hstack([A[:pos], B[pos:]])
b = hstack([B[:pos], A[pos:]])
print 'pos =', pos
print 'parent A =', A
print 'parent B =', B
print 'child a =', a
print 'child b =', b
and is illustrated below
In [61]:
DrawCrossover(A, B, a, b, pos)
where we see that the descendants will have bits of genetic data from both parents.
With a small probability ($p_m$), from time to time, bases from the chromosome are arbitrarily changed. This represents the mutation phenomenon in genetics and it actually helps the discoverty of new data that could be useful (optimal). But it also lead to worse cases. This is not a problem since natural selection will take care of the worse cases.
The effect of mutation is to widen the search space and therefore is an essential component to avoid getting stuck in local minima.
Mathematically, anything could be used as a mutation strategy. Nonetheless, some analysis for each particular problem must be made.
Here, we model mutation by simply and randomly adding or subtracting one base in the chromosome by 10% of the largest value.
In [62]:
# apply mutation to new individual 'a'
print 'before: a =', a
pos = randint(0, nbases) # position
bmax = max(a)
if FlipCoin(0.5): a[pos] += bmax * 1.1
else: a[pos] -= bmax * 1.1
print 'after: a =', a
print 'pos =', pos, ' bmax =', bmax
Finally, the steps explained above are combined in a loop known as evolution, where new populations are successively created by replacing the old one.
After a number $n_{gen}$ of generations (a.k.a "years"), the individuals will be closer to the optimal solution. Nonetheless, as in Nature, there is no way to tell when "evolution has to be stopped". Therefore, it is hard to decide on the maximum number of generations.
Some strategy to analyse the speed of convergence could be adopted for the best objective value; however, here we simply choose a "large" number of generations.
The crossover and mutation steps are usually called the reproduction step.
We also note that the algorithm performs better if, at each generation, the best individual from the previous population is copied into the place of the worst individual of the new population. In this way, optimal solutions found already are never lost. This step is commonly known as elitism.
The resulting Simple Genetic Algorithm (SGA) is given below
In [63]:
def Evolve(C, xFcn, objFcn, ngen=10, elite=True, pc=0.8, pm=0.01, verb=False, showC=False):
# objective values
Y = array([objFcn(c) for c in C]) # objective values
ninds = len(C)
nbases = len(C[0])
# fitness and probabilities (sorted)
F = Fitness(Y)
C, Y, F = SortPop(C, Y, F)
P = F / sum(F)
M = cumsum(P)
# results
OV = zeros(ngen+1)
OV[0] = Y[0] # best first objective value
# evolution
for gen in range(ngen):
# best individual
bestC = C[0].copy()
bestY = Y[0]
# print generation
if gen==0 or verb:
print
PrintPop(C, Y, xFcn, F, showC=showC)
# selection
S = RouletteSelect(M, ninds)
idxA, idxB = FilterPairs(S)
# reproduction
Cnew = [] # new chromosomes
for k in range(ninds/2):
# parents
A, B = C[idxA[k]], C[idxB[k]]
# crossover
if FlipCoin(pc):
pos = randint(1, nbases-1)
a = hstack([A[:pos], B[pos:]])
b = hstack([B[:pos], A[pos:]])
else:
a, b = A.copy(), B.copy()
# mutation
if FlipCoin(pm):
pos = randint(0, nbases)
bmax = max(a)
if FlipCoin(0.5): a[pos] += bmax * 1.1
else: a[pos] -= bmax * 1.1
# new individuals
Cnew.append(a)
Cnew.append(b)
# new population
C = array(Cnew)
Y = array([objFcn(c) for c in C]) # objective values
F = Fitness(Y)
# elitism
if elite:
I = F.argsort()[::-1] # the [::-1] is a trick to reverse the sorting order
best = I[0]
worst = I[ninds-1]
if bestY < Y[best] and bestY < Y[worst]:
C[worst] = bestC
Y[worst] = bestY
F = Fitness(Y)
# probabilities (sorted)
C, Y, F = SortPop(C, Y, F)
P = F / sum(F)
M = cumsum(P)
# objective values
OV[gen+1] = Y[0] # best current objective value
# results
return C, Y, OV
and can be used as follows:
In [64]:
# input data
ninds = 10 # number of individuals: population size
nbases = 5 # number of bases in chromosome
ngen = 10 # number of generations
pc = 0.8 # probability of crossover
pm = 0.01 # probability of mutation
elite = 1 # use elitism
verb = False # verbose
# run GA
C, Y, OV = Evolve(C, xFcn, objFcn, ngen, elite, verb=verb, showC=False)
X = [xFcn(c) for c in C]
# print and plot
subplot(2, 1, 1)
plot(xcurve, y(xcurve), label='y(x)')
print '\nFinal population:'
PrintPop(C, Y, xFcn)
sol = '\nSolution: x=%g y=%g' % (X[0], Y[0])
print sol
# plot results
subplot(2, 1, 1)
plot(X, Y, 'k*', label='final population')
Gll('x', 'y', 'upper left')
# plot convergence graph
subplot(2, 1, 2)
G = range(ngen+1)
plot(G, OV, 'b.-', label=sol)
Gll('generation', 'y(x)')
show()
Pedroso DM and Williams DJ (2011) Automatic calibration of soil-water characteristic curves using genetic algorithms. Computers and Geotechnics. 38(3):330-340 10.1016/j.compgeo.2010.12.004