Calculation of Global Maxima of a Stalagmite Function using Genetic Algorithm

Concept of Genetic Algorithm :-

Genetic algorithm is the method of finding optimum solutions for the problems based on the theory of evolution by Charles Darwin. It reflects the process of natural selection of the fittest individuals (parents) from the current generation (population) to produce an individual (child) of the next generation (population).

Phases in genetic algorithm are :-

1) Generation of initial population - Initially, individuals are randomly selected by algorithm and current population is generated which become parents for the next generation. In genetic algorithm, the set of genes of an individual is represented by a string. Usually, string is a set of binary values (1\'s and 0\'s).

2) Evaluation using Fitness Function - Fitness of every individual of the current population is evaluated using the fitness function f(x, y).

3) New Population - The process 1 and 2 are repeated until the required no of individulas are selected and the current population is generated.

4) Selection - Individuals from the current population are selected on the basis of their fitness value.

5) Crossover - It is the most significant phase in genetic algorithm. For each pair of parents, a crossover point is chosen at random. Everything before this point is copied from first parent and then everything after a crossover point from the second parent.

6) Mutation - Mutation is performed to prevent falling all solutions in population into a local optimum of solved problem. In mutation, some bits of the string are flipped (0 change to 1 and vice versa).

7) Termination - When two conditions are satisfied either fixed number of generations reached or the highest ranking solution\'s fitness is reaching or such that successive iterations no longer produce better results or both.

Genetic Algorithm (GA) in Matlab :-

Genetic algorithm is run in matlab using solver name \'ga\'. In Matlab, genetic algorithm gives minimum value of the function, but it can be used to calulate both minimum/maximum value of the function.

Syntax for ga - ga(function, nvars, A, b, Aeq, beq, lb, ub, nonlcon, IntCon, options)

ga - solver inbuilt in matlab

function - objective function (i.e function to be minimized or maximized)

nvars - number of design

A and b are matrix which represents linear inequality constraint such that A*x `<=` b (if linear inequality constraints are not mentioned in the problem then defined as A=[] and b=[])

Aeq and beq are matrix which represents equality constraint such that Aeq*x = b (if equality constraints are not mentioned in the problem then defined as Aeq=[] and beq=[])

lb - lower bound/limit for the design variables (in matrix form)

ub - upper bound/limit for the design variables (in matrix form)

nonlcon - sets nonlinear constraints, if any

IntCon - sets integer constraint, if any

options - these are the attributes with which we can change the default functioning of ga, for example population size, which if not mentioned, ga takes random value.

Code to calculate global maxima of stalagmite function using genetic algorithm:-

%PART I (MAIN PROGRAM)
%Program to calculate the global maxima using genetic algorithm
clear all
close all
clc

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

%Creating a 2D mesh
[X Y] = meshgrid(x, y);

%Evaluation of stalagmite function
for i = 1:length(X)
    for j = 1:length(Y)
        input_vector(1) = X(i, j);
        input_vector(2) = Y(i, j);
        F(i, j) = stalagmite(input_vector);
    end
end

%Number of iterations
num_cases = 100;

%Optimization of Stalagmite Function for Global Maxima 
% Study 1 - Statistical Behaviour 
tic 
for i = 1:num_cases
    [inputs1, fopt1(i)] = ga(@stalagmite, 2);
    xopt1(i) = inputs1(1);
    yopt1(i) = inputs1(2);
end
study1_time = toc

figure(1)
subplot(2, 1, 1)
hold on
surfc(X, Y, (-F))
xlabel(\'X-value\')
ylabel(\'Y-value\')
shading interp
plot3(xopt1, yopt1, (-fopt1), \'marker\', \'o\', \'markersize\', 5, \'markerfacecolor\', \'r\')
title(\'Study 1 : Unbounded Inputs\')
subplot(2, 1, 2)
plot((-fopt1))
xlabel(\'Iterations\')
ylabel(\'Function Maximum\')

% Study 2 - Statistical Behaviour with Bounded Inputs
tic 
for i = 1:num_cases
    [inputs2, fopt2(i)] = ga(@stalagmite, 2, [], [], [], [], [0; 0], [0.6; 0.6]);
    xopt2(i) = inputs2(1);
    yopt2(i) = inputs2(2);
end
study2_time = toc

figure(2)
subplot(2, 1, 1)
hold on
surfc(X, Y, (-F))
xlabel(\'X-value\')
ylabel(\'Y-value\')
shading interp
plot3(xopt2, yopt2, (-fopt2), \'marker\', \'o\', \'markersize\', 5, \'markerfacecolor\', \'r\')
title(\'Study 2 : Bounded Inputs\')
subplot(2, 1, 2)
plot((-fopt2))
xlabel(\'Iterations\')
ylabel(\'Function Maximum\')

%Study 3 - Statistical Behaviour with Bounded Inputs and increased
%Pouplation Size
options = optimoptions(\'ga\');
options = optimoptions(options, \'PopulationSize\', 200);

tic 
for i = 1:num_cases
    [inputs3, fopt3(i)] = ga(@stalagmite, 2, [], [], [], [], [0; 0], [0.6; 0.6], [], [], options);
    xopt3(i) = inputs3(1);
    yopt3(i) = inputs3(2);
end
study3_time = toc

figure(3)
subplot(2, 1, 1)
hold on
surfc(X, Y, (-F))
xlabel(\'X-value\')
ylabel(\'Y-value\')
shading interp
plot3(xopt3, yopt3, (-fopt3), \'marker\', \'o\', \'markersize\', 5, \'markerfacecolor\', \'r\')
title(\'Study 3 : Bounded Inputs with increased Population Size\')
subplot(2, 1, 2)
plot((-fopt3))
xlabel(\'Iterations\')
ylabel(\'Function Maximum\')

%Analytical Value of global maxima to know the real maximum value of stalagmite function (extra work)
max_value = max(abs(fopt3));
maxima = sprintf(\'Value of Global Maxima for Case 3 = %e\', max_value)


%PART II (SUBPROGRAM DEFINING STALAGMITE FUNCTION)
function [F] = stalagmite(input_vector)

f1x = ((sin((5.1*pi.*input_vector(1)) + 0.5)).^6);
f2x = (exp(((-4)*log(2).*(input_vector(1) - 0.0667).^2)./0.64));
f1y = ((sin((5.1*pi.*input_vector(2)) + 0.5)).^6);
f2y = (exp(((-4)*log(2).*(input_vector(2) - 0.0667).^2)./0.64));

F = -(f1x.*f2x.*f1y.*f2y)
end

Stepwise Explanation of Code :-

PART I - MAIN PROGRAM

1) First Step - Use of basic commands

Description-

clear all %This command clears the workspce

close all %This command closes all the plots and figures that are opened earlier

clc %This command clears the workspace

Output of the step - The workspace and command window gets cleared. All the plots get closed.

2) Second Step - Defining search space and creating a 2D mesh

Description -

x = linspace(0, 0.6, 150) %Row vector defined using linspace command to define search space

y = linspace(0, 0.6, 150) %Row vector defined using linspace command to define search space

[X Y] = meshgrid (x, y) %A mesh created using meshgrid command which is given as a input to the stalagmite function

Output of the step - Search space gets defined using linspace command and a mesh gets created using meshgrid command and passing x and y vectors as input arguments.

3) Third Step - Evaluation of stalagmite function

Description -

for i = 1:length(X) %First for loop defined with counter variable \'i\' varying from 1 to length(X) which is equal to 150 as defined earlier

for i = 1:length(Y) %First for loop defined with counter variable \'j\' varying from 1 to length(X) which is equal to 150 as defined earlier

input_vector(1) = X(i, j) %Denotes elements of 1st column of matrix name \'input_vector\' for every particular iteration (i.e for i and j values)

input_vector(2) = Y(i, j) %Denotes elements of 2nd column of matrix name \'input_vector\' for every particular iteration (i.e for i and j values)

F(i, j) = stalagmite(input_vector) %Value of stalagmite function for all elements of matrix \'input_vector\' (i.e for all values of x and y)

end %Command to end the second for loop

end %Command to end the first for loop

Output of the step - A matrix name as \'input_vector\' gets defined and stalagmite function gets evaluated for elements of matrix \'input_vector\' (i.e for all values of x and y).

4) Fourth Step - Evaluation and plotting of statistical behaviour of stalagmite function without any bounds and without defining population size

Description -

tic %Command for initiating timer

for i = 1:num_cases %Loop defined with counter variable \'i\' varying from 1 to num_cases which is equal to 100 to obtain the points at which the value of stalagmite function is maximum amd the maximum value

inputs1 %Variable (output argument) used to store matrix containing the values of x and y at which value of stalagmite function is maximum for all 100 iterations

fopt1 %Variable (output argument) used to store a row vector containing the maximum values of stalagmite function for all 100 iterations

ga(@stalagmite, 2) %ga solver used to perform genetic algorithm by defining input arguments like function name and no of design variables

xopt1(i) = inputs1(1) %Variable used to store the optimum values of x (i.e elements of 1st column of the matrix name \'inputs1\' for all 100 iterations

yopt1(i) = inputs1(2) %Variable used to store the optimum values of y (i.e elements of 2nd column of the matrix name \'inputs1\' for all 100 iterations

end %Command to end the for loop

study1_timer = toc %Variable used to store value time required to perform the genetic algorithm run using \'toc\' command (to stop the timer also)

subplot %Command used to create two stacked plots in one figure

hold on %Command used to include all other plots in one figure

surfc(X, Y, (-F)) %Command used to create a 3D contour plot (In this case, the vertices of the surface faces are (X(i), Y(i), F(i,j)) triples. To create X and Y matrices for arbitrary domains, use the meshgrid function.)

shading interp %Command used to obtain smoother plots without grid lines

plot3 %Plots one or more lines in three-dimensional space through the points whose coordinates are the elements of xopt1, yopt1 and -fopt1 (to get inverted plot)

title %Command used to give title to the plot

plot((-fopt1)) %To plot the (-fopt1) (to get inverted plot) vs iterations graph

xlabel %Command used to name the x-axis

ylabel %Command used to name the y-axis

Output of the step - Matrices inputs1 and fopt1 gets created after passing function and no of design variables as input arguments. Also optimum values of x and y (i.e xopt1 and yopt1) are obtained in the matrix form with 100 columns (1 element = 1 iteration). Time required to perform ga gets calulated. Plot of 3D contour and optimum value of function vs no of iterations are obtained in a single figure.

5) Fifth Step - Evaluation and plotting of statistical behaviour of stalagmite function by defining bounds and without defining population size

Description -

tic %Command for initiating timer

for i = 1:num_cases %Loop defined with counter variable \'i\' varying from 1 to num_cases which is equal to 100 to obtain the points at which the value of stalagmite function is maximum amd the maximum value

inputs2 %Variable (output argument) used to store matrix containing the values of x and y at which value of stalagmite function is maximum for all 100 iterations

fopt2 %Variable (output argument) used to store a row vector containing the maximum values of stalagmite function for all 100 iterations

ga(@stalagmite, 2, [], [], [], [], [0; 0], [0.6; 0.6]) %ga solver used to perform genetic algorithm by defining input arguments like function name, no of design variables, lower and upper bound

xopt2(i) = inputs2(1) %Variable used to store the optimum values of x (i.e elements of 1st column of the matrix name \'inputs2\' for all 100 iterations

yopt2(i) = inputs2(2) %Variable used to store the optimum values of y (i.e elements of 2nd column of the matrix name \'inputs2\' for all 100 iterations

end %Command to end the for loop

study2_timer = toc %Variable used to store value time required to perform the genetic algorithm run using \'toc\' command (to stop the timer also)

subplot %Command used to create two stacked plots in one figure

hold on %Command used to include all other plots in one figure

surfc(X, Y,(- F) ) %Command used to create a 3D contour plot (In this case, the vertices of the surface faces are (X(i), Y(i), F(i,j)) triples. To create X and Y matrices for arbitrary domains, use the meshgrid function.)

shading interp %Command used to obtain smoother plots without grid lines

plot3 %Plots one or more lines in three-dimensional space through the points whose coordinates are the elements of xopt2, yopt2 and -fopt2 (to get inverted plot)

title %Command used to give title to the plot

plot((-fopt2) %To plot the (-fopt2)(to get inverted plot) vs iterations graph

xlabel %Command used to name the x-axis

ylabel %Command used to name the y-axis

Output of the step - Matrices inputs2 and fopt2 gets created after passing function, no of design variables and bound limits as input arguments. Also optimum values of x and y (i.e xopt2 and yopt2) are obtained in the matrix form with 100 columns (1 element = 1 iteration). Time required to perform ga gets calulated. Plot of 3D contour and optimum value of function vs no of iterations are obtained in a single figure. Plots are more acurate as compared to earlier one.

6) Sixth Step - Evaluation and plotting of statistical behaviour of stalagmite function by defining bounds and population size

Description -

options = optimoptions(\'ga\') %This command returns the set of default options for the solver name

options = optimoptions(options, \'PopulationSize\', 200) %This command returns a copy of old options with the named parameters (i.e population size) altered with the specified values (i.e 200).

tic %Command for initiating timer

for i = 1:num_cases %Loop defined with counter variable \'i\' varying from 1 to num_cases which is equal to 100 to obtain the points at which the value of stalagmite function is maximum amd the maximum value

inputs3 %Variable (output argument) used to store matrix containing the values of x and y at which value of stalagmite function is maximum for all 100 iterations

fopt3 %Variable (output argument) used to store a row vector containing the maximum values of stalagmite function for all 100 iterations

ga(@stalagmite, 2, [], [], [], [], [0; 0], [0.6; 0.6], options) %ga solver used to perform genetic algorithm by defining input arguments like function name, no of design variables, lower and upper bound

xopt3(i) = inputs3(1) %Variable used to store the optimum values of x (i.e elements of 1st column of the matrix name \'inputs2\' for all 100 iterations

yopt3(i) = inputs3(2) %Variable used to store the optimum values of y (i.e elements of 2nd column of the matrix name \'inputs2\' for all 100 iterations

end %Command to end the for loop

study3_timer = toc %Variable used to store value time required to perform the genetic algorithm run using \'toc\' command (to stop the timer also)

subplot %Command used to create two stacked plots in one figure

hold on %Command used to include all other plots in one figure

surfc(X, Y, (-F)) %Command used to create a 3D contour plot (In this case, the vertices of the surface faces are (X(i), Y(i), F(i,j)) triples. To create X and Y matrices for arbitrary domains, use the meshgrid function.)

shading interp %Command used to obtain smoother plots without grid lines

plot3 %Plots one or more lines in three-dimensional space through the points whose coordinates are the elements of xopt3, yopt3 and -fopt3 (to get inverted plot)

title %Command used to give title to the plot

plot((-fopt3)) %To plot the (-fopt3)(to get inverted plot) vs iterations graph

xlabel %Command used to name the x-axis

ylabel %Command used to name the y-axis

Output of the step - Matrices inputs3 and fopt3 gets created after passing function, no of design variables, bound limits and options (i.e for modification of default functioning of ga, such as defining population size is equal to 200) as input arguments. Also optimum values of x and y (i.e xopt3 and yopt3) are obtained in the matrix form with 100 columns (1 element = 1 iteration). Time required to perform ga gets calulated. Plot of 3D contour and optimum value of function vs no of iterations are obtained in a single figure. Plots are more acurate as compared to earlier two.

7) Seventh Step - Calculation and printing of global maxima (An extra work)

Description -

max_value = max(abs(fopt3)) %Used to store the maximum value of the stalagmite function which is obtained by using \'max\' command and passing absolute value of \'fopt3\' as the input argument

maxima %Variable used to print value of variable \'max_value\' in the command window using \'sprintf\' command

PART II - SUBPROGRAM DEFINING STALAGMITE FUNCTION

Since, in MATLAB, genetic algorithm calculates the minimum value of the function (i.e minimizes function), we need to pass the negative value of the stalagmite function. As the value of negative of stalagmite function minimizes, the value of original stalagmite function maximizes.

1) First Step - Definition of stalagmite function

function %Default keyword in matlab to define a function

F %Output argument used to store the value of the function

input_vector %Input argument of the function

Output of the step - Input and output arguments of the function gets defined.

2) Second Step - Defining terms of the function

Description -

f1x, f2x, f1y, f2y %Terms required to calculate the value of the function

F %Variable used to store the value of function as the negative of product of f1x, f2x, f1y, f2y

Output of the step - The value of the terms f1x, f2x, f1y, f2y and the function gets calculated as negative of the product of terms f1x, f2x, f1y and f2y.

3) Third Step - Ending the function

Description -

end %Command to end the function

Output of the step - The function gets ended.

Output of the Code :-

Study 1 - Statistical behaviour with unbounded inputs

Study 2 - Statistical behaviour with bounded inputs

Study 3 - Statistical behaviour with bounded inputs and modified population size

Analytical Value of Global Maxima -

By seeing the above two figures, we can say that the analytical and graphical values of global maxima are very close and approximately equal. This validates the results.

Analytical Value of Global Maxima = 0.9999998

Graphical Value of Global Maxima = 1