In Chapter 7, we discussed many optimization algorithms and inbuilt MATLAB functions, and we saw examples of their use. However, almost all of these algorithms are gradient-based and return only local minima. If the objective function has many local minima and we start from an initial guess close to one such minimum, the algorithm will most probably output the local minimum closest to the initial guess. In such cases, we would never reach the global minimum. Hence the success of the method depends highly on the initial guess. Most real world objectives are multivariate functions and contain multiple minima. Therefore, new heuristic search-based algorithms have been devised. These algorithms start with multiple initial guesses and result in solutions which evolve with time. Hence these methods are known as evolutionary computation methods.

Every initial guess initiates its own optimization chain containing values updated at each step. At any step, the value from a chain is known as a candidate solution. The set of all candidate solutions from a single step is known as a population. Since the updates are performed at every step for all the chains, we need to use simple and less extensive computations at each step rather than performing steps involving gradients or Hessians. In this chapter, we will study two algorithms of this kind and learn how we can implement them. These algorithms are inspired by nature and natural processes, motivated by the belief that nature is the best optimizer.

The Rastrigin Function

In this chapter, we will consider a special benchmark function, known as the Rastrigin function, which is used to evaluate the performance of optimization methods. It is a non-convex function with multiple local minima and one global minimum. This is a perfect example which illustrates why a gradient-based optimization algorithm will not always work. The Rastrigin function in two dimensions is given by

which is plotted in Figure 8-1 for x∈[-1.2 ,1.2]. As is clear from the figure, this function has multiple minima and one global minimum. In the next two sections, we will discuss two evolutionary computation algorithms, namely PSO and GA, which are able to solve these kinds of difficult optimization problems. Before moving forward, let us first define this function in MATLAB in the following way

function f=rastrigin(x)
s=x.^2-10*cos(2*pi*x);
f=20+sum(s);

The above implementation of the function accepts a vector x representing a candidate solution and computes the value of the Rastrigin function, treating each element of the vector as one dimension/variable. Let us define a vectorized version of the Rastrigin function which accepts a matrix Xpop representing the whole population. Here each row of Xpop is a candidate solution. The following MATLAB function implements the same.

function f=rastrigin_vec(Xpop)
s=Xpop.^2-10*cos(2*pi*Xpop);
f=20+sum(s,2);

Here the output f is a vector with each element representing the function value at the corresponding candidate solution.

Particle Swarm Optimization

Particle swarm optimization (PSO) is an optimization algorithm from the class of evolutionary computation methods and is inspired by the procedure by which a swarm of birds finds its food or optimum habitat.

In PSO, the population is known as a swarm and each candidate solution is called a particle. The objective is to minimize a function f(x) where x denotes a particle. We start with an initial population. For any particle, a position x is said to be better than a position y if f(x)<f(y) Now the population is updated by the movement of particles around the search space. Each particle i remembers the best position P_i attained by it in the past. Also the swarm remembers the global best value G attained by any particle in the past. Each particle's movement is determined by its current velocity, its own best known position and the global best value of the swarm. The swarm is expected to move towards the global best solution or at least towards a satisfactory solution for which the function value is close to the global minimum.

Algorithm

The algorithm is given as follows:)

Initialization

1. Define N=size of population.
2. Generate N particles randomly in the search space.
3. Denote the location of the i^th particle by x_i.
4. Randomly generate the velocity v_i of the particles.
5. Define the velocity weight factor w = 0.5.
6. For each particle, set P_i=x_i.
7. Compute the value f_i as f_i=f(x_i) for each particle.
8. Find the particle that has the minimum  f_i among all particles. Set G to it.

   Iterations and Updates

9. Now for each particle i
   1. Generate two uniform random variables r₁ and r₂ .
   2. Compute v_(i,next)=wv_i+C₁r₁(P_i-x_i)+C₂r₂(G-x_i).
   3. Compute x_(i,next)  =x_i+v_i.
   4. If f(x_i,next)<f(x_i)
      1. Update P_i=x_i,next.
   5. Update the velocity vi of the particle by v_i,next.
   6. Update the location of the particle by x_i,next.
10. Compute the minimum of all f(P_i) and set the corresponding particle to G.
11. If the terminate conditions are met, go to step 11, otherwise go to step 8.

    Output

12. Output G as the solution of the optimization.

Here C₁ and C₂  are two parameters which define the weighting assigned to the guidance given by the local best value and the global best value.

Implementation

In this section, we will see how we can implement PSO as a generic function. Our goal is to create a function which accepts a handle to the objective function and a few configuration parameters as inputs, performs the steps of the PSO algorithm and outputs the optimal solution. Let us first implement it as a script and assume that the following three variables are given to us)

objfunction %handle to function
nvars %number of independent variables
opt %structure with following configuration fields
%opt.C1    % Cognitive parameter
%opt.C2    % Social parameter
%opt.N     % Size of swarm
%opt.UB    % Upper limit of search space, must be a 1xnvars vector
%opt.LB    % Lower limit of search space, must be a 1xnvars vector
%opt.MaxIter %maximum number of iterations
%opt.MinIter % minimum number of iteration to terminate when
                               % global best is not changing much
%opt.TolG  % tolerance to compare the change in global best

Now we will implement each part of the algorithm.

Initialization

N=opt.N;
X= repmat(opt.LB,N,1)+repmat(opt.UB-Out.LB,N,1).*rand(N,nvars);
Vrange=2*opt.UB;
V=-repmat(Vrange,N,1)+2*repmat(Vrange,N,1).*rand(N,nvars);
f=objfunction(X);
P=X; %initialize pbest with the initial particles locations
[Gbestval G]=min(f);
termination=0; %termination flag to decide when to terminate
iter=1;

Note that we initialized each particle’s velocity by creating a vector with uniform random variables between 0 and 1 and then scaled it so that it falls inside the range [2opt.LB, 2opt.UB].

Iteration Loop

while (termination ==0)

Updates

C1=opt.C1;C2=opt.C2;
               for i=1:N
                       %compute next step velocities and locations
                       vnext=0.5*V(i,:)+...
                               C1*rand(1,2).*(P(i,:)-X(i,:))...
                               +C2*rand(1,2).* (G-X(i,:));
                      xnext=X(i,:)+vnext;
                       %update the pbest if the next computed values are
                       %better than the current values
                       if objfunction(xnext)<objfunction(P(i,:)
                               P(i,:)=xnext;
                       end
                       V(i,:)=vnext;
                       X(i,:)=xnext;
               end
               f=objfunction(P);
               %compute the global best particle
               [Gbestval indexG]=min(f);
               G=P(indexG,:);
               Sol(iter,:)=G;
               iter=iter+1;

Termination Conditions

if iter>opt.MaxIter
        termination=1;
end
if iter>opt.MinIter+1
        %check if the value of G hasn't changed from past
        %MinIter number of iterations.
        e=abs(Sol(iter-1,:)-Sol(iter-opt.MinIter-1,:));
        if e<opt.TolG
                termination=1;
        end
end

Iteration Ends and Outputs

end
optSol=Sol(end,:);
optValue=objfunction(optSol);
optTrace =Sol;

Function Definition

Now we will transform the script into a primary function by adding the following line at the start of the script file.

function [optSol optValue optTrace]=solvePSO(objfunction,nvars,opt);

Example

Now we can define the objective function that we want to optimize along with configuration parameters in a separate script/function file and call the solvePSO in the following way:

objfunction =@rastrigin_vec;
nvars =2;
opt.C1    = 2;
opt.C2    = 2;
opt.N     = 30;
opt.UB    = [1.2 1.2];
opt.LB    = [-1.2 -1.2];
opt.MaxIter =60;
opt.MinIter =30;
opt.TolG  = 1e-5;
[optS optV trace]=solvePSO(objfunction,nvars,opt);

Then the following code can plot the objective function value at the global best particle at each step (see Figure 8-2)

plot(1:size(trace,1),objfunction(trace),'LineWidth',2);

We can also plot the movement of the global best particle with respect to iteration time (see Figure 8-3) in the following way:

plot(trace(:,1),trace(:,2),'LineWidth',2);

The Genetic Algorithm

The genetic algorithm (GA) is a particular example of an evolutionary algorithm (EA) that uses techniques inspired by evolutionary biology such as inheritance, mutation, selection, and crossover (also called recombination).  It imitates the natural reproduction and selection process where genes combine to generate a better offspring. We will discuss different steps of the genetic algorithm and the way to implement each step. In GA, the objective function is also called the fitness function as it determines the capability of any candidate solution to survive and reproduce.

Representation

In a natural evolution process, there are multiple chromosomes (each representing one individual) at any generation and they make up the population. In GA, we also have multiple candidate solutions at any step (a step is also known as a generation) which make up the solution population. Each candidate solution imitates a single chromosome or genotype (hence represents an individual) and is represented by a binary vector where each binary value represents a gene (see Figure 8-4). We will fix the length of the vector to simplify the key operations. For example, if there are five candidate solution values (x= 5,6,7,2,10) in any population, we can represent the population by the following matrix

Xpop=   [0 0 0 0 0 1 0 1;  ...
         0 0 0 0 0 1 1 0;  ...
                        0 0 0 0 0 1 1 1;  ...
                        0 0 0 0 0 0 1 0; ...
                        0 0 0 0 1 0 1 0];

If there is more than one variable, we can append the binary representation of all the variables to make a big vector. For example, if there are three candidate solutions (x= (1,5),(6,10),(7,2)) in any population, we can represent the population as

Xpop=   [0 0 0 0 0 0 0 1 0 0 0 0 0 1 0 1;  ...
         0 0 0 0 0 1 1 0 0 0 0 0 1 0 1 0;  ...
         0 0 0 0 0 1 1 1 0 0 0 0 0 0 1 0];

Initialization

Let N denote the population size, nvars denote the number of variables and nbits denote the number of bits needed to represent each variable. We can see that each chromosome can be represented by nvars*nbits bits. To initialize the population, we will first generate nvars*nbits random values from a uniform distribution for each chromosome and then convert each of these random values to a binary value.

N=opt.N;
nbits=opt.nbits;
XInitPop=rand(N,nvars*nbits)<0.5;
XPop    = XInitPop;

Selection

At each generation, some candidates are selected for reproduction. This selection is made according to their fitness value. In this implementation, we would select the top 50% of the population to reproduce. Let fitnessfunction be the handle of the fitness function.

f=fitnessfunction(XPop);
selectionratio =0.5;

%sort the chromosomes in terms of their fitness values
[fsort indexsort]=sort(f,1,'descend');
XPop=XPop(indexsort,:);

%select the top selectionratio % of the total population
indexselect= 1:floor(N*selectionratio);
nselect = length(indexselect);

Crossover

In crossover operations, two chromosomes attach to each other at some body location and exchange their body with each other as illustrated in Figure 8-5. The following code divides the selected pool into two partitions (say males and females), takes one chromosome from each pool, randomly generates an attaching point in them, performs the crossover to generate two new chromosomes and repeats this process until all selected chromosomes have participated in reproduction.

%shuffle the order of selected chromosomes
Indexselect     =indexselect(randperm(nselect));

%divide into two pools
indexmales      =indexselect(1:floor(nselect/2));
indexfemales    =indexselect(1+floor(nselect/2):end);
nchildren= floor(nselect/2)*2;
for i=1:length(indexmales)
       father= XPop(indexmales(i),:);
       mother= XPop(indexfemales(i),:);
       %find random attaching bit position
       attachbit=1+floor(rand(1,1)*nbits);
       %perform crossover
       child1= [father(1:attachbit-1) mother(attachbit:end)];
       child2= [mother(1:attachbit-1) father(attachbit:end)];

       children(i*2-1,:)       = child1;
       children(i*2,:) = child2;
end

Note to generate the random bit position for attachment, we first generate a random uniform variable between 0 and 1. We then scale it so that it falls between 0 and nbits. Then we take the floor of this number which will result in a random integer uniformly distributed between 0 and nbits. We finally add 1 so that it becomes a random integer uniformly distributed between 1 and nbits.

Mutation

In the mutation operation, a randomly picked gene of a randomly picked chromosome gets inverted, leading to a child which may have very different properties than the parent. The following code takes 10% of the selected population (which is 5% of the complete population), selects a random bit position at each selected chromosome and flips the bit to generate a new child.

%find random 10% of the mutants from the selected population
indexmutants=indexselect(randperm(length(indexselect)));
nmutants=floor(nselect*0.1);
indexmutants=indexmutants(1:nmutants);

%generate the nmutants number of mutation bit locations, one for %each mutant parent chromosome
mutationbits=1+floor(rand(nmutants,1)*nbits);

%select the mutant values
valuemutants = XPop(indexmutants,: );

%flip the bit
valuemutants(:,mutationbits)=1-valuemutants(:,mutationbits);
mutantchildren=valuemutants;

New Generation

Finally we let the bottom 55% (50+5) of the population die and fill it with children to create the next generation.

XPopNew=XPop;
XPopNew(end-nmutants+1:end,:)=mutantchildren;
XPopNew(end-nmutants+1-nchildren:end-nmutants,:)=children;
XPop=XPopNew;

Store the Best Chromosome of the Generation

We will store the best chromosome at each generation as a trace to observe the convergence behavior of our implementation.

f=fitnessfunction(XPop);
[fbest indexbestchromosome]=max(f);
optTrace(iter,:)=XPop(indexbestchromosome,:);

Termination Conditions

iter=iter+1;
if iter>opt.MaxIter
        termination=1;
end

Iterations

The above discussed steps involving selection, crossover, mutation and formation of a new generation must be repeated using a while loop until the termination condition is met.

iter=1;
while(termination==0)
       % code follows here

end

Output

optSol=optTrace(end,:);
[optV optSolReal]=fitnessfunction(optSol);

Note that for the last line of the above code, the fitness function must adhere to a syntax. See the example section to learn the exact syntax required to run this implementation.

Function Definition

Let us convert the above written code into a function by adding the following lines at the beginning of the file.

function [optSol optSolReal optV optTrace] =solveGA(fitnessfunction,nvars,opt)

Example

We will now apply our solveGA function to solve the Rastrigin optimization problem. Let us first define an intermediate interface function which can process the chromosome representation of the population and generate the fitness of each candidate solution from the rastrigin_vec function.

First we need to define how the binary representation is linked to the actual value. Suppose the range of the variable is [a,b] and the number of bits assigned to each variable is nbits. Note that a binary string with nbits will represent a value between 0 to 2^nbits–1 which can be scaled to a value between 0 to 1 by dividing the original value by 2^nbits–1. Finally, we can scale the result again to transform it to the corresponding value between a and b.

function [f y]=fitness_interface(XPop,nvars,nbits,range)
a=range(1);
b=range(2);
for i=1:size(XPop,1)
        x=XPop(i,:);
        for j=1:nvars
                bitrange= (j-1)*nbits+(1:nbits)
                y(i,j)=bin2dec(num2str(x(bitrange)));
        end
end
y=a+(b-a)*y/(2^nbits-1);
f=-rastrigin_vec(y);

Notice the negative sign. It is there because GA maximizes the fitnessfunction whereas our goal here is to minimize the Rastrigin function.

The following code defines the final fitness function that can be passed to solveGA along with configuration parameters and calls the solveGA:

nvars=2;
nbits=10;
range=[-1.2 1.2];
fitnessfunction=@(x) fitness_interface(XPop,nvars,nbits,range);
opt.MaxIter=70;
opt.N=50;
opt.nbits=nbits;
[optSol ptSolReal optVal optTraceoptTrace]= ...
                               solveGA(fitnessfunction,nvars,opt);

We can plot the optimum best fitness value achieved at each generation (see Figure 8-6) using the following code:

[optTraceValue optTraceReal]=fitnessfunction(optTrace);
plot(optTraceReal(:,1),optTraceReal(:,2),'-','LineWidth',2);

The Inbuilt Function ga

MATLAB provides an inbuilt method ga to solve any optimization problem using the genetic algorithm. Note that you need to have the Genetic Algorithm Toolbox installed to use this inbuilt function. Like other inbuilt optimization functions, we need to first define the objective function or the fitness function. For the example we considered, we can implement the fitness function in the following way

function f=fitness_rastrigin(x)
s=x.^2-10*cos(2*pi*x);
f=20+sum(s);

Note that the inbuilt function ga will try to minimize the fitness function. This is opposite to the previous convention used in our GA implementation, where we wanted to maximize the fitness function.

We can also implement the function in the vectorized fashion so that it accepts the complete population XPop as input where each row in XPop represents a single candidate solution. This can help ga to run faster, which will be investigated later.

Now in a separate function/scriptfile, we define nvars by the total number of parameters, which is 2 here. To use the genetic algorithm at the command line, we can now call the genetic algorithm function ga with the following syntax

nvars=2;
[xopt fval] = ga(@fitness, nvars)

The function ga will return the following outputs

xopt      Optimal point at which the final value is attained

fval      Value of the fitness function at xopt

Simulation options

The genetic algorithm has many configuration parameters and its performance can be fine-tuned by adjusting parameters. The function ga also accepts an extra input options in the following syntax.

 [x fval] = ga(@fitnessfun, nvars, options)

Here options is a structure containing various simulation parameters and configuration settings for the genetic algorithm. If we do not pass this argument, ga uses the default options. To create the options structure, we should use the inbuilt function gaoptimset. With its help, we need to pass only the options which we want to change and the rest of the options which are not passed are automatically set to their default values. The structure of the variable options is given below.

options =
        PopulationType: 'doubleVector'
          PopInitRange: [2x1 double]
        PopulationSize: 20
            EliteCount: 2
     CrossoverFraction: 0.8000
        ParetoFraction: []
    MigrationDirection: 'forward'
     MigrationInterval: 20
     MigrationFraction: 0.2000
           Generations: 100
             TimeLimit: Inf
          FitnessLimit: -Inf
         StallGenLimit: 50
        StallTimeLimit: Inf
                TolFun: 1.0000e-006
                TolCon: 1.0000e-006
     InitialPopulation: []
         InitialScores: []
        InitialPenalty: 10
         PenaltyFactor: 100
          PlotInterval: 1
           CreationFcn: @gacreationuniform
     FitnessScalingFcn: @fitscalingrank
          SelectionFcn: @selectionstochunif
          CrossoverFcn: @crossoverscattered
           MutationFcn: {[1x1 function_handle]  [1]  [1]}
    DistanceMeasureFcn: []
             HybridFcn: []
               Display: 'final'
              PlotFcns: []
            OutputFcns: []
            Vectorized: 'off'
           UseParallel: 'never'

Consider an example where we want to increase the population size to 200 and leave the other options unchanged. The following code will generate the desired options structure

Op=gaoptimset('PopulationSize',200);

and then we can call the ga function with the above options structure

X=ga(@fitness,2,Op);

We can see a configuration parameter ‘Vectorized’ which is set to off by default. This means that MATLAB expects a fitness function that is not vectorized and therefore ga will call fitnessfunction once for each candidate solution. Instead, we can change the fitness function in the following way so that the complete population can be passed at once to the fitnessfunction and it should be able to evaluate the fitness for all candidate solutions from that population.

function f=fitness_rastrigin_vec(Xpop)
s=Xpop.^2-10*cos(2*pi*Xpop);
f=20+sum(s,2);

Then we should let GA know that fitnessfunction is now vectorized by setting the ‘Vectorized’ option to on in the options structure and passing this options structure to ga

Op=gaoptimset('Vectorized','on');
X=ga(@fitnessvec,2,Op);

GUI Interface

MATLAB provides a graphical user interface to solve various kinds of optimization problems. This is helpful to those researchers who don’t want to write a program and are more comfortable working with GUI-based software. To open the GA graphical user interface, type the following in the command window

optimtool('ga');

This will open a GA GUI. Now in the GUI, we can set different options and constraints. Once we are done, we can run the genetic algorithm method by clicking the run button. We can also see real time plots of convergence as the algorithm advances.

Summary

Evolutionary computation methods are meta-heuristic methods that help us to solve complex problems. They are able to learn the optimal pattern during execution. Two important classes of evolutionary computation methods are swarm intelligence and evolutionary algorithms. In this chapter, we have learnt two such algorithms, one from each class. The first algorithm, PSO, is a swarm intelligence algorithm and mimics the way birds find their food. The second algorithm, GA, is an evolutionary algorithm which mimics the process of natural evolution to generate the best species. We implemented both algorithms from scratch. As a side product, this chapter has also taught us in great detail how we can implement a complex algorithm.