MATLAB function to optimize the stalagmite function using genetic algorithm and to find the global maxima.

I. AIM: To optimise the stalagmite function and find the global maxima.

II. OBJECTIVES:

• - To understand how genetic algorithm (GA) works
• - To find the global maxima of the stalagmite function

III. THEORY:

1.) Stalagmite function: A stalagmite is a type of speleothem that rises from the floor of a limestone cave due to the dripping of mineralized solutions and the deposition of calcium carbonate.[1] Scientists have modeled this formation by developing a two variable function called the stalagmite function. The mathematical function can be expressed as,

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

where,

`f_(1,x) = [sin(5.1.pi.x + 0.5)]^6`

`f_(2,x) = exp[-4.ln(2).(x-0.0667)^2/0.64]`

`f_(1,y) = [sin(5.1.pi.y + 0.5)]^6`

`f_(2,y) = exp[-4.ln(2).(y-0.0667)^2/0.64]`

2.) Genetic Algorithm (GA): Genetic algorithm (GA) is a method for solving both constrained and unconstrained optimization problems based on a natural selection process that mimics biological evolution. The algorithm repeatedly modifies a population of individual solutions. At each step, the genetic algorithm randomly selects individuals from the current population and uses them as parents to produce the children for the next generation. Over successive generations, the population "evolves" toward an optimal solution. You can apply the genetic algorithm to solve problems that are not well suited for standard optimization algorithms, including problems in which the objective function is discontinuous, nondifferentiable, stochastic, or highly nonlinear. [2]

a.) working:

The following outline summarizes how the genetic algorithm works:

figure (1.a) Working of Genetic Algorithm. Source: Genetic Algorithms By Anas Amjad Obeidat Advanced Algorithms 02 Semester 2 - 2008/2009 March 18 - 2009

figure (1.b): Working of GA.

b.) Syntax:

X = ga([1],[2],[3],[4],[5],[6],[7],[8],[9],[10])

[1]: Fitness Function (FITNESSFCN). FITNESSFCN accepts a vector X of size 1-by-NVARS, and returns a scalar evaluated at X.

[2]: Number of design variables (NVARS). NVARS is the dimension of the FITNESSFCN.

[3] and [4]: Linear inequalities (A & b). Finds a local minimum X to the function FITNESSFCN, subject to the linear inequalities A*X <= B. Linear constraints are not satisfied when the PopulationType option is set to 'bitString' or 'custom'.

[5] and [6]: Linear constraints (Aeq & beq). Finds a local minimum X to the function FITNESSFCN, subject to the linear equalities Aeq*X = beq as well as A*X <= B. (Set A=[] and B=[] if no inequalities exist.) Linear constraints are not satisfied when the PopulationType option is set to 'bitString' or 'custom.

[7] and [8]: Lower and upper bound respectively (ub & lb). Defines a set of lower and upper bounds on the design variables, X, so that a solution is found in the range lb <= X <= ub. Use empty matrices for lb and ub if no bounds exist. Set lb(i) = -Inf if X(i) is unbounded below; set ub(i) = Inf if X(i) is unbounded above. Linear constraints are not satisfied when the PopulationType option is set to 'bitString' or 'custom'.

[9]: Non-linear constraints (NONLCON). Subjects the minimization to the constraints defined in NONLCON. The function NONLCON accepts X and returns the vectors C and Ceq, representing the
nonlinear inequalities and equalities respectively. ga minimizes FITNESSFCN such that C(X)<=0 and Ceq(X)=0. (Set lb=[] and/or ub=[] if no bounds exist.) Nonlinear constraints are not satisfied when the PopulationType option is set to 'bitString' or 'custom'.

[10]: Options. Minimizes with the default optimization parameters replaced by values in OPTIONS. OPTIONS can be created with the OPTIMOPTIONS function.

Complete Syntax: X = ga(FITNESSFCN, NVARS, A, b, Aeq, beq, lb, ub, NONLCON, options) [4]

IV. PROGRAM:

1.) Stalagmite function:

1. We store the two rows of the input vector ('input_vector') into two variables 'x' and 'y'.
2. The two inputs are passed through the stalagmite function and the output is stored in the variable 'f'.

function f = stalagmite(input_vector)

x = input_vector(1);
y = input_vector(2);
f = -(sin(5.1*pi*x + 0.5)^6)*exp(-4*log(2)*(((x - 0.0667)^2)/0.64))*(sin(5.1*pi*y + 0.5)^6)*exp(-4*log(2)*(((y - 0.0667)^2)/0.64));

end

Note: Here, a minus sign is given at the starting of the function as GA only finds the minima, hence if we use GA it finds the local minima of the stalagmite function which is flipped and once this is reverted, the calculated local minima becomes the local maxima.

2.) Optimization:

1. An array of inputs are stored in the variable 'x' and 'y'. The vector is then converted into a matrix by using the 'meshgrid' function.
2. A nested loop is created containing two 'for' loops which stores each element of the matrices 'xx' and 'yy' in the 2 row vector 'input_vector'.
3. This vector is then passed as an argument into the 'stalagmite.m' function and the outputs are stored in the variable 'f'. 'f' is a matrix of dimensions 'length(xx)*length(yy)'.
4. Now three studies are conducted:
   • Study 1: The inputs that are given to GA are not bounded, i.e, a upper and lower bound hasn't been mentioned. Since, GA conducts the optimization by Random Sampling the answers calculated will be different when the program is executed again. The variation of the minimum of the function is very significant. See figure (2.a).
   • Study 2: The inputs that are given to GA, now have a lower and upper bound. This limits GA to find the local maxima within the function range. Now, the variation of the minimum of the function is reduced but it is still significant as the local maxima doesn't converge with the global maxima. See figure (2.b).
   • Study 3: GA by default takes 50 as the population size. Here, we increase the population size to 170 along with upper and lower bounds. Now, the variation of minimum of the function is negligible and the local maxima converges with the global maxima. The population size is increased by using 'optimoptions'. See figure (2.c).
5. The time of each study is recorded in three different variables.
6. Now, two subplots are created for each study. Sub plot-1 shows the stalagmite function and the point of local maxima. Here, we can see wether the local maxima converges with the global maxima. And Sub plot-2 show the variation of the minimum of the function with respect to the number of iterations (here the number of iterations is set to 50).
7. The results are recorded, analysed and plotted.

clear all
close all
clc

%% Inputs
x = linspace(0,0.6,500);
y = linspace(0,0.6,500);
num_cases = 50;                     % Number of iterations
[xx yy] = meshgrid(x,y);

%% Calculate the stalagmite function and store it in the variable f
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


%% Study-1 | Unbound inputs
tic
for i = 1:num_cases
    [inputs, fopt(i)] = ga(@stalagmite, 2);
    xopt(i) = inputs(1);
    yopt(i) = inputs(2);
end
study_time1 = toc;

figure(1)
subplot(2,1,1)
hold on
surfc(xx,yy,-f)
shading interp
plot3(xopt,yopt,-fopt,'marker', 'o', 'markersize', 5, 'markerfacecolor', 'r')
title('Unbound Inputs');
subplot(2,1,2)
plot(-fopt)
xlabel('Number of Iterations');
ylabel('Function maximum');
hold on

%% Study-2 | Bound inputs
tic
for i = 1:num_cases
    [inputs, fopt(i)] = ga(@stalagmite, 2, [], [], [], [], [0,0],[0.1,0.1]);
    xopt(i) = inputs(1);
    yopt(i) = inputs(2);
end
study_time2 = toc;

figure(2)
subplot(2,1,1)
hold on
surfc(x,y,-f)
shading interp
plot3(xopt,yopt,-fopt,'marker', 'o', 'markersize', 5, 'markerfacecolor', 'r')
title('Bound Inputs');
subplot(2,1,2)
plot(-fopt)
xlabel('Number of Iterations');
ylabel('Function maximum');

%% Study-3 | Bound Inputs plus increase in population size
options = optimoptions('ga');
options = optimoptions(options, 'PopulationSize', 170);

tic
for i = 1:num_cases
    [inputs, fopt(i)] = ga(@stalagmite, 2, [], [], [], [], [0,0],[0.1,0.1], [], [], options);
    xopt(i) = inputs(1);
    yopt(i) = inputs(2);
end
study_time3 = toc;

figure(3)
subplot(2,1,1)
hold on
surfc(x,y,-f)
shading interp
plot3(xopt,yopt,-fopt,'marker', 'o', 'markersize', 5, 'markerfacecolor', 'r')
title('Bound Inputs with increase in population size');
subplot(2,1,2)
plot(-fopt)
xlabel('Number of Iterations');
ylabel('Function maximum');

V. OUTPUTS:

VI. PLOTS:

figure (2.a) Unbound inputs.

figure (2.b) Bound inputs.

figure (2.c) Bound inputs with increase in population size.

VII. RESULTS:

1. Global maxima of the stalagmite function is found. The value is 1 and the global maxima is found at [0.066831780971505,0.066830233700099].
2. The time taken for each study is shown below,
   • Study 1: 6.6996
   • Study 2: 8.1094
   • Study 3: 21.6179

VIII. REFERENCES:

[1]-https://www.sciencedaily.com/terms/stalagmite.htm#:~:text=Stalagmite,is%20known%20as%20a%20stalactite.

[2]-https://in.mathworks.com/discovery/genetic-algorithm.html

[3]-https://in.mathworks.com/help/gads/how-the-genetic-algorithm-works.html#:~:text=The%20following%20outline%20summarizes%20how,to%20create%20the%20next%20population.

[4]-Using the 'help' command in the MATLAB command window. [ help ga ]