Optimization of Stalagmite function using Genetic Algorithm in MATLAB
OBJECTIVE:
To optimize the Stalagmite Function and find the global maxima of the function using the Genetic Algorithm in MATLAB.
THEORY:
- GLOBAL OPTIMIZATION:
Global optimization is a branch of applied mathematics and numerical analysis that attempts to find the global minima or maxima of a function or a set of functions on a given set. A few application examples include acoustics equipment design, cancer therapy planning, chemical process modeling, data analysis, classification and visualization, economic and financial forecasting, environmental risk assessment and management, industrial product design, laser equipment design, model fitting to data (calibration), optimization in numerical mathematics, optimal operation of "closed" (confidential) engineering or other systems, packing and other object arrangement problems, portfolio management, potential energy models in computational physics and chemistry, process control, robot design and manipulations, systems of nonlinear equations and inequalities, and waste water treatment systems management.
- GENETIC ALGORITHM:
A genetic algorithm is a discover for something by themselves which was inspired by Charles Darwin’s theory of natural evolution.This algorithm reflects the process of natural selection where the fittest individuals are selected for reproduction in order to produce offspring of the next generation.
Five phases are considered in a genetic algorithm.
- 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 idea of selection phase is to select the fittest individuals and let them pass their genes to the next generation.
- 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. Mutation occurs to maintain diversity within the population and prevent premature convergence.
- ga FUNCTION IN MATLAB:
Basically, what a ga function in matlab does is calculates the minimum of the function by default.
The general syntax for the ga function is
x = ga(fun,nvars,A,b,[],[],lb,ub,nonlcon,IntCon,options)
where
x= minimum value of resultant function
fun= fitness function
nvars= no of design variables
A= Linear inequality constraints, specified as real matrix A.
b= Linear inequality constraints, specified as real vector b.
[]= Linear equality constraints, specified as real matrix and real vector.
lb= lower bound in the search space of the algorithm
ub= upper bound in the search space of the algorithm
nonlcon= non linear constraints which specified as function name.
Int con= Integer variables, specified as a vector of positive integers taking values from 1 to nvars.
options= Create options by using optimoptions or by exporting options from the Optimization app. For details, see Importing and Exporting Your Work (Optimization Toolbox).
- GLOBAL MAXIMA AND MINIMA:
Global maxima or minima of f(x) in [a, b] is basically the greatest or least value of f(x) in [a, b]. Mathematically, it is written as:
The function f(x) has a global maximum at the point ‘a’ in the interval I if f (a) ≥f(x), for all x∈I.
Similarly, f(x) has a global minimum at the point ‘a’ if f (a) ≤f (x), for all x∈I.
Global maxima or minima in [a, b] will always occur either at the critical points of f(x) within [a, b] or at the end points of the interval.
- STALAGMITE FUNCTION
Stalagmite represents the 3D pointed formations like a cone formed on ploting the given data values. It basically, creates local maximas and minimas for the given data. Thus, this helps us finding design characteristics of a product like drag coefficient of the aircraft wing or a car body.
PROGRAM:
- EXPLANATION OF THE CODE:
EXPLANATION OF BASIC CODE:
clear all
close all
clc
%Defining our search space
x= linspace(0,0.6,150);
y= linspace(0,0.6,150);
% Creating 2D mesh
[xx yy]= meshgrid(x,y);
%Evaluating the stalagamite 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
figure(11)
surfc(xx,yy,f)
xlabel('x value')
ylabel('y value')
shading interp
[inputs, fval]= ga(@stalagmite, 2);1) First, we provide input data of the program i.e. the x values and y values using the linspace command.
2) Then, we perform the meshgrid operation to form the matrix of xx and yy or rather creating a 2D mesh.
3) We then evaluate the stalagmite function using nested for loops for the 2D mesh data formed.
EXPLANATION OF FUNCTION CODE:
function [f]= stalagmite(input_vector);
a= input_vector(1);
b= input_vector(2);
f1_x= [sin(5.1*pi*a+0.5)]^6;
f1_y= [sin(5.1*pi*b+0.5)]^6;
f2_x= exp(-4*log(2)*(((a-0.0667)^2)/(0.64)));
f2_y= exp(-4*log(2)*(((b-0.0667)^2)/(0.64)));
f = -(f1_x*f1_y*f2_x*f2_y);
end4) In Stalagmite function, we provide the input vectors which are basically the 2D mesh xx and yy formed before we calculate f1x,f2x,f1y and f2y values.
5) After we calculate the values of f1x,f2x,f1y and f2y we find the resultant function value f.
6) We are finding global maxima hence, we invert the function to the global maximum value.
7) After we find resultant we plot the graph of 3D data using surfc command and we also give shading interp to give a smooth shading on the surface of the stalagmite.
EXPLANATION OF ALL THE STUDIES TO FIND THE GLOBAL MAXIMA:
%Study1- Statistical Behaviour
num_case=50
tic
for i= 1:num_case
[inputs,fopt(i)]= ga(@stalagmite, 2);
xopt(i)= inputs(1);
yopt(i)= inputs(2);
end
study_time_1= toc
figure(1)
subplot(2,1,1)
hold on
surfc(x,y,f)
shading interp
plot3(xopt, yopt, fopt,'marker','o','markersize',5,'markerfacecolor','g')
title('unbounded inputs')
subplot(2,1,2)
plot(fopt)
xlabel('iterations')
ylabel('function maximum')
%Study2- Statistical Behaviour- with upper and lower bounds
tic
for i= 1:num_case
[inputs,fopt(i)]= ga(@stalagmite, 2,[],[],[],[],[0;0],[1;1]);
xopt(i)= inputs(1);
yopt(i)= inputs(2);
end
study_time_2=toc
figure(2)
subplot(2,1,1)
hold on
surfc(x,y,f)
shading interp
plot3(xopt, yopt, fopt,'marker','o','markersize',5,'markerfacecolor','g')
title('bounded inputs')
subplot(2,1,2)
plot(fopt)
xlabel('iterations')
ylabel('function maximum')
%Study3- Statistical Behaviour- increasing the ga iterations
tic
options= optimoptions('ga');
options= optimoptions(options,'populationsize',350);
for i= 1:num_case
[inputs,fopt(i)]= ga(@stalagmite, 2,[],[],[],[],[0;0],[1;1],[],[],options);
xopt(i)= inputs(1);
yopt(i)= inputs(2);
end
study_time_3= toc
figure(3)
subplot(2,1,1)
hold on
surfc(x,y,f)
shading interp
plot3(xopt, yopt, fopt,'marker','o','markersize',5,'markerfacecolor','g')
title('bounded inputs')
subplot(2,1,2)
plot(fopt)
xlabel('iterations')
ylabel('function maximum')
8) Here we perform the study 1 for unbounded inputs meaning we don't provide axes or limits for inputs, wherein we are providing num_case= 50 which the no of iterations the function should perform to find the global maxima. The default ga iterations are set to 50.
9) Here tic and toc commands are used to find the time taken to give global maximas. here we use the same for loop to find stalagmite function. Here we use the subplot command to plot two graphs or plots in the same figure window subplot(m,n,p) where m is the no of rows, n is the no of columns and p is the position of the subplot.
10) In study 2, we take bounded inputs which basically is the constraints or boundry conditions given by us.
11) In study 3, we take bounded inputs and control the population size to obtain a global maxima.The optimoptions command returns options with the named parameters altered with the specified values.
- FULL CODE:
clear all
close all
clc
%Defining our search space
x= linspace(0,0.6,150);
y= linspace(0,0.6,150);
% Creating 2D mesh
[xx yy]= meshgrid(x,y);
%Evaluating the stalagamite 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
figure(11)
surfc(xx,yy,f)
xlabel('x value')
ylabel('y value')
shading interp
[inputs, fval]= ga(@stalagmite, 2);
%Study1- Statistical Behaviour
num_case=50
tic
for i= 1:num_case
[inputs,fopt(i)]= ga(@stalagmite, 2);
xopt(i)= inputs(1);
yopt(i)= inputs(2);
end
study_time_1= toc
figure(1)
subplot(2,1,1)
hold on
surfc(x,y,f)
shading interp
plot3(xopt, yopt, fopt,'marker','o','markersize',5,'markerfacecolor','g')
title('unbounded inputs')
subplot(2,1,2)
plot(fopt)
xlabel('iterations')
ylabel('function maximum')
%Study2- Statistical Behaviour- with upper and lower bounds
tic
for i= 1:num_case
[inputs,fopt(i)]= ga(@stalagmite, 2,[],[],[],[],[0;0],[1;1]);
xopt(i)= inputs(1);
yopt(i)= inputs(2);
end
study_time_2=toc
figure(2)
subplot(2,1,1)
hold on
surfc(x,y,f)
shading interp
plot3(xopt, yopt, fopt,'marker','o','markersize',5,'markerfacecolor','g')
title('bounded inputs')
subplot(2,1,2)
plot(fopt)
xlabel('iterations')
ylabel('function maximum')
%Study3- Statistical Behaviour- increasing the ga iterations
tic
options= optimoptions('ga');
options= optimoptions(options,'populationsize',350);
for i= 1:num_case
[inputs,fopt(i)]= ga(@stalagmite, 2,[],[],[],[],[0;0],[1;1],[],[],options);
xopt(i)= inputs(1);
yopt(i)= inputs(2);
end
study_time_3= toc
figure(3)
subplot(2,1,1)
hold on
surfc(x,y,f)
shading interp
plot3(xopt, yopt, fopt,'marker','o','markersize',5,'markerfacecolor','g')
title('bounded inputs')
subplot(2,1,2)
plot(fopt)
xlabel('iterations')
ylabel('function maximum')
- ANSWERS TO QUESTIONS:
1) Why the num case is considered as 50?
The num case is taken as 50 because we are running the ga for 50 times and we can change the values also for the num case.
2) Why do we provide Upper and Lower bounds?
The upper and lower bounds are provided to set a constraint or provide boundry conditions from which initial value to which final value should the inputs be taken to find the global maxima.
3) Why do we increase the population size?
We increase the population size to find the a single global maxima point which is our needed results instead of find a bunch of local maxima which got in study 2. While, increasing the population size.
4) What are Stopping Criteria of Genetic Algorithm?
The genetic algorithm can be stopped using the following conditions
1. Generations
2. Time limit
3. Fitness limit
4. Stall generations
5. Stall time limit
6. Stall test
7. Function tolerance
8. Constraint tolerance
The genetic algorithm stops if any one of these conditions are met.
- ERRORS WHILE PERFORMING THE CODE AND CORRECTION:
1) I had found that the time for each study as wasn't printed at the end of every study. Since, I had put a semicolon(;) for study time the output wasn't displayed. I had corrected it and ran the program again I got the study time.
2) I had got another error too while running program i.e. the value of Aeq and beq wasn't inputted.
This was corrected by providing the inputs for Aeq and beq as [],[].
OUTPUT:
- PLOT FOR GENERAL STALAGMITE FUNCTION:
- PLOT FOR STALAGMITE WITH UNBOUNDED INPUTS:
The study time for this is also computed.
- PLOT FOR STALAGMITE WITH BOUNDED INPUTS:
The study time for this is also computed.
- PLOT FOR STALAGMITE WITH BOUNDED INPUTS AND INCREASED POPULATION SIZE:
The study time for this is also computed.
Hence, here we are getting a single global maxima value as -1 when the population size was set to 350.
REFERENCES:
2) https://www.askiitians.com/iit-jee-applications-of-derivatives/global-maxima-and-minima/
3) https://in.mathworks.com/help/gads/ga.html?s_tid=srchtitle
4) https://www.mathworks.com/help/releases/R2019b/gads/how-the-genetic-algorithm-works.html#f6243