MATLAB Project: Genetic Algorithm

Aim:

To write a code in MATLAB to optimise the stalagmite function and find the global maxima of the function.

Theory:

A genetic algorithm (GA) is a method for solving both constrained and unconstrained optimization problems based on a natural selection process that mimics biological evolution [1]. That is, GA use principles of natural Selection and Genetics in timetabling and scheduling problems to optimise it to a better solution. GAs are a subset of much larger branch of computation known as Evolutionary Computation. It is frequently used to find the optimal or near-optimal solutions to difficult problems which otherwise would take a lifetime to solve [2].

Optimization is the process of making something better which is basically finding the values of the inputs in such a way that we get the best output values. In mathematical terms, the best output refers to the maximising or minimising one or more objective functions by varying the input parameters. The set of all possible solutions or values which the inputs can take make up the search space (population). In the search space, lies a point or a set of points which gives the optimal solution. The aim of optimisation is to find that point or set of points in that search space.

In MATLAB, the major function of the genetic algorithm is to find the global minima for highly nonlinear problems. It is only possible to find minima and not maxima points in MATLAB. The GA repeatedly modifies the population of individual solutions. At each step, the GA randomly selects individuals from the current population and uses them as parents to produce the children for the next generation which is known as recombination and mutation and this process is repeated over a various generation. Each individual is assigned a fitness value (based on the objective function value) and the fitter individuals are given a higher chance to mate and yield more fitter individuals. Over successive generations, the population evolves toward a optimal solution.

We can apply the genetic algorithm to solve problems that are not well suited for standard optimization, including problems in which the objective functions is discontinuous, nondifferentiable, stochastic, or highly nonlinear [1].  

Stalagmite function
A stalagmite is a type of rock formation that rises from the floor of a cave due to the accumulation of material deposited on the floor from the ceiling drippings. 

A similar geometrical shape can be created by a mathematical function. This function when applied to two variables in two dimensions would yield the stalagmite. This function can be mathematically expressed as below:

`f(x,y)=f_(1,x)*f_(2,x)*f_(1,y)*f_(2,y)`

where,
`f_(1,x)=[sin(5.1pix+0.5)]^6`
`f_(1,y)=[sin(5.1piy+0.5)]^6`
`f_(2,x)=exp[-4ln(2)(x-0.0667)^2/0.64]`
`f_(2,y)=exp[-4ln(2)(y-0.0667)^2/0.64]`

Objectives:

1. Understanding the concept of Genetic Algorithm and syntax for ga in MATLAB.
2. Writing a code to optimise the stalagmite function and finding the global maxima of the function.
3. Plotting the graphs for all different studies and its evaluation.

Program to optimise the Stalagmite function and finding the global maxima:

The programming is done in MATLAB.

The program initially plots the stalagmite function. The genetic algorithm in MATLAB is made use to find the global minimum (here maximum, how: is explained after function definition of the function 'stalagmite_func'.

clear all
close all
clc

�fining search space
x = linspace(0,0.6,150);
y = linspace(0,0.6,150);
cases = 50; % used for 3 studies for 50 iterations

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

% Evaluating 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_func(input_vector);
    end
end

% to change the negative term in the stalagmite function
F = -f; % in order to remove the -ve sign of f_s in the function 'stalagmite_func'
figure (1) % creating figure
surfc (xx,yy,F);
colorbar;
xlabel('x value');
ylabel('y value');
zlabel('Stalagmite Function Values');
shading interp; % in order to remove nodal lines
grid on;
grid minor;

% Study 1: Statistical behaviour
tic % for calculating the time taken for study1 (50 iterations)
for i=1:cases
    
[coordinates1,min_val1(i)] = ga(@stalagmite_func,2); % find coordinates and find the min value using GA
x_coordinate1(i) = coordinates1(1); % finding x co-ordinate of each max point found by GA
y_coordinate1(i) = coordinates1(2); % finding y coordinate of each max point found by GA

end
t1 = toc;

figure (2) % creating figure
subplot(6,1,[1:3]) % for plotting stalagmite function
hold on;
surfc(x,y,-f) % -ve sign for f to remove -ve sign of f_s in 'stalagmite_func' function
xlabel('x value');
ylabel('y value');
zlabel('Function Values');
shading interp; % to remove nodal lines
max_val1 = -min_val1; %finding max value which is -ve of min value
plot3(x_coordinate1,y_coordinate1,max_val1,'color','r','marker','*') % for marking each max points found by GA on stalagmite function plot
grid on;
grid minor;
title('Study 1');

subplot(6,1,[5,6]) % for plotting each max value obtained by GA for each iteration
plot(max_val1,'marker','*')
xlabel('Iterations');
ylabel('Function maximum');
axis([0 50 0 1]);
grid on;
grid minor;

% Study 2: Statistical behaviour of GA with upper and lower bounds
tic
for i=1:cases
    
[coordinates2,min_val2(i)] = ga(@stalagmite_func,2,[],[],[],[],[0;0],[0.6;0.6]);
x_coordinate2(i) = coordinates2(1);
y_coordinate2(i) = coordinates2(2);

end
t2 = toc;

figure (3)
subplot(6,1,[1:3]) % plotting stalagmite function
hold on;
surfc(x,y,-f)
xlabel('x value');
ylabel('y value');
zlabel('Function Values');
shading interp; % to remove nodal lines
max_val2 = -min_val2; % finding max value
plot3(x_coordinate2,y_coordinate2,max_val2,'color','r','marker','*') % for plotting max points found by GA on the stalagmite surface curve
grid on;
grid minor;
title('Study 2');

subplot(6,1,[5,6]) % for plotting max points found by GA for each iteration
plot(max_val2,'marker','*')
xlabel('Iterations');
ylabel('Function maximum');
axis([0 50 0 1]);
grid on;
grid minor;

% Study 3: Statistical behaviour of GA by increasing the no. of iterations
options = optimoptions('ga')
options = optimoptions(options,'populationsize',400);
tic
for i=1:cases

[coordinates3,min_val3(i)] = ga(@stalagmite_func,2,[],[],[],[],[0;0],[0.6;0.6],[],[],options);
x_coordinate3(i) = coordinates3(1);
y_coordinate3(i) = coordinates3(2);

end
t3 = toc;

figure (4)
subplot(9,1,[1:3]) % for plotting stalagmite function 
hold on;
surfc(x,y,-f)
xlabel('x value');
ylabel('y value');
zlabel('Function Values');
shading interp;
max_val3 = -min_val3;
plot3(x_coordinate3,y_coordinate3,max_val3,'color','r','marker','*') % for plotting max points found by GA
grid on;
grid minor;
title('Study 3');

subplot(9,1,[5,6]) % for plotting max value found by GA for each iteration
plot(max_val3,'marker','*')
xlabel('Iterations');
ylabel('Function maximum');
axis([0 50 0 1]);
grid on;
grid minor;

subplot(9,1,[8,9]) % for plotting max value found by GA for each iteration (zoomed view)
plot(max_val3,'marker','*')
xlabel('Iterations');
ylabel('Function maximum');
grid on;
grid minor;

txt = [' TIME TAKEN' newline 'For study 1:  ' num2str(t1) ' s' newline 'For study 2:  ' num2str(t2) ' s' newline 'For study 3:  ' num2str(t3) ' s']

Each step of the program is explained with comments.

The x and y values vary from 0 to 0.6 and it is equally spaced into 150 elements using the linspace command.

The function definition of the function 'stalagmite_func' is provided below:

function [f_s] = stalagmite_func(input_vector)

        
f_1x = (sin(5.1*pi*input_vector(1)+0.5))^6;
f_1y = (sin(5.1*pi*input_vector(2)+0.5))^6;
f_2x = exp(-4*log(2)*((input_vector(1)-0.0667)^2/0.64));
f_2y = exp(-4*log(2)*((input_vector(2)-0.0667)^2/0.64));
        
f_s = -1*f_1x*f_2x*f_1y*f_2y;

% Note
% Genetic algorithm can find only the minimum optimised value of the
% function. But here, our task is to find the maximum, which is not
% possible with GA in MATLAB. So, in order to find the maximum, the
% variable, f_s is given -ve sign, so that the GA can find the minimum
% value which is actually the maximum value of the vector f_s.
        
        

The meshgrid command is used to transform one dimensional array into two-dimensional array, which can be used for evaluation of the functions of two variables and 3D surface plots.
The ga command is used with different parameters for each study in order to improve the quality of the answer as well as accuracy.

Output of the program:


The output plot of the stalagmite function is provided above.

Three studies are conducted on the genetic algorithm so as to find the what's happening for each syntax of genetic algorithm and to check whether it is giving same result as global maximum for 50 iterations.

Study 1: Statistical behaviour
The basic syntax of genetic algorithm '[x,fval] = ga(fun,nvars)', which finds the local unconstrained minimum to the objective function, fun and nvars is the number of design variable of fun. Here, GA does not really find the global minimum rather local minimum. Only basic parameter about the number of variables of the function is mentioned. The GA is not provided any other parameters than this which would result in low accuracy or wrong answer as the GA is taking 50 random points which could be out of the range leading to error.

The output of the study 1 is shown above. Since no upper and lower is provided for the variable x and y, we could see that, the points which was found as maximum using ga is different for different iteration, that is, it is giving different value at each time. This reveals that, there exist a requirement for us to add parameters to ga command for making the optimization process smoother which leads to the study 2.

Study 2: Statistical behaviour with upper and lower bounds
We know that our x and y values ranges from 0 to 0.6. Therefore, in order to find the point at which the function is maximum, there is no need to randomly select the point outside the domain area. The upper and lower limit for the ga command is provided [x,fval] = ga(fun,nvars,A,b,Aeq,beq,lb,ub). Since there does not exist any linear inequalities exist, we can replace A, b, Aeq and beq with [] and provide lb with 0 value and ub with 0.6 as [0;0] and [0.6;0.6] for both x and y variables. Hence, the values of x and y will be within the defined search space in inputs.

After executing the program, we get the output as:

From the plot, it is clear that, with the introduction of lower and upper bounds, the GA improved its performance resulting in only 3 coordinate positions to show the maximum function value.

Study 3: Statistical behaviour by increasing the no. of iterations with upper and lower bounds
For the other two studies, the population was set as 50 by default. So, in order to increase the number of iterations for improving the solution, the ga command was updated. [x,fval] = ga(fun,nvars,A,b,Aeq,beq,lb,ub,nonlcon,IntCon,options) where the x, fval, A, b, Aeq, beq, lb and ub are same as Study 1 and 2. nonlcon refers to the nonlinear constraints, specified as a function handle or function name and IntCon refers to the integer variable, specified as a vector of positive integers. Since we don't have both, we can replace both with []. options refer to the optimization options, specified as the output of optimoptions or a structure. Using options command, the population size is designated as 400 which could definitely improve the solution but could lead to increase in the time taken for calculating it.

The output plot of the Study 3 is provided below:

From the plot, it is obvious that there is only one point with maximum value for the stalagmite function. In order to clarify it, the third plot from the study 3 can be viewed where the max value will only vary from 0.99999985 to 0.9999995 which can be approximated to 1 which is seen on second plot.

The time taken for different study is found using the tic toc command and below shows the time taken for each study:

With this result, we can conclude that the time would increases with increase in accuracy of the result which is obvious. The third study has taken 19.8732 seconds because its population is 400 while the other two studies have only 50 (by default). 

Conclusion:

The MATLAB program to optimise the stalagmite function and finding the global maxima of the function has been successfully created. the concept of Genetic Algorithm and syntax for ga in MATLAB has been understood properly. Three studies were conducted and graphs for each study has been plotted and evaluated. Also, the time taken for each study has also been found.

Errors:

1. No particular errors were there except typo mistakes. 
2. For study 3, the population size was increased till a particular value (here, 400) so that the required output was obtained.

References:
[1] https://in.mathworks.com/discovery/genetic-algorithm.html
[2] https://www.tutorialspoint.com/genetic_algorithms/genetic_algorithms_quick_guide.htm
[3] https://in.mathworks.com/help/gads/ga.html
[4] https://in.mathworks.com/help/gads/genetic-algorithm.html
[5] https://www.researchgate.net/profile/Khalid_Al-Hussaini/post/How_to_write_genetic_algorithms_codes_in_Matlab/attachment/59d6339b79197b8077991473/AS:375438057721856@1466522816038/download/Genetic-Algorithm-Implementation-Using-Matlab.pdf
[6] https://www.youtube.com/watch?v=esGZtl3gIlY