Week 4 - Genetic Algorithm

Aim:

To write a MATLAB code for optimizing and thereby finding global maxima of stalagmite function.

Objectives:

To plot the graph to see how maximum value of function varies with x, y and iterations for 3 cases:

• Unbounded case
• Bounded case
• Increased population size.

Theory:

Genetic Algorithm (GA) is a population-based search and optimization tool or technique. It works on the mechanism of Natural genetics and Darwin’s Principle of Natural selection.

It mainly employs population-based searching, where the number of solutions at any iterations are more than one. In GA, solutions are represented in coded form like chromosomes. Chromosomes consists of genes, similarly in GA the coded solutions are represented in string, each string consists of number of bits. In GA, the solutions which are superior have better chance to get selected for next generation (Survival of the fittest), fitness is calculated using objective function.

General outline of GA:

• Firsts of all, it generates initial population of solution randomly. One string of solution is equivalent to one chromosome, each string consists of number of bits. The length of string also depends on accuracy we want. E.g. let 4 bits are used to represent each sub-string, so size of string = 4+4 = 8………………… (for two variables).

Then those substrings are further decoded to get values of x1, x2 as in above example.

• Fitness Evaluation: To determine fitness values of each solutions real values of each solutions are taken into account as done in first step. Then fitness values of that particular solution can be calculated by putting respective variable value in objective function (here, in our case objective function is stalagmite function).

• Reproduction: It is also called as solution of mating pool. In reproduction process, the solutions which are good in terms of fitness values are selected for mating pool. Reproduction process is probabilistic in nature i.e. the fitter solution has higher chance of being selected. There are four types of reproduction schemes as:

1. Roulette wheel selection.
2. Ranking selection.
3. Tournament selection.
4. Elitism.

• After reproduction, next step is crossover. In cross over process, two parents solutions exchange their properties to create two new children solutions. To perform this step, first mating pairs/parents are selected randomly from the mating pool. If there are N no. of solutions then N/2 mating pairs are created. Many crossover operators are available in GA such as:

1. Single point crossover.
2. Two-point crossover.
3. Multi point crossover.
4. Uniform crossover.

• Mutation: Mutation is performed to get diverse solution in a population. Mutation is a bit wise operator i.e. Mutation is performed bit by bit & is performed on each child solution. Mutation changes the bit value from 0 to 1 and vice-versa.

Stalagmite function used in the program is given as:

`f1x = (sin(5.1pix + 0.5))^6;`

`f1y = (sin(5.1piy + 0.5))^6;`

`f2x = exp(-4ln2xx(x-0.0667)^2/0.64);`

`f2y = exp(-4ln2xx(y-0.0667)^2/0.64);`

and `f = (f1x*f1y*f2x*f2y)`

Body of content:

First of all, we need to code the stalagmite function in a separate MATLAB File. One thing we must take care of is that the function name and the file name we are using is the same.

We defined the stalagmite function as a function of input vectors and those input vectors are nothing but the x and y arrays. And then we further specified the required formulas as already covered in the theory part. As GA solver minimizes the function by default, and our main objective is to maximize. So, we need to multiply the whole function by minus sign.

function [f] = stalagmite(input_vector)
x = input_vector(1);
y = input_vector(2);

f1x = (sin((5.1*pi*x) + 0.5))^6;
f1y = (sin((5.1*pi*y) + 0.5))^6;
f2x = exp((-4*log(2))*(x-0.0667)^2/0.64);
f2y = exp((-4*log(2))*(y-0.0667)^2/0.64);

f = -(f1x*f1y*f2x*f2y);

end

Now coming to main program, I have included all the three cases in a single program. first of all, we need to define the search space within which we would like to find the solutions. I this case I have divided 0 to 0.6 into 150 intervals using linspace command and that is saved in x & y array. The number of iterations is set to 50 for first two cases and 200 for the third case. After this we have used meshgrid command so as to generate 2D grid of xx, yy (150x150).

Our next step is to evaluate the stalagmite function using the meshgrid generated and defined function. As both xx and yy are 150x150 matrix, we have to use nested for loop where we stored xx(i, j) and yy(i, j) into input vectors 1 and 2 respectively. Then using these values, we called our stalagmite function whose value will be stored in ‘f’.

clear all
close all
clc

%% defining search space 
x = linspace (0, 0.6, 150); %dividing 0 to 0.6 into 150 intervals
y = linspace (0, 0.6, 150);
itr = 50;  %no. of iterations
itr2 = 200;
%2D mesh
[xx, yy] = meshgrid (x, y);

%% Evaluating the Stalagmite Function
f = zeros(150);
for i = 1:length(xx)
    for j = 1:length(yy)
        input_vector(1) = xx(i,j);  %value of xx at (i,j) will be stored
        input_vector(2) = yy(i,j);
        f(i,j) = stalagmite(input_vector);  %calling stalagmite function
    end
end

Now coming to coding for different cases, case1 is unbounded i.e. we are not setting the limit on where the GA solver should take values of x and y from. Even though we defined our search space from 0 to 0.6 in initial part of our program, but GA solver can take any value to optimize the function.

We started with command- tic which basically initiates the timer & it is complimented by tok command at the end of the loop. This allows us to evaluate the time taken to execute case1 loop. In the for loop we have used standard syntax for GA solver as:

[x, fval] = ga(fun, nvars ) where, fval stands for optimum value of function at input x. we have to enter our objective function name in place of fun and nvars asks for number of variables. After executing this GA solver, we will get optimum value of function and corresponding inputs at which that optimum value is reached. We further stored those inputs in xopt and yopt array by performing this for itr = 50 times using for loop. So xopt and yopt is actually 1x50 vector.

For plotting the data, we need to have 2 plots in single window, for which we have used subplot command to plot stalagmite function f against x and y we have used surfc command. To see the optimum values of function against respective values of x and y we have used plot3 command in which we have specified the marker shape, size & color. The other subplot is Function maxima vs number of iterations for which again we have used plot command.

Case2 is with pre specified upper and lower bound. In case1 the optimum values that we got are randomly distributed and ambiguous. Therefore, to bring a good level of certainty we are limiting the values of x and y where GA-solver can search for functional optima. The coding logic here is similar to case1, only difference being the syntax of GA solver. Here the syntax used for GA solver is:

[x, fval] = ga(fun,nvars,A,b,Aeq,beq,lb,ub) here we are not using any inequality, lb and ub stands for lower bound and upper bound respectively. So here we have defined [0;0], [0.6;0.6] for lb and ub and other unknown fields are left blank. Then again, we have stored optimum inputs into xopt and yopt arrays as in case1. Approach for plotting is same as that of case1 as explained above.

In case3 we are increasing number of iterations and population just to see how well our GA performs for large population size. Here it is logical to note that, due to large population size GA has higher chances at selecting best optimum solution. To implement this case, we used optimoptions command with solver as GA. Below that line in the code we have specified population size as 200. Now in for loop syntax, we have included options after upper bound limit while leaving unknown field blank. Rest of the logic for that part is same as that of earlier cases. Again, we have used same logic to plot the graphs as in earlier cases.

%% case1-unbounded
tic
fopt = zeros;
for i = 1 : itr
    [inputs, fopt(i)] = ga(@stalagmite, 2);
    xopt(i) = inputs(1);
    yopt(i) = inputs(2);
    
end

case1_time = toc %to stop the timer and store the time required.


% Plotting for case1
figure(1)
subplot (2, 1, 1)
colorbar
hold on
surfc(x, y, -f)     %to create 3D surface plot
shading interp
plot3(xopt, yopt, -fopt, 'marker', 'p', 'markersize', 5, 'markerfacecolor', 'c')
title('Unbounded case')
xlabel('x value')
ylabel('y value')
zlabel('f')
subplot(2, 1, 2)
plot(-fopt)
xlabel('No. of Iterations')
ylabel('Function Maxima')

%% case2-with upper and lower bounds
tic
for i = 1 : itr
    [inputs, fopt(i)] = ga(@stalagmite, 2, [], [], [], [], [0;0], [0.6;0.6 ]);
        xopt(i) = inputs(1);
        yopt(i) = inputs(2);       
end
case2_time = toc


% Plotting
figure(2)
subplot (2, 1, 1)
colorbar
hold on
surfc(x, y, -f)
shading interp
plot3(xopt, yopt, -fopt, 'marker', 'p', 'markersize', 5, 'markerfacecolor', 'c')
title('Bounded Inputs')
xlabel('x value')
ylabel('y value')
zlabel('f')
subplot(2, 1, 2)
plot(-fopt)
xlabel('No. of Iterations')
ylabel('Function Maxima')

%% case3-Increasing the number of iterations and population
options = optimoptions('ga');
options = optimoptions(options, 'PopulationSize', 200);     %assigning population size as 200 in our GA.
tic
for i = 1 : itr2
    [inputs, fopt(i)] = ga(@stalagmite, 2, [], [], [], [], [0;0], [0.6;0.6], [], [], options);
    xopt(i) = inputs(1);
    yopt(i) = inputs(2);
end
case3_time = toc



% Plotting
figure(3)
subplot(2,1,1)
colorbar
hold on
surfc(x, y, -f)
shading interp
plot3(xopt, yopt, -fopt, 'marker', 'p', 'markersize', 5, 'markerfacecolor', 'c')
title('With higher No. of interations')
xlabel('x value')
ylabel('y value')
zlabel('f')
subplot(2,1,2)
plot(-fopt)
xlabel('No. of Iterations')
ylabel('Function Maxima')

Results:

The plots that we got are as:

command window output:

Output from the workspace window:

Conclusions:

• We get more clarity about optimum values when we limit our search space.
• The certainty in the solution has further improved for higher population size with considerable increase in computation time as well.

Errors:

I got the error where it warned that fopt appears to change size on every loop iteration.

To correct this error I need to wrote f = zeros(150) before stalagmite function for loop, here according to the syntax, 150 is included for the fact that f is 150x150 matrix.

We have also included fopt = zeros just before the for loop in case1. That sorted out our problem.

References:

• https://www.iitg.ac.in/rkbc/CE602/CE602/Genetic%20Algorithms.pdf
• https://en.wikipedia.org/wiki/Genetic_algorithm
• https://in.mathworks.com/help/gads/how-the-genetic-algorithm-works.html
• https://nptel.ac.in/courses/112/105/112105235/

hh