Optimizing the Stalagmite function using Genetic algorithm

Code for optimizing a function to give a global maximum using genetic algorithm.

%Stalagmite Mathematical Model
clc
clear all
close all

%Setting up the search space
x=linspace(0,0.6,150);
y=linspace(0,0.6,150);

%creating a 2D mesh
[xx yy]=meshgrid(x,y);

%Calling the stalagmite function
for i=1:length(xx)
    for j=1:length(yy)
        input_vector(1)=xx(i,j);
        input_vector(2)=yy(i,j);
        f(i,j)=stalagmite(input_vector);
    end
end

%Plotting the function
surfc(xx,yy,f)
shading interp

%using generic algorithm 
[results fval] = ga(@stalagmite,2);
maximum = -fval

Stalagmite Function:

function r = stalagmite(inputvec)
t1 = sin ((5.1*pi*inputvec(1))+0.5);
t2 = sin ((5.1*pi*inputvec(2))+0.5);
t3 = (-4*log(2))/0.64;
t4 = inputvec(1) - 0.0667;
t5 = inputvec(2) - 0.0667;

f1x = t1^6;
f1y = t2^6;
f2x = exp(t3*t4^2);
f2y = exp(t3*t5^2);
f = (f1x*f1y*f2x*f2y);
r = -f;
end

Explanation of the code:

1. We create 2 arrays x and y and with the help of meshgrid command, we convert them into 2D arrays. The meshgrid command repeats all the row elements of the first array into subsequent rows while for the second array, it repeats the elements of the row as columns. These two arrays will allow us to traverse the 2D plane on which we will test our function.
2. We use two for loops to access each and every element of the plane and store the x coordinate and y coordinate in an array called input_vector and send the location points to stalagmite function which uses them to calculate the value of function according to the formula given below
3. Since we finally intend to find the global maxima, we return the negative value of function.
4. We plot the function as a surface using the surf command and use shading interp to disable grid lines and shade the surface.
5. We finally call the genetic algorithm and pass it the stalagmite function as it is the function to be optimized. We also indicate the number of variables in the stalagmite function(2).
6. The ga returns the coordinates of the maximum value and also the maximum value itself, which are stored in array results and fval.
7. Since we had given our function a negative value, the GA also returns negative maxima(which is the mimina according to ga). So we finally convert it back by negating it again and storing it in \'maximum.

Result:

However, the above code gives different results for each execution of the code and hence this is not our final code. We will improve and produce a stable code later.

Genetic Algorithm:

The genetic algorithm is an optimization technique used to optimise functions based on the philosophy of Darwin\'s theory of evolution. It is used to solve non-linear and non-differentiable optimization problems.

The algorithm works by first starting with a random set of population and tests these random data points against the objective function. Based on the results a set of better points is selected. These points will serve as both, Next-generation points as well as parents to the next generation. Using them some new points are created by Mutation and Cross over.

The above process is repeated until the algorithm converges. A Genetic algorithm can converge by several means. The most popular ones include: Iteration until a fixed number of generations or until the best value of the objective function is not changing.

Some other ways a generic algorithm stops are:

• Time Limit: Algorithm stops after it has run for a specific time period.
• Stall Time limit: Algorithm stops if it fails to improve the function value for a specified time limit known as stall time
• Fitness Limit: Algorithm stops if the value of the function comes out to be less than or equal to the fitness limit.

Genetic Algorithm Syntax:

%Study 1 - Statistical Behaviour

clear all
close all
clc

x =linspace(0,0.6,150);
y =linspace(0,0.6,150);

[xx yy] = meshgrid(x,y);


for i = 1:length(xx)
    for j = 1:length(yy)
        input_vector(1) = xx(i,j);
        input_vector(2) = yy(i,j);
        f(i,j) = stalagmite(input_vector);
    end
end

iterations = 50;
tic
for i=1:iterations
    [input fopt(i)] = ga(@stalagmite, 2);
    xopt(i) = input(1);
    yopt(i) = input(2);
end
study1_time = toc
figure(1)
subplot(2,1,1)
hold on
surfc(x,y,f)
shading interp
plot3(xopt,yopt,fopt,\'marker\',\'o\',\'markersize\',5,\'markerfacecolor\',\'r\')
title(\'Unbounded inputs\')
subplot(2,1,2);
plot(fopt)
xlabel(\'Iterations\')
ylabel(\'Function minimum\')

In this study, we execute the genetic algorithm several times and observe the output of each iteration by plotting them on graph.

Explanation:

1. Here we start similarly and call the stalagmite function to provide us the value of function for various locations.
2. We then set \'iterations\' to 50 as that is the number of times we want to run and compare our GA. For multiple iterations, we use the for loop and call ga in each iteration. As Ga return the location coordinates and optimum value of function, we store them in arrays xopt/yopt and fopt.
3. We then open a new figure window and use the subplot command to have multiple plots inside the same figure window.
4. We then create the surface plot using surface command and then use plot3 command which will allow us to plot 3D graph of x and y coordinates along with the optimum values for all the 50 iterations of GA.
5. The plot3 command will give us an idea of the variation in the optimum value provided by GA. But to get a better idea of the variation, we create another plot. For this we again use subplot command and plot all the values of \'fopt\' using plot command.
6. Suitable labels and title is also provided. The tic toc command provides us with the time that the iteration of a certain part of a program took. We have used it on the genetic algorithm for loop to calculate the time 50 executions of GA take.

Result:

The following results were achieved after the execution:

1. We can see that the values predicted by GA vary greatly and are not consistent. This is because GA uses random sampling and hence every time comes up with a new optimum value.

2. Also the x and y coordinate go beyond our range of 0 to 0.6. Because of such vast variation in results, they are not usable to us.

Study 2: Statistical behaviour(With defined upper and lower bounds)

In this study, we define the upper and lower bounds so that x and y values remain restricted (between 0 and 1) and this, in turn, will restrict our results from GA.

%Study 2
tic
for i=1:iterations
    [input fopt(i)] = ga(@stalagmite, 2,[],[],[],[],[0;0],[1;1]);
    xopt(i) = input(1);
    yopt(i) = input(2);
end
study2_time = toc
figure(2)
subplot(2,1,1)
hold on
surfc(x,y,f)
shading interp
plot3(xopt,yopt,fopt,\'marker\',\'o\',\'markersize\',5,\'markerfacecolor\',\'r\')
title(\'Bounded inputs\')
subplot(2,1,2);
plot(fopt)
xlabel(\'Iterations\')
ylabel(\'Function minimum\')
 

The GA syntax is used to provide the upper and lower bounds. Like the last study, we again execute the GA 50 times and store the location and optimum values in the respective arrays.

Result:

We can see, that the variation in optimum values provided by several executions of GA has reduced drastically.

Also majority of results are lying in a small range of about 0.2 except for a few peaks.

Thus by reducing the search space by providing lower and upper bounds, we are able to greatly improve the accuracy of GA.

Study 3:

For the third study, apart from giving the bounds to x and y, we also increase the population size for the random sampling of GA. By default, this size is set to 50, but as we increase it, the accuracy of our results increases, at the cost of time consumed for processing.

%Study 3
options = optimoptions(\'ga\');
options = optimoptions(options,\'PopulationSize\', 150)

tic
for i=1:iterations
    [input fopt(i)] = ga(@stalagmite, 2,[],[],[],[],[0;0],[1;1],[],[],options);
    xopt(i) = input(1);
    yopt(i) = input(2);
end
study3_time = toc
figure(3)
subplot(2,1,1)
hold on
surfc(x,y,f)
shading interp
plot3(xopt,yopt,fopt,\'marker\',\'o\',\'markersize\',5,\'markerfacecolor\',\'r\')
title(\'Bounded inputs\')
subplot(2,1,2);
plot(fopt)
xlabel(\'Iterations\')
ylabel(\'Function minimum\')

In this code, we the \'options\' to change our population size. The \'options\' feature allows us to change the optimization parameters of the GA, and hence we have used the appropriate syntax to call the GA.

Result:

The rage of results provided by GA reduces even more now(3). Also, all the results lie in a very small range. Only on two occasions, the ga did not provide us with accurate results.

This problem can be solved by increasing the population size even more.

The following results were achieved after the population size was increased to 250 from 150.

As we can see from the plot, the variation in the optimised value predicted by GA his highly stable and consistent (upto 6 decimal places)

Another trial was done by keeping the population size as 200, but the search space was further refined to exact values of 0 to 0.6. The following results were achieved:

Hence we can finally conclude that for getting accurate and stable results from genetic algorithm, we should define the search space and also increase the population size.

Code to optimize the stalagmite function with GA:

The following code is our final code to provide us with stable results for optimizing the stalagmite function.

%Stalagmite Mathematical Model
clc
clear all
close all

%Setting up the search space
x=linspace(0,0.6,150);
y=linspace(0,0.6,150);

%creating a 2D mesh
[xx yy]=meshgrid(x,y);

%Calling the stalagmite function
for i=1:length(x)
    for j=1:length(y)
        input_vector(1)=x(j);
        input_vector(2)=y(i);
        f(i,j)=stalagmite(input_vector);
    end
end

%Plotting the function
surfc(x,y,f)
shading interp

%Initiating the genetic algorithm

options = optimoptions(\'ga\');
options = optimoptions(options,\'PopulationSize\', 200)

[input fopt] = ga(@stalagmite, 2,[],[],[],[],[0;0],[0.6;0.6],[],[],options);
maximum = -fopt

Plot and Result: