Genetic Algorithm in MATLAB to find GLOBAL Maxima using Stalagmite function

AIM : To Optimize the Stalagmite function and to find Global Maxima of the function using Genetic Algorithm in MATLAB.

OBJECTIVES : The objective of this report is to optimize the Stalagmite function by coding the program with help of Genetic algorithm in MATLAB and explain the concept of Genetic algorithm and also explain the syntax of “ga” in MATLAB .

STALAGMITE :

A Stalagmite is a type of rock formation that rises from the floor of a cave due to the accumulation of the material deposited on the floor from ceiling drippings. Stalagmites are typically composed of calcium carbonate, but may consist of lava, mud, peat, pitch, sand, sinter and amberat (crystallized urine of pack rats.)

Since our function gives a shape of stalagmite it is named as stalagmite function . Stalagmite function is the mathematical model that helps us to generate minimas and maximas. It is 3-Dimensional function and is designate to generate local maximas.

THEORY OF GENETIC ALGIORITHM:

A Genetic Algorithm is a search heuristic that is inspired by Charles Darwin’s Theory of natural evolution. This algorithm reflects the process of natural selection where the fittest individual is selected for reproduction in order to produce offspring of the next generation.

In computer science, operations research , a genetic algorithm are commonly used to generate high-quality solutions to optimation and search problems by relying on biologically inspired operators such as mutations, crossover and selection . John Holland introduced genetic algorithm in 1960 based on the concept of Darwin’s theory Od evolution, and his student David E.Goldberg further extended GA in 1989.

Notion of Natural selection

The process of natural selection starts with the selection of fittest individual from a population. They produce offspring which inherit the characteristics of the parents and will be added to the next generation .If parents have better fitness ,their offspring’s will be better than parents and have better chance of survival. This process keeps on iterating until it generates the fittest individual.

This notion can be applied for search problem .We considers set of solution for a problem and select the set of best one out of them.

Five set are consider in a Genetic Algorithm.

1. Initial population
2. Fitness function
3. Selection
4. Crossover
5. Mutation

• Initial Population :

The process begins with the set of individual which is called as population. Each individual is solution to the problem you want to solve.

An individual is characterized by a set of parameters (variables) known as Genes .Genes are joined into strings to form a Chromosome (solution).

In genetic algorithm, the set of genes of individual is represented using string, in term of an alphabet; usually binary values are used (strings of 1s and 0s). We say that we encode the genes in chromosome.

• FITNESS FUNCTION :

The fitness function determines how fit an individual is (the ability of an individual to compete with individuals). It gives the fitness score to each individual .The probability that an individual will be selected for reproduction is based on it fitness score.

• SELECTION :

The idea of selection phase is to select the fittest individuals and let them pass on their genes to next generation.

Two pairs of individuals (parents) are selected based on fitness scores. Individual with high fitness have more chance to be selected for reproduction.

• CROSSOVER :

Crossover is the most significant phase in a genetic algorithm .For each pair of parents to be mated a crossover point is chosen at random from within the genes.

For example, consider the crossover points to be 3 as shown below.

Offspring’s are created by exchanging the genes of parents among themselves until the crossover point is reached.

The new offspring are added to the population

• MUTATION :

In certain new offspring formed, some off their genes are subjected to a mutation with a low random probability. This implies that some of the bots in the bit string can be flipped.

Mutation occurs to maintain diversity within the population and prevent premature convergence.

• TERMINATION :

The algorithm terminates if the population has converged. Then it is said that the genetic algorithm has provided the set of solution to our problem.

Using this approach of Genetic Algorithm we are going to optimize the stalagmite function and find the optimized value of stalagmite function.

GOVERNING EQUATION:

To generate the stalagmite function following function is used.

`f_1 (x)= [sin⁡(5.1*π*x+0.5 ) ]^6`

`f_1 (y)= [sin⁡(5.1*π*y+0.5 ) ]^6`

`f_2 (x)=exp (-4.log2 .(x - 0.0667)^2)/0.64`

`f_2 (y)=exp (-4.log2 .(y - 0.0667)^2)/0.64`

`f(x,y)= f_1 (x).f_1 (y).f_2 (x) 〖.f〗_2 (y)`

Maximizing f (x,y) is the same as minimizing - f (x,y)

PROGRAM:

CODE FOR FINDING GLOBAL MAXIMA :-

clear all 
close all
clc

% Defining our Search Space
x = linspace(0,0.6,150);
y = linspace(0,0.6,150);

% Creating 2-d Array
[xx,yy] = meshgrid(x,y)

% Calculating Fitting Function
for i = 1:length(xx);
    for j = 1:length(yy);
        input(1) = xx(i,j);
        input(2) = yy(i,j);
        f(i,j) = stalagmite_func(input);
    end
end

% Converting Fitting Function to change the graph upwards
F = -f;

% Plotting Stalagmite Function
surface(xx,yy,F);
% shading interp;
hold on;

% Performing the Genetic Algorithm
[output , fval] = ga(@stalagmite_func,2,[],[],[],[],[0,0],[0.6,0.6])

% Plotting Global Maxima
x_max = output(:,1);
y_max = output(:,2);
f_max = -fval;
plot3(x_max,y_max,f_max,'marker','o','markersize',5,'markerfacecolor','r');
xlabel ('X-value');
ylabel ('y-value');
title('GLOBAL MAXIMA');
legend('STALAGMITE','GLOBAL MAXIMA')

• Here the command “Close all”, “Clear all”, “ Clc” are use to close any previous plot and to clear the work space and command window respectively.
• Then the value of ‘x’ and ‘y’ are entered using ‘linspace ()’ command. This is our search space, that means for which value of ‘x’ and ‘y’ that have defined in ‘linspace ()’ are we getting an optimized value for stalagmite function.
• Now using “meshgrid” command this 1-D array is converted into 2-D array. These 2-Dimensional arrays are named as ‘xx’ &’yy’ in the code.
• Now to calculate the stalagmite function we have to use “FOR LOOP”, and in this loop we have define ‘i’ & ’j’ as each value of ‘x’ and ‘y’ of 2-dimensional array.
• Now for each value of ‘x’ & ‘y’ we are creating an array named as “inputs” and thus this input will be row vector i.e. 2 column matrix.
• Now first column of input vector is defined as value of ‘xx’ matrix and second column of input vector are defined as ‘yy’ matrix values.
• Now using these values the stalagmite function is calculated. Function for stalagmite is entered separately as a function code using the equations as mentioned above.

FUNCTION CODE:

• Now this function takes the ‘input’ values as the input for the function and returns the value of ‘f’.
• Now this above function is called in the “FOR LOOP” for different value of ‘x’ and ‘y’ and different value of ‘f’ are calculated.
• These two “FOR LOOP” generates the data for the plot that we need to generate a stalagmite graph.
• Then ‘surface ()’ command is used which creates a primitive, 3-dimensional surface plot .If the command is entered as “surface(X,Y,Z)” Then the function plots the value in matrix Z as height above the grid in the X-Y plane defined by X & Y .The color of the surface varies according to the height specified by Z. Here in this code (X,Y,Z) is replaced by (xx,yy,F) respectively.
• Then “shading interp” command is used to remove the grid that are form on plot ,which will give us a smooth and better looking plot.
• Now to optimize the above function we use a inbuilt MATLAB command called “ga ()” .This command perform the genetic algorithm on the function which we have define and gives us output as optimized results. It will automatically perform the genetic algorithm calculation and give us the optimized solution for the fitting function as output.
• Syntax for “ga ()” command is as follows:

X = ga(fun,nvars) find the local unconstrained minimum x , to the objective function (fun); and nvars is the dimension (no. of design variable) of fun.

X = ga(fun,nvars,A,b,Aeq,beq,lb,ub) defines the set of lower and upper bounds on the design variables x, so that the solution is found in the range ub,{ Set A = [] and b = [] if no linear inequalities exist } . { Set Aeq = [] and beq = [] if no linear equalities exist.}

[x , fval] = ga(__), for any previous input argument , also return fval , the value of the fitness function at x.

• The above command are syntaxes for the ”ga ()” command used un above code. Now the 1^st syntax is the basic algorithm command used where ‘ga ()’ takes the inputs as fitting equation that has to be evaluated and the number of variable present in the fitting function .So each time this command is run it take the random value for the function variables initially and then use those value for further calculations.
• The next syntax of ‘ga ()’ allows us to provide the bounds to the ‘ga ()’ command that will limit it, to choose the value within the specified value for the function variable. We can also give linear inequalities and equalities as input and ‘ga ()’ will calculate accordingly to these inequalities and equalities.
• The 3^rd syntax of ‘ga ()’ gives the value of ‘x’ as well as function minimum value. So we can get value of ‘x’ where the function is minimum and that minimum value of function is also calculated and given as output.
• Now drawback of ‘ga ()’ command is, it always calculates the value for ‘x’ where the function is minimum ,thus gives us the global minima. But we need to calculate global maxima in our code ,So to achieve that we are minimizing the fitting function with a
  “ – “ (minus) symbol, because Maximizing f(x,y) is the same as minimizing –f(x,y).
• That is the reason we have put” – “ (minus) symbol in the equation of “f” in the function code.
• Now since the value of “f” in the function code will be negative, So our plot will be downwards, in order to make it upwards, after the ‘FOR LOOP’ the value of ‘f’ are again multiplied with “ – “ (minus) symbol, thus giving us positive value and making the plot upwards. These new value of are assigned to a variable “F” this is done only to get the graph upwards so we get the better look as it won’t affect the ‘ga ()’ command and its calculations.
• Now after the “shading interp” command ,the ‘ga ()’ is entered in the code and the output of ‘ga ()’ are mentioned as ‘output’ and ‘fval’ i.e. for which value of ‘x’ and ‘y’ which are ’output’ in the code, the function value which is “ fval “ in the code is minimum and it will give value of this function’s minimum value also.
• Then for the inputs of ‘ga () ‘ ,the fitting function which is named as “stalagmite_func” which is called with the function handler, and the no. of variable in the fitting code are entered and the bounded value (limits of ’x’ & ‘y’)are entered so that ‘ga ()’ only gives value between these bonds are entered.
• Since we have no linear inequalities and equalities so it is just entered “[]”in the ’ga ()’ command.

• Now once we run the code we will get the row vector for “output” and the value for “fval” in command window.1^st value of outputs vector is the value of ‘x’ and 2^nd value of the outputs vector is the value of ‘y’ for which value of fitting function is minimum.
• Now to plot the Global Maxima point in the graph we have define ‘x’ value as “x_max” which is 1^st value for output vector & ‘y’ value as “y_max” which is 2^nd value of outputs vector.
• As mention above Maximizing f(x,y) is the same as Minimizing –f(x,y), so value of ‘fval’ is the global maxima of the function “+f”, Since now our graph will appear upwards.
• Then “plot3()” as we need to plot point in 3-D Graph ,Global Maxima is plotted using “x_max” & “y_max” .

• In figure we can observe that there are many local Maximas in the plot and the plot gives us the global maxima of the function.

To demonstrate how “ga ()” works we have done 3 type of studies and coded it. These studies show how to use genetic algorithm to obtain best results.

STUDY _ 1 : STATISTICAL BEHAVIOUR

CODE :

% Statistical Behaviour

clear all 
close all
clc

% Defining our Search Space
x = linspace(0,0.6,150);
y = linspace(0,0.6,150);
num_case = 50;

% Creating 2-d Array
[xx,yy] = meshgrid(x,y)

% Calculating Fitting Function
for i = 1:length(xx);
    for j = 1:length(yy);
        input(1) = xx(i,j);
        input(2) = yy(i,j);
        f(i,j) = stalagmite_func(input);
    end
end

% Converting Fitting Function to change the graph upwards
F = -f;


% Performing the Genetic Algorithm
tic
for k =1:num_case
[inputs , fopt(k)] = ga(@stalagmite_func,2);
xopt(k) = inputs(1);
yopt(k) = inputs(2);
end

study_time = toc

% Plotting Global Maxima And Function Values
figure(1)
subplot(2,1,1)
hold on
surfc(x,y,F)
shading interp
plot3(xopt,yopt,-fopt,'marker','o','markersize',5,'markerfacecolor','r');
xlabel ('X-value');
ylabel ('y-value');
title('Unbounded Inputs');
legend('STALAGMITE','GLOBAL MAXIMA')
subplot(2,1,2)
hold on
plot(fopt)
xlabel ('Iterations');
ylabel ('Function Maxima');
title('Variation in F');
legend('Change in F value')

EXPLANATION :

• This study shows us the statistical behavior of genetic algorithm in MATLAB.
• The code is similar to that of above code, the only change here we are performing “ga ()” for some number of iteration and see how the value are changing for each iteration.
• In the above code, we have defined ‘x’ & ‘y’ values and then we have decide the number of times the ‘ga ()’ has to be performed as a variable “num_case” and we have given it a value.
• Then using ‘x’ and ‘y’ we create 2-D array using “mesh grid()” command.
• Now using the For Loop value of “f” are determined, it is same as the above code used in this “FOR LOOP”.
• Now since the value of ‘f’ is negative and we get the downward graph ,so in order to get upward graph ,we have multiplied ‘f’ with “-“(minus) symbol and assigned it to variable “F”.
• Now here another FOR LOOP is used to evaluate “ga ()” for 50 iterations we have defined “l” as the number of iteration and then we have entered the “ga ()" function with the normal syntax that we have mention in previous code. Only the fitting function and number of variable are entered as inputs for the “ga ()” function.
• The 1^st value of input is assigned to a variable “xopt” and the 2^nd variable “yopt” . Every time for loop is executed it give the new value opf [“fopt”,”xopt”,”yopt”]
• Then the “fopt” value are multiplied with “-“ (minus) symbol and assigned to the variable “Fopt”.This is done in order to point the global maxima since our graph will be upward . this is only done for plotting purpose.
• Now for this for loop there is a command used “Tic and Tok” this command helps to know the time taken to calculate this “FOR LOOP” ,”tic” is entered at start of “FOR LOOP” & “toc” is entered at end of “FOR LOOP” and for the “toc” a variable is assigned so the time in calculated under that variable.
• Then the figure () command is entered to open a new figure window, then the subplot () command is used to open the first graph in the figure window. Subplot (m , n , p) divides the current figure into an M-by-N grid and creates axes in the position specified by p. Accordingly the value are entered in the code for “subplot()”
• Then the plotting is done using the “subplot() “ command and then “Shading interp “ is entered in the code to remove grid on the graph in the 1^st “subplot”
• Then “plot3 ()” is used to generate the point in a 3-D graph, this command is used to plot the global maxima point. This help us how we are getting the global maxima each and every times “a ()” function is run.
• Then “title ()” command is use to give a title to the 1^st plot, “xlabel” &”ylabel” command are entered to name x-axis and y-axis respectively.
• Next 2^nd subplot command is entered to get 2^nd graph between the optimize “f” value obtained from”ga ()” function and the number of iteration that “ga ()” is run.
• Then “plot()” command is used to plot a graph between “Fopt” values and number of iterations.
• Thus in the same figure windows we will get two graph using “subplot ()” command.

• If we run the code we get the above results, here we can observed that we have got different optimized values of function “f” each the “ga ()” was run. We have got 50 different value of “f” when 50 times “ga ()” was run. So in this way we can’t use “ga ()” for our calculation in reality.
• We can observe the 2^nd graph in the figure that how values of “f” are changing for each iteration.
• Genetic Algorithm always take random population sample and then perform the mutation and crossover process to get a new generation. Due to this random inputs taken by genetic algorithm we get the different value when we run the “ga ()” each time .Here genetic algorithm does not know what is our search space i.e. limits of ”x” and “y” .
• To get better optimized value of “f” that are within limits of “x” and “y” defined by us, we have used a modified “ga ()” syntax, which we will use in study 2.

• Now the “tic” & “toc” command entered in the code has calculated the time taken to finish whole “FOR LOOP” of the “ga ()” syntax ,which we will use in study_2.

STUDY _ 2 : STATISTICAL BEHAVIOUR WITH UNBOUNDED INPUTS AND INCREASED POPULATION.

CODE :

% statistical Behaviour with unBounded Inputs

clear all 
close all
clc

% Defining our Search Space 1-D Array
x = linspace(0,0.6,150);
y = linspace(0,0.6,150);
num_case = 50;

% Creating 2-d Array
[xx,yy] = meshgrid(x,y)

% Calculating Fitting Function
for i = 1:length(xx);
    for j = 1:length(yy);
        input(1) = xx(i,j);
        input(2) = yy(i,j);
        f(i,j) = stalagmite_func(input);
    end
end

% Converting Fitting Function to change the graph upwards
F = -f;


% Performing the Genetic Algorithm
tic
for m =1:num_case
[inputs , fopt(m)] = ga(@stalagmite_func,2,[],[],[],[],[0;0],[0.6;0.6]);
xopt(m) = inputs(1);
yopt(m) = inputs(2);
end
Fopt = - fopt;
study2_time = toc

% Plotting Global Maxima And Function Values
figure(1)
subplot(2,1,1)
hold on
surface(x,y,F)
shading interp
plot3(xopt,yopt,Fopt,'marker','o','markersize',5,'markerfacecolor','r');
xlabel ('X-value');
ylabel ('y-value');
title('Unbounded Inputs');
legend('STALAGMITE','GLOBAL MAXIMA')
subplot(2,1,2)
hold on
plot(Fopt)
xlabel ('Iterations');
ylabel ('Function Maxima');
title('Variation in F');
legend('Change in F value')

EXPLANATION:

1. The above code is same that’s of previous code, only the change is in the syntax of “ga ()” command.
2. Now in this code the syntax used for “ga ()” command is

X = ga(fun,nvars,A,b,Aeq,beq,lb,ub) define a set of lower and upper bounds on the design variable,x,so solution is found between the range (Set A = [] and b = [] if no linear inequalities exist.) (Set Aeq = [] and beq = [] if no linear equalities exist)

3. The above syntax help us to enter upper and lower bound value as input to get the better solution than previous code. This means that when “ga()” calculates the value of ”f” ,it will display on the value of “f” that lies between the bounded values of “x” and “y”.
4. Now it is not compulsory to take the values of “x” and ”y” values only as your bounded values , these values are taken according to the problem we are working on , sometimes for the bounded values to be specified need some practical experience for some situations .So bounded values are taken according to our problem.
5. When we run the code we get the following figure,

6. In the above figure we can observe that now the graph is much better clear then previous one and now –ga() predict that the global maxima is somewhere between three points that are in first plot.
7. When we look at the second plot, we can observe that our second value is between 0.6 to 1
8. But here we still didn’t get the optimized value ,i.e results are yet not clear. To overcome this problem , further we are using modified “ga ()” command.

6. In above figure we can see that time taken by “ga ()” to calculate ,is more ,compared to previous study . that means me have entered more inputs “ga” is taking more taking more time to generate the solution.

STUDY_3 : STATITICAL BEHAVIOUR WITH BOUNDED INPUTS AND INCREASED POPULATION

CODE :-

% statistical Behaviour with Bounded Inputs

clear all 
close all
clc

% Defining our Search Space 1-D Array
x = linspace(0,0.6,150);
y = linspace(0,0.6,150);
num_case = 50;

% Creating 2-d Array
[xx,yy] = meshgrid(x,y)

% Calculating Fitting Function
for i = 1:length(xx);
    for j = 1:length(yy);
        input(1) = xx(i,j);
        input(2) = yy(i,j);
        f(i,j) = stalagmite_func(input);
    end
end

% Converting Fitting Function to change the graph upwards
F = -f;

% Study 3 - Increasing Iteration
options = optimoptions ('ga');
options = optimoptions(options,'PopulationSize',600);


% Performing the Genetic Algorithm
tic
for p =1:num_case
[inputs , fopt(p)] = ga(@stalagmite_func,2,[],[],[],[],[0;0],[0.1;0.1]);
xopt(p) = inputs(1);
yopt(p) = inputs(2);
end
Fopt = - fopt;
study3_time = toc

% Plotting Global Maxima And Function Values
figure(1)
subplot(2,1,1)
hold on
surface(x,y,F)
shading interp
plot3(xopt,yopt,Fopt,'o','markersize',5,'markerfacecolor','r');
xlabel ('X-value');
ylabel ('y-value');
title('Bounded Inputs and Increasing Population');
legend('STALAGMITE','GLOBAL MAXIMA')
subplot(2,1,2)
hold on
plot(Fopt)
xlabel ('Iterations');
ylabel ('Function Maxima');
title( 'F Maximum');
legend('Change in F value')

EXPLANATION:-

1. The code is same that of the previous two code , the only change is in the “ga ()” command and one extra command is added to code.
2. In this code after the fitting function is calculated , we have a command called “optimoptions()”, the syntax for this command is

Options = optimoptions(SolverName)returns a set of default options for the SolverName solver.

Options = optimoptions(SolverName,Name,Value) return options with name parameters altered with the specified values.

3. This particular code helps us to write the default values of any inbuilt MATLAB function . for the 1^st syntax when we run it we will get all the default value of the solver function which is mentioned in the bracket. Here in our code we have given “ga” in the bracket so this command will return default value of “ga” function in MATLAB.

3. In the above figure we can see what the 1^st syntax does ,Here we can see that the default value for “Population Size “ is 50 . That means, when “ga” takes random samples as inputs ,it takes those value from a population of 50 , so the number of option to choose for “ga” is less .
4. If we increase the value of “PopulationSize” then “ga” will have more samples tp choose and perform mutation and cross-over then it can give us more accurate optimized value.
5. To change this value we used 2^nd syntax of “optimoptions()” i.e. we have define which solver the value have to changed, then the parameter in that solver which has to be changed and then what value it has to be replaced.
6. Thus by doing this we have changed the “PopulationSize” to 300
7. Then, in the “ga()” function we have used the following syntax to give this change as one of the input.

X = ga(fun,nvars,A,b,Aeq,beq,lb,ub,nonlcon,options)minimizes with the default optimization parameter replaced by values in option .

9. Here in the above syntax we have enter all the inputs as that of previous code,additionally we can also give nonlinear constraints and “options” i.e. (change of default values as inputs in “ga ()” function.
10. Since here we don’t have any nonlinear constraints, we have entered [] in its place and we have entered the “options” as inputs to tell the “ga ()” that there is change in default values.
11. Then rest of the code is same as the previous code and when we run this code we get following results.

12. From the above fig we can observe that by increasing the population size we have got a single optimized value in the 1^st graph . If we increase the population size more further, the value will be more accurate .
13. In the 2^nd plot we can observe that maximum value of for many iteration is near to 0.9999985 and only in 3 iteration “ga ()” predicts that value of “f” are less that of other iteration values, So the maximum value is somewhere near to this value.
14. This problem can be further rectified by increasing the population size and our graph will be much more clear.

15. Now we see in the command window that the time taken by “ga ()” function is more than that for previous codes. So as we increase the inputs “ga ()” takes more time to calculate. This is one of the drawback of “ga ()” command.

WORKSPACE:

CONCLUSION:

• The idea behind the above studies is to observe how the “ga ()” commands works with different inputs and how to obtain a same value for a function even if “ga” is run several times.
• Thus Global Maxima is found for the stalagmite function.
• Genetic Algorithm is one of the most used method for obtained optimized value.
• Genetic Algorithm is mostly used to find the minimum value of a function , for example to find minimum drag force that will be created by the Aerofoil
• Similar approach as above code can be used to find maximum value using “ga” command

REFERENCE :

1. https://www.youtube.com/watch?v=8DZH20i52jw&t=2258s
2. https://www.youtube.com/watch?v=8DZH20i52jw&list=TLPQMTgwMTIwMjE80zYapPSC5g&index=1
3. https://www.youtube.com/watch?v=nas28_u-dJw&list=PLYwpaL_SFmcDHUTN26NXKfjg6wFJKDO9R&index=2
4. https://www.youtube.com/watch?v=skqn8nAwt5w
5. https://www.youtube.com/watch?v=HM6RhxWcXnY&t=57s
6. https://www.youtube.com/watch?v=LREPlWJmvAY