Optimization of stalagmite function using Genetic Algorithm

Problem Statement :-

1. Write a code in MATLAB to optimize the stalagmite function and find the global maxima of the function.
2. Explain the concept of genetic algorithm and also explain the syntax for "ga" in MATLAB.
3. Make sure that your code gives the same output, even if it is made to run several times.
4. Plot graphs for studies and for F maximum v/s no. of iterations. Title and axes labels are a must, legends could be shown if necessary.

Theory :-

Genetic Algorithm

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 individuals are selected for reproduction in order to produce offspring of the next generation.

The genetic algorithm is a method for solving both constrained and unconstrained optimization problems that is based on natural selection, the process that drives biological evolution. The genetic algorithm repeatedly modifies a population of individual solutions. At each step, the genetic algorithm selects individuals at random from the current population to be parents and uses them to produce the children for the next generation. Over successive generations, the population "evolves" toward an optimal solution. The genetic algorithm can apply to solve a variety of optimization problems that are not well suited for standard optimization algorithms, including problems in which the objective function is discontinuous, nondifferentiable, stochastic, or highly nonlinear. The genetic algorithm can address problems of mixed integer programming, where some components are restricted to be integer-valued.

The genetic algorithm uses three main types of rules at each step to create the next generation from the current population:

• Selection rules select the individuals, called parents, that contribute to the population at the next generation.
• Crossover rules combine two parents to form children for the next generation.
• Mutation rules apply random changes to individual parents to form children.

Genetic Algorithm Terminology

Fitness Functions

The fitness function is the function you want to optimize. For standard optimization algorithms, this is known as the objective function. The toolbox software tries to find the minimum of the fitness function. Write the fitness function as a file or anonymous function, and then pass it as a function handle input argument to the main genetic algorithm function.

Individuals

An individual is any point to which you can apply the fitness function. The value of the fitness function for an individual is its score. For example, if the fitness function is

`f(x_1,x_2,x_3)=(2x_1+1)^2+(3x_2+4)^2+(x_3−2)^2`,

the vector (2, -3, 1), whose length is the number of variables in the problem, is an individual. The score of the individual (2, –3, 1) is f(2, –3, 1) = 51.

An individual is sometimes referred to as a genome and the vector entries of an individual as genes.

Populations and Generations

A population is an array of individuals. For example, if the size of the population is 100 and the number of variables in the fitness function is 3, then the population represented by a 100-by-3 matrix. The same individual can appear more than once in the population. For example, the individual (2, -3, 1) can appear in more than one row of the array. At each iteration, the genetic algorithm performs a series of computations on the current population to produce a new population. Each successive population is called a new generation.

Diversity

Diversity refers to the average distance between individuals in a population. A population has high diversity if the average distance is large; otherwise it has low diversity. In the following figure, the population on the left has high diversity, while the population on the right has low diversity. Diversity is essential to the genetic algorithm because it enables the algorithm to search a larger region of the space.

Fitness Values and Best Fitness Values

The fitness value of an individual is the value of the fitness function for that individual. Because the toolbox software finds the minimum of the fitness function, the best fitness value for a population is the smallest fitness value for any individual in the population.

Parents and Children

To create the next generation, the genetic algorithm selects certain individuals in the current population, called parents, and uses them to create individuals in the next generation, called children. Typically, the algorithm is more likely to select parents that have better fitness values.

Genetic Algorithm Works : -

Following are steps through the genetic algorithm works.

1. The algorithm begins by creating a random initial population.
2. The algorithm then creates a sequence of new populations. At each step, the algorithm uses the individuals in the current generation to create the next population. To create the new population, the algorithm performs the following steps:
1. Scores each member of the current population by computing its fitness value. These values are called the raw fitness scores.
2. Scales the raw fitness scores to convert them into a more usable range of values. These scaled values are called expectation values.
3. Selects members, called parents, based on their expectation.
4. Some of the individuals in the current population that have lower fitness are chosen as elite. These elite individuals are passed to the next population.
5. Produces children from the parents. Children are produced either by making random changes to a single parent—mutation—or by combining the vector entries of a pair of parents—crossover.
6. Replaces the current population with the children to form the next generation.
3. The algorithm stops when one of the stopping criteria is met.

‘ga’ Syntax in MATLAB

In MATLAB ‘ga’ is use for find minimum of function using genetic algorithm. Following syntaxes are use for ‘ga’ in MATLAB.

1. x = ga(fun,nvars) [This finds a local unconstrained minimum, x, to the objective function, fun. nvars is the dimension (number of design variables) of fun.]
2. x = ga(fun,nvars,A,b) [ This finds a local minimum x to fun, subject to the linear inequalities A*x ≤ b. ga evaluates the matrix product A*x as if x is transposed (A*x').]
3. x = ga(fun,nvars,A,b,Aeq,beq) [ This finds a local minimum x to fun, subject to the linear equalities Aeq*x = beq and A*x ≤ b. (Set A=[] and b=[] if no linear inequalities exist.) ga evaluates the matrix product Aeq*x as if x is transposed (Aeq*x').]
4. x = ga(fun,nvars,A,b,Aeq,beq,lb,ub) [ This 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. (Set Aeq=[] and beq=[] if no linear equalities exist.)]
5. x = ga(fun,nvars,A,b,Aeq,beq,lb,ub,nonlcon) [ This subjects the minimization to the constraints defined in nonlcon. The function nonlcon accepts x and returns vectors C and Ceq, representing the nonlinear inequalities and equalities respectively. ga minimizes the fun such that C(x) ≤ 0 and Ceq(x) = 0. (Set lb=[ ] and ub=[ ] if no bounds exist.)]
6. x = ga(fun,nvars,A,b,Aeq,beq,lb,ub,nonlcon,options) [ This subjects the minimization to the constraints defined in nonlcon. The function nonlcon accepts x and returns vectors C and Ceq, representing the nonlinear inequalities and equalities respectively. ga minimizes the fun such that C(x) ≤ 0 and Ceq(x) = 0. (Set lb=[ ] and ub=[ ] if no bounds exist.)]
7. x = ga(fun,nvars,A,b,[],[],lb,ub,nonlcon,IntCon) or x = ga(fun,nvars,A,b,[],[],lb,ub,nonlcon,IntCon,options [ This requires that the variables listed in IntCon take integer values.]
8. x = ga(problem) [ This finds the minimum for problem, where problem is a structure.]
9. [x,fval] = ga(___) [ for any previous input arguments, also returns fval, the value of the fitness function at x.]
10. [x,fval,exitflag,output] = ga(___) [ This also returns exitflag, an integer identifying the reason the algorithm terminated, and output, a structure that contains output from each generation and other information about the performance of the algorithm.
11. [x,fval,exitflag,output,population,scores] = ga(___) [ This also returns a matrix population, whose rows are the final population, and a vector scores, the scores of the final population.]

[In MATLAB Code used ga Syntax is explain in comment.]

INPUT ARGUMENTS:-

1. fun – Objective function (Given by function handle or by function name)
2. nvars – Number of Variables (Given by positive integer)
3. A – Linear inequality constraints (Given by real matrix)
4. b – Linear inequality constraints (Given by real matrix)
5. Aeq – Linear equality constraints (Given by real matrix)
6. beq – Linear equality constraints (Given by real matrix)
7. lb – Lower Bounds (Real vector or array)
8. ub – Upper Bounds (Real vector or array)
9. nonlcon – Nonlinear Constraints (function handle or function name)
10. options – Optimization options
11. IntCon – Integer variables (Vector of positive integers)
12. problem – Problem description(structure)

OUTPUT ARGUMENTS:-

1. x – Solution(real vector)
2. fval – Objective function value at solution(real number)
3. exitflag – Reason ga stopped (integer)
4. output – Information about the optimization process(structure)
5. population – Final population(matrix)
6. scores – Final scores(column vector)

STOPPING CONDITIONS FOR THE ALGORITHM

The genetic algorithm uses the following conditions to determine when to stop:

• Generations — The algorithm stops when the number of generations reaches the value of Generations.
• Time limit— The algorithm stops after running for an amount of time in seconds equal to Time limit.
• Fitness limit — The algorithm stops when the value of the fitness function for the best point in the current population is less than or equal to Fitness limit.
• Stall generations — The algorithm stops when the average relative change in the fitness function value over Stall generations is less than Function tolerance.
• Stall time limit — The algorithm stops if there is no improvement in the objective function during an interval of time in seconds equal to Stall time limit.
• Stall test — The stall condition is either average change or geometric weighted. For geometric weighted, the weighting function is 1/2ⁿ, where n is the number of generations prior to the current. Both stall conditions apply to the relative change in the fitness function over Stall generations.
• Function tolerance — The algorithm runs until the average relative change in the fitness function value over Stall generations is less than Function tolerance.
• Constraint tolerance — The Constraint tolerance is not used as stopping criterion. It is used to determine the feasibility with respect to nonlinear constraints. Also, max(sqrt(eps),ConstraintTolerance)determines feasibility with respect to linear constraints.

The algorithm stops as soon as any one of these conditions is met.

Stalagmite Function

Stalagmite function is given by following equations.

`f(x,y) = f_{1,x} f_{2,x} f_{1,y} f_{2,y}`

Where :-

`f_{1,x} = [sin(5.1pix+0.5)]^6`

`f_{2,x} = exp [-4 ln(2) frac{(x-0.0667)^2}{0.64}]`

`f_{1,y} = [sin(5.1piy+0.5)]^6`

`f_{2,y} = exp [-4 ln(2) frac{(y-0.0667)^2}{0.64}]`

By these equations stalagmite function is define. We using this equation to define Stalagmite function in MATLAB.

MATLAB Code : -

1. MATLAB Code for Stalagmite Function :-

% Define function f_xy is output, function name is stalagmite and input_vector is input variable.
function  f_xy = stalagmite(input_vector)

% Provide value of x and y through input_vector
A = input_vector(1);   % x is taken as A
B = input_vector(2);   % y is yaken as B

% Define f_1x and f_2x function according to formula of Stalagmite function
f_1x = (sin((5.1*pi*A) + 0.5))^6 ; 
f_2x = exp(-4* log(2)*(((A-0.0667)^2) / 0.64));

% Define f_1y and f_2y function according to formula of Stalagmite function
f_1y = (sin((5.1*pi*B) + 0.5))^6 ; 
f_2y = exp(-4* log(2)*(((B-0.0667)^2) / 0.64));

% Define output function f(x,y)
f_xy = f_1x * f_2x * f_1y * f_2y;   % This is maximum case

% Minimize function from maximization
f_xy = -f_xy; % Because when calculating "ga" that time ga take only minimum values
end

% Noted:- this file should be save as "stalagmite.m" same as function name

2. MATLAB Code for Optimize Stalagmite Function and Find Global Maxima

clear all   % Clear the Workspace
close all   % Clear all figures and plot
clc         % Clear the Command Window

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

% Input the number how many time genetic algorithm run
iteration = 50;

% Creating a 2-Dimensional Mesh 
[xx, yy] = meshgrid(x,y);

% Evaluate the stalagmite 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(1) 
% Generate stalagmite function as a surfc plot
surfc(xx,yy,-f) 

% Title and Label
xlabel('x-value')
ylabel('y-value')
zlabel('Function Value')  
title('Stalagmite Function')

shading interp % Use for Generate more smoother plot

rng default % This command is random number generator which is use to get same value of ga for every run

% Study 1 :- Unbounded Statistical Behavior

tic % Calculating time taken in study 
for i=1:iteration
    % Passing 'ga' function
    % Here syntax use for 'ga' is      [x,fval] = ga(fun,nvars)
    % Where fun is objective function(Here stalagmite function is use as objective function)
    % nvars is number of variable. fval is value of fitness function at x.
    [inputs1,fopt1(i)]=ga(@stalagmite,2);
    
    xopt1(i)=inputs1(1); % give optimum value to x
    yopt1(i)=inputs1(2); % give optimum value to y
end
study1_time=toc;

figure(2)
% Use subplot function to generate multiple figure in same window
subplot(2,1,1) 
hold on
surfc(xx,yy,-f)  % Generate surfc plot 
shading interp    % Generate more smoother plot

% By using plot3 we can plots multiple sets of coordinates on the same set of axes
plot3(xopt1,yopt1,-fopt1,'marker','o','markersize',5,'markerfacecolor','r')
% labelling and title
xlabel('x-value')
ylabel('y-value')
zlabel('Function value')
title('Unbounded Statistical Behvaiour')

% Second subplot
subplot(2,1,2)
% Plot optimum value
plot(fopt1,'linewidth',1.5)  
xlabel('Iterations')
ylabel('Function Maximum')

% Study 2 - Statistical behavior with Upper and Lower Bounds
tic
for i=1:iteration
    % Passing 'ga' function
    % Here syntax use for 'ga' is      [x,fval] = ga(fun,nvars,A,b,[],[],lb,ub)
    % This 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.
    % Where fun is objective function(Here stalagmite function is use as objective function)
    % nvars is number of variable. fval is value of fitness function at x.
    % lb is lower bound and ub is upper bound
    
    [inputs2,fopt2(i)]=ga(@stalagmite,2,[],[],[],[],[0;0],[1;1]);
    
    xopt2(i)=inputs2(1); % give optimum value to x
    yopt2(i)=inputs2(2); % give optimum value to y
end
study2_time=toc;

figure(3)
% Use subplot function to generate multiple figure in same window
subplot(2,1,1) 
hold on 
surfc(xx,yy,-f)  % Generate surfc plot 
shading interp    % Generate more smoother plot

% By using plot3 we can plots multiple sets of coordinates on the same set of axes
plot3(xopt2,yopt2,-fopt2,'marker','o','markersize',5,'markerfacecolor','r')
% labelling and title
xlabel('x-value')
ylabel('y-value')
zlabel('Function value')
title('Bounded Statistical Behvaiour')

% Second subplot
subplot(2,1,2)
% Plot optimum value
plot(-fopt2,'linewidth',1.5)  
xlabel('Iterations')
ylabel('Function Maximum')

% Study 3 Increasing Population size for Statistical behavior calculation 

% By using optimoptions we create option for optimize 'ga'.
options = optimoptions('ga');
options = optimoptions(options,'PopulationSize',200);

for i=1:iteration
    % Passing 'ga' function
    % Here syntax use for 'ga' is      [x,fval] = ga(fun,nvars,A,b,[],[],lb,ub,nonlcon,IntCon,options)
    % This 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.
    % This Syntax minimize with the default optimization parameters replaced by the values in options.
    % Where fun is objective function(Here stalagmite function is use as objective function)
    % nvars is number of variable. fval is value of fitness function at x.
    % lb is lower bound and ub is upper bound
    
    [inputs3,fopt3(i)]=ga(@stalagmite,2,[],[],[],[],[0;0],[1;1],[],[],options);
    
    xopt3(i)=inputs3(1); % give optimum value to x
    yopt3(i)=inputs3(2); % give optimum value to y
end
study3_time=toc;

figure(4)
% Use subplot function to generate multiple figure in same window
subplot(2,1,1) 
hold on
surfc(xx,yy,-f)  % Generate surfc plot 
shading interp    % Generate more smoother plot

% By using plot3 we can plots multiple sets of coordinates on the same set of axes
plot3(xopt3,yopt3,-fopt3,'marker','o','markersize',5,'markerfacecolor','r')
% labelling and title
xlabel('x-value')
ylabel('y-value')
zlabel('Function value')
title('Bounded Statistical Behvaiour With Modified Population Size')

% Second subplot
subplot(2,1,2)
% Plot optimum value
plot(-fopt3,'linewidth',1.5)  
xlabel('Iterations')
ylabel('Function Maximum')


% Consolidated Results

fprintf ('\n Study 1 :- Unbounded Statistical Behavior')
fprintf ('\n The Time taken for the study 1 is : %f ',study1_time)
fprintf ('\n The Co-ordinates of x and y for the Global Maxima of the Function is : %f,%f',inputs1(1),inputs1(2))
fprintf ('\n The Global Maxima of the Function is : %f',-(stalagmite(inputs1)))  

fprintf ('\n\n Study 2 - Statistical behavior with Upper and Lower Bounds')
fprintf ('\n The Time taken for the study 2 is : %f',study2_time)
fprintf ('\n The Co-ordinates of x and y for the Global Maxima of the Function is : %f,%f',inputs2(1),inputs2(2))
fprintf ('\n The Global Maxima of the Function is : %f',-(stalagmite(inputs2)))

fprintf ('\n\n Study 3 Increasing Population size for Statistical behavior calculation ')
fprintf ('\n The Time taken for the study 3 is : %f',study3_time)
fprintf ('\n The Co-ordinates of x and y for the Global Maxima of the Function is : %f,%f',inputs3(1),inputs3(2))
fprintf ('\n The Global Maxima of the Function is : %f',-(stalagmite(inputs3)))

Important point in code making:-

• First a function stalagmite is created. This function is created by using above equations of stalagmite function.
• This equation give the maximum value of function.
• But it is known that genetic algorithm optimize minimum value of function for that we converted this function in minimum function by a simple rule : `max f = -min f`
• Now the optimization is take care in program.
• For optimization predefine MATLAB function 'ga' is use. First unbounded statistical behaviour was studied after that bounded statistical behaviour was studied. For get more accurate result more bounded study with increased population size is taken into account.
• It is known that this result generate minimum function value but for solution maximum value is requaired. For achive maximum value, minimum value is multiolied by -1.
• By all this study of stalagmite function is achived.

Result:-

1. Figure 1

2. Figure 2

3. Figure 3

4. Figure 4

5. Output : (In command window)

Study 1 :- Unbounded Statistical Behavior
The Time taken for the study 1 is : 73.564228
The Co-ordinates of x and y for the Global Maxima of the Function is : 0.278594,0.065924
The Global Maxima of the Function is : 0.679370

Study 2 - Statistical behavior with Upper and Lower Bounds
The Time taken for the study 2 is : 52.041066
The Co-ordinates of x and y for the Global Maxima of the Function is : 0.066842,0.066837
The Global Maxima of the Function is : 1.000000

Study 3 Increasing Population size for Statistical behavior calculation
The Time taken for the study 3 is : 182.171378
The Co-ordinates of x and y for the Global Maxima of the Function is : 0.066833,0.066831
The Global Maxima of the Function is : 1.000000

References:-

1. https://in.mathworks.com/help/gads/ga.html
2. https://www.researchgate.net/profile/Mohamed_Mourad_Lafifi/post/How_can_we_use_Genetic_Algorithm_in_Matlab_such_that_it_returns_same_answer_every_time_we_run_it/attachment/5aa19b874cde266d589024da/AS%3A602004858626049%401520540551775/download/An+Introduction+to+Genetic+Algorithms+Carrjk.pdf
3. https://in.mathworks.com/help/gads/some-genetic-algorithm-terminology.html
4. https://in.mathworks.com/help/gads/how-the-genetic-algorithm-works.html
5. https://towardsdatascience.com/introduction-to-genetic-algorithms-including-example-code-e396e98d8bf3?gi=4b3f2cf1bdf4
6. https://in.mathworks.com/help/gads/what-is-the-genetic-algorithm.html