Implementation of Genetic Algorithm to determine Global Maxima of Stalagmite Function using MATLAB
AIM: Implementation of Genetic Algorithm to determine Global Maxima of Stalagmite Function using MATLAB.
OBJECTIVES:
- To understand Genetic Algorithm and Stalagmite Function.
- To identify and distinguish between different syntax available for Genetic Algorithm in MATLAB.
- To code a MATLAB programme for implementation of Genetic Algorithm to determine Global Maxima of Stalagmite function for (a) unbounded inputs (b) bounded inputs and (c) bounded inputs with increased population size.
GENETIC ALGORITHM:
Before starting with the introduction of Genetic Algorithm let’s see what is optimization? And how it is useful to engineering systems.
What is Optimization?
When it comes to engineering systems such as IC engines or air craft wings there are certain things which needs perform very well for example aerofoil needs to provide the proper lift at different angle of attack and at the same time needs to minimize drag. It means that for the different flight conditions the drag value should be minimum. When we say the drag value should be minimum, “how do we modify the designs to get there?” surely, we can modify all the design parameters and then see which design parameter combination gives lowest drag coefficient but this process is quite exhausting it’s going to take a lot of time to use a brute force method to calculate the best design point and hence we use optimization. Optimisation basically an area of study that deals with how to minimize or maximize a particular function which depends on multiple input parameters.
There are several ways to go for optimization today we are just going to introduce this area. Stalagmite function is a mathematical model that actually helps us to generate local minima and maxima. It is basically a 3D function. The way the function is designed, it has lots of local maxima, local maxima refers to tiny peaks the top of little mountain there is a local maxima. Point at the highest peak is global maxima. Now the objective of the optimization study is to determine that point, without the prior knowledge of this particular field. In this case we are going to start with the known function because if the function is known we know optimization points and we can actually march towards it. For this particular problem we are going to use algorithm called as genetic algorithm. Genetic algorithm is based upon evaluation theory that was proposed by Charles Darwin and MATLAB provides inbuilt functions that actually helps us to do Genetic Algorithm quite comfortably.
What is Genetic Algorithm?
Genetic Algorithms are the heuristic search and optimization techniques that mimic the process of natural evolution. The principle of natural selection is can be best described by the quote “Select The Best, Discard The Rest”
The process of natural selection starts with the selection of fittest individuals 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 will be better than parents and have a better chance at surviving. This process keeps on iterating and at the end, a generation with the fittest individuals will be found.
GA Operators and Parameters
- Initial Population: The process begins with a set of individuals which is called a Population. Each individual is a solution to the problem you want to solve.
- Fitness Function: The fitness function determines how fit an individual is (the ability of an individual to compete with other individuals). It gives a fitness score to each individual. The probability that an individual will be selected for reproduction is based on its fitness score.
- Selection: The primary objective of the selection operator is to emphasize the good solutions and eliminate the bad solutions in a population while keeping the population size constant.“Selects the best, discards the rest”
- 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.
- Mutation: In certain new offspring formed, some of their genes can be subjected to a mutation with a low random probability. This implies that some of the bits in the bit string can be flipped.
- Termination: The algorithm terminates if the population has converged (does not produce offspring which are significantly different from the previous generation). Then it is said that the genetic algorithm has provided a set of solutions to our problem.
Syntax of Genetic Algorithm in MATLAB & its description
In MATLB, the function for a genetic algorithm is 'ga'. It always tries to minimise the given fitness function. Therefore, Genetic Algorithm can be used when a certain function needs to be minimised. For example, minimising some error.
The syntax for the 'ga' is as:
[x, fval]= ga (fun, nvars, A, b, Aeq, beq, lb, ub, nonlcon, IntCon, options)
Where:
'x' returns the input values.
'fval' returns the value of the fitness function.
'fun' is the objective function which is required to be minimise.
'nvars' are the number of variables associated with the function.
'A' and 'b' are the linear inequality constraint of the type A*x ≤ b, if any.
'Aeq' and 'beq' are the linear equality constraint of the type Aeq*x = beq, if any.
lb and ub set the lower and upper bound for the GA. These provide the GA the search-space. Otherwise, the GA will try to look at points beyond the function's domain.
'nonlcon' sets the non-linear constraint, if any.
'IntCon' sets the integer constraint, if any.
'options' are used whenever we need to change the way of the MATLABs default GA functioning to the desired optimisation requirements. For example, MaxGenerations is used to change the default number of generations that the GA uses before stopping the search. (Generally 100*number of variables). Similarly, the 'PopulationSize' is used when we need to change the size of the population to be considered (default size = 50 for 5 or less no of input variables).
Note: When there are no constraints so to change the upper and lower bounds, we need to provide [] (set of empty square brackets) for each entry that comes before the ‘ub’ and ‘lb’ in the syntax of GA function. For example in the given problem to change the upper and lower bounds only. The syntax becomes:
[inputs, fopt(i)] = ga(@stalagmite, 2, [], [], [], [], [0;0], [1;1])
Stalagmite Function:
- It is a mathematical model that helps to generate local minima and maxima.
- It is a three-dimensional function and the way the function is designed is that it has a lot of local maxima or minima and one global maxima or minima.
- The name stalagmite is actually a type of rock formation that rises from the floor of a cave due to the accumulation of material deposited on the floor from ceiling drippings.
PROGRAMME CODE IN MATLAB:
Software used: MATLAB R2020a trial version
Programme for Stalagmite function: I have modified the stalagmite function by adding the negative sign. Therefore, when the MATLAB performs GA it calculates the miminum value of negative product, which in turn is the point of maximum value of stalagmite function.
function [F] = stalagmite(input_vector) % Function Defination that creates a stalagmite function with base equations
% [F] = Output
% input_vector = input
x=input_vector(1);
y=input_vector(2);
% The equation for a stalagmite function depends upon two input variables x and y
% The output 'F' is calculated with the following formulaes.
f1_x = (sin(5.1*pi*x+0.5))^6;
f1_y = (sin(5.1*pi*y+0.5))^6;
f2_x = exp((-4*log(2))*((x-0.0667)^2)/0.64);
f2_y = exp((-4*log(2))*((y-0.0667)^2)/0.64);
[F] = -(f1_x.*f1_y.*f2_x.*f2_y);
% F will give a Stalgmite function which is an upward mounting peaks function
% multiplying the product with (-1) will invert the values towards the negative side.
% This indeed is the Stalagmite function with peaks (troughs) facing downwards.
% This is required in order to find the minimum of the function as the Genetic Algorithm function 'ga' of MATLAB only gives us the minimum values.
endProgramme for Genetic Algorithm
%Programme for Implementation of Genetic Algorithm to determine Global Maxima of Stalagmite Function
clear all % clears all workspace variables
close all % closes all plot and figures
clc % clears the command window screen
% Defining our search space % 1-Dimensional input arry, this does not define search space for GA
% To set that, we need to use the syntax inputs of the ga function.
x = linspace(0,0.6,150); % linespace for the X-axis, i.e., x values ranges from 0 to 0.6 in 150 steps
y = linspace(0,0.6,150); % linespace for the Y-axis, i.e., y values ranges from 0 to 0.6 in 150 steps
% Calculating a 2Dimentionl mesh
[xx,yy]= meshgrid(x,y); % created 2-D array based on 1-D array.
%Calculation of Stalagmite function
for i= 1:length(xx)
for j= 1:length(yy)
input_vector(1)= xx(i,j); % for each value of xx creating an array called as input_vector
input_vector(2)= yy(i,j); % Contents of input_vector are x value and y value which are put together in an array
f(i,j)=stalagmite(input_vector); % stalagmite takes the input_vector and calculates F
end % two for loops are basically calculating the data for plot
end
%Total Number of Cases or Iterations
num_cases = 50; % If we use ga only one time or num_cases = 1, we may not get the desired value that we look for.
% Each iteration may give different results from the ga and thus increasing the number of
% iterations will increase the changes ofthe occurrence (converge to the actual) of the best solution.
% STUDY 1: STATISTICAL BEHAVIOUR WITH UNBOUNDED INPUTS
tic
for i=1:num_cases
[inputs,fopt(i)]=ga(@stalagmite,2); % function ga is inbuilt in the matlab & we are passing two inputs
xopt(i)=inputs(1); % @stalagmite is the function handle
yopt(i)=inputs(2); % at 'inputs' ga will give maximum value 'fopt'
end % ga(@stalagmite, 2) in this 2 stands for number of varaiables function have
study_1_time=toc % to measure how long iteration is going to take to disply result
figure(1)
subplot(2,1,1) % to show multiple plot in same figure window
hold on
surfc(xx,yy,-f) % first plot to be created in figure (1)
shading interp
plot3(xopt,yopt,-fopt,'marker','o','markersize',5,'markerfacecolor','r') % plot3 command to plot xopt,yopt & fopt values that GA predicts for 50 studies
% this helps us to visualize how alogorithm is performing
xlabel('X values')
ylabel('Y values')
zlabel('Z/ function values')
title('STATISTICAL BEHAVIOUR WITH UNBOUNDED INPUTS')
subplot(2,1,2) % second plot to be created in figure (2)
plot(-fopt) % to make above plot clear we plot optimum value that the function is predicting as a function of iteration
xlabel('Iterations')
ylabel('Function Minimum')
title('Minimum Value At Different Iterations')
% STUDY 2: STATISTICAL BEHAVIOUR WITH UPPER & LOWER BOUNDS
for i = 1:num_cases
[inputs , fopt(i)] = ga(@stalagmite, 2,[],[],[],[],[0;0],[1;1]);
% 'inputs' stores the x and y values that the ga algorithm returns
% fopt stores the optimum value that each iteration of the ga returns. This stores the minimum z (function) values
% The [] square brakets are necessary to provide in the function handle even if they are empty.
% They are for linear equality and non equality constrains.
% The [0;0] and [1;1] are the lower and upper bounds of the function.
% If this is not properly set, values beyond the function domain will be returned by ga.
% In other words, these provide the linespace for the ga function.
xopt(i) = inputs(1); % stores the values of x returned by ga at each iteration
yopt(i) = inputs(2); % stores the values of y returned by ga at each iteration
end
figure(2)
subplot(2,1,1)
hold on
surfc(xx,yy,-f)
shading interp
plot3(xopt,yopt,-fopt,'marker','o','markerfacecolor','r')
% plot3 plots a line in 3D-space through the points xopt, yopt and fopt
% plotting the point on the same figure with xopt, yopt and zopt coordinates
% with a red marker i.e., the minimum value returned by the ga optimisation.
xlabel('X values')
ylabel('Y values')
zlabel('Z/ function values')
title('STATISTICAL BEHAVIOUR WITH UPPER & LOWER BOUNDS')
legend('Function','Minima') %provides visualisation on the figure
subplot(2,1,2)
plot(-fopt)
xlabel('Iterations')
ylabel('Function Minimum')
title('Minimum Value At Different Iterations')
% STUDY 3: INCREASING GA ITERATIONS
% In addition to specifing the bounds, we are also going to change initial population size
% as we know GA starts with random sample so by defualt GA takes 50 as a number of sample
% since there are two varaible we are going to have 50 sets of two random
% varaibles i.e. going to be population, if initial population is too
% small, it can actully not find the global maxima but it can actully find
% local maxima which is what happned study 2
% increasing population size to 300, it is going to slow down GA becouse of
% more number of samples, more indivuduals and they need to be evalauted
options = optimoptions('ga');
options = optimoptions(options,'PopulationSize',300);
tic
for i=1:num_cases
[inputs, fopt(i)]= ga(@stalagmite,2,[],[],[],[],[0;0],[1;1],[],[],options);
xopt(i)=inputs(1);
yopt(i)=inputs(2);
end
study3_time= toc;
figure(3)
subplot(2,1,1)
hold on
surfc(xx,yy,-f)
shading interp
plot3(xopt,yopt,-fopt,'marker','o','markersize',5,'markerfacecolor','r')
title('BOUNDED INPUTS WITH INCREASED POPULATION SIZE')
xlabel('X values')
ylabel('Y values')
zlabel('Z/ function values')
subplot(2,1,2)
plot(-fopt)
xlabel('Iterations')
ylabel('Function Minimum')
title('Minimum Value At Different Iterations')
% To find the maximum value of function and the point at which it occurs
% To extract a single minimum value from the fopt and to convert it into the maximum value
maxvalue = -(min(fopt));
display('The max value of Stalagmite function is :')
display(maxvalue)
display('And the point at which this occurs is :')
display(inputs)
N = ['The study took ',num2str(study3_time),' Seconds']
disp(N)
Study 1 Output: STATISTICAL BEHAVIOUR WITH UNBOUNDED INPUTS
In this study we perform the Genetic Algorithm without setting any bounds and hence we get absurd results for x and y ouside the doman as well as many values because of which any single value cannot be determined.
Output of Study 2: STATISTICAL BEHAVIOUR WITH UPPER & LOWER BOUNDS
In this study, we set the lower and upper bounds of x and y as shown in the code and syntax. As we set the bounds, the most likely answer is coming near 3 peaks as shown in figure and the minimum value ranges between 0.5 and 1. To get more refined answer we will increase the population size in the next study.
Output of Study 3: BOUNDED INPUTS WITH INCREASED POPULATION SIZE
The default population size in GA fuction is considers as 50, we are getting absurd results as in above study. To get accurate and considerable results increase the population size of GA using optimoptions so that more generations and offsprings are produced that will lead the result towards convergence. I set the population size to 300. Now, in the result the minimum value occurs between 0.9999997 and 0.99999985 which approximates to 1 and the point of minima is at the highest peak,i.e. the global minima and that occurs at (0.0668, 0.0668).
In my study, the global maximum is at (x,y) = (0.0668, 0.0668). The value is 1 as shown in figure below
REFRENCES: