Evaluating the Maximum value of Stalagmite function and the point at which occurs using GENETIC ALGORITHM in MATLAB

OPTIMISATION USING GENETIC ALGORITHM

1. Problem Statement/Objective

Optimizing Stalagmite function and finding the maximum value of it and the point at which it occurs using Genetic Algorithm in MATLAB.

1. Introduction

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. You can apply the genetic algorithm 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 rulesselect the individuals, called parents, that contribute to the population at the next generation.
• Crossover rulescombine two parents to form children for the next generation.
• Mutation rulesapply random changes to individual parents to form children.

• Intuition behind Genetic Algorithm

Let’s start by a famous quote given by Charles Darwin;

“It is not the strongest species that survives, nor the most intelligent, but the one that is more responsive to change”

The entire concept of GA is based on the above line. Let us understand this with the help of an example;

Let’s take a hypothetical situation where, you are head of a country, and in order to keep your safe from bad things, you implement a policy like this.

• You select all the good people, and ask them to extend their generation by having their children.
• This repeats for a few generations.
• You will notice that now you have an entire population of good people.

Now, that may not be entirely possible, but this example was just to help you understand the concept. So, the basic idea was that we changed the input (i.e. population) such that we get better output (i.e. better country). Therefore, the concept of a genetic algorithm is somewhat related to biology.

• Biological Inspiration

“Cells are the basic building block of all living things”

Therefore, in each cell, there is the same set of chromosomes. Chromosome are basically the strings of DNA.

Traditionally, these chromosomes are represented in binary as strings of 0’s and 1’s.

A chromosome consists of genes, commonly referred as blocks of DNA, where each gene encodes a specific trait, for example hair color or eye color. Now, we have recalled these basics concept of biology, now let’s get back and understand what actually is a genetic algorithm?

2.3 Genetic Algorithm

Let’s get back to the example we discussed above and summarize what we did.

1. Firstly, we defined our initial population as our countrymen.
2. We defined a function to classify whether is a person is good or bad.
3. Then we selected good people for mating to produce their off-springs.
4. And finally, these off-springs replace the bad people from the population and this process repeats.

This is how genetic algorithm actually works, which basically tries to mimic the human evolution to some extent. Thus, to formalize a definition of a genetic algorithm, we can say that it is an optimization technique, which tries to find out such values of input so that we get the best output values or results. The working of a genetic algorithm is also derived from biology, which is as shown in the image below.

2.4 Steps of Genetic Algorithm

Here, to make things easier, let us understand it by the famous problem. If you haven’t come across this problem, let me introduce my version of this problem.

Let’s say, you are going to spend a month in the wilderness. Only thing you are carrying is the backpack which can hold a maximum weight of 30 kg. Now you have different survival items, each having its own “Survival Points” (which are given for each item in the table). So, your objective is maximizing the survival points.

Here is the table giving details about each item.

2.5 Initialization

To solve this problem using genetic algorithm, our first step would be defining our population. Therefore, our population will contain individuals, each having their own set of chromosomes.

We know that, chromosomes are binary strings, where for this problem 1 would mean that the following item is taken and 0 meaning that it is dropped.

This set of chromosome is considered as our initial population.

2.6 Fitness Function

Let us calculate fitness points for our first two chromosomes.

For A1 chromosome [100110],

Similarly, for A2 chromosome [001110],

So, for this problem, our chromosome will be considered as more fit when it contains more survival points.

Therefore, chromosome 1 is more fit than chromosome 2.

2.7 Selection

Now, we can select fit chromosomes from our population which can mate and create their off-springs.

General thought is that we should select the fit chromosomes and allow them to produce off-springs. But that would lead to chromosomes that are closer to one another in a few next generations, and therefore less diversity. Therefore, we generally use Roulette Wheel Selection method.

Consider a wheel, and let’s divide that into m divisions, where m is the number of chromosomes in our populations. The area occupied by each chromosome will be proportional to its fitness value.

Based on these values, let us create our roulette wheel.

So, now this wheel is rotated and the region of wheel which comes in front of the fixed point is chosen as the parent. For the second parent, the same process is repeated.

Sometimes we mark two fixed point as shown in the figure below.

So, in this method we can get both our parents in one go. This method is known as Stochastic Universal Selection method.

2.8 Cross-Over

In this previous step, we have selected parent chromosomes that will produce off-springs. Thus, in biological terms, crossover is nothing but reproduction.

Now, let us find the crossover of chromosome 1 and 4, which were selected in the previous step. Take a look at the image below.

This is the most basic form of crossover, known as one-point crossover. Here we select a random crossover point and the tails of both the chromosomes are swapped to produce a new off-springs.

If you take two crossover point, then it will be called as multi point crossover which is as shown below.

2.9 Mutation

Now if you think in the biological sense, are the children produced have the same traits as their parents? The answer is NO. During their growth, there is some change in the genes of children which makes them different from its parents. This process is known as mutation, which may be defined as a random tweak in the chromosome, which also promotes the idea of diversity in the population.

A simple method of mutation is shown in the image below.

Thus, the entire process is summarized as shown in the figure.

The off-springs thus produced are again validated using our fitness function, and if considered fit then will replace the less fit chromosomes from the population.

But the question is how we will get to know that we have reached our best possible solution?

So basically, there are different termination conditions, which are listed below:

• There is no improvement in the population for over x iterations.
• We have already predefined an absolute number of generations for our algorithm.
• When our fitness function has reached a predefined value.

1. GA Syntax in MATLAB

3.1 Syntax

x = ga(fun,nvars)

Description: Finds a local unconstrained minimum, x, to the objective function, fun. nvars is the dimension (number of design variables) of fun.

x = ga(fun,nvars,A,b)

Description: 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\').

x = ga(fun,nvars,A,b,Aeq,beq)

Description: 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\').

x = ga(fun,nvars,A,b,Aeq,beq,lb,ub)

Description: 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.)

x = ga(fun,nvars,A,b,Aeq,beq,lb,ub,nonlcon)

Description: Subjects the minimization to the constraints defined in nonlcon. The function nonlcon accepts x and returns vectors Cand 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.)

x = ga(fun,nvars,A,b,Aeq,beq,lb,ub,nonlcon,options)

Description: Minimizes with the default optimization parameters replaced by values in options. (Set nonlcon=[] if no nonlinear constraints exist.) Create options using optimoptions.

x = ga(fun,nvars,A,b,[],[],lb,ub,nonlcon,IntCon) and

x = ga(fun,nvars,A,b,[],[],lb,ub,nonlcon,IntCon,options)

Description: Requires that the variables listed in IntCon take integer values.

x = ga(problem)

Description: Finds the minimum for problem, where problem is a structure.

[x,fval] = ga(___)

Description: For any previous input arguments, also returns fval, the value of the fitness function at x.

[x,fval,exitflag,output] = ga(___)

Description: Also returns exit flag, 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.

[x,fval,exitflag,output,population,scores] = ga(___)

Description: Also returns a matrix population, whose rows are the final population, and a vector scores, the scores of the final population.

3.2 Options and Outputs

To run the genetic algorithm with the Default options, call GA with the syntax

[x,fval] = ga(@fitnessfun, nvars)

The input arguments to GA are

• @fitnessfun — A function handle to the file that computes the fitness function.
• nvars — The number of independent variables for the fitness function.

The output arguments are

• x — The final point
• fval — The value of the fitness function at x

You can specify any of the options that are available for ga by passing Options as an input argument to GA using the syntax

[x,fval] = ga(@fitnessfun,nvars,[],[],[],[],[lb],[ub],[],options)

NOTE: This syntax does not specify any linear equality, linear inequality, or nonlinear constraints.

You create options using the function optimoptions.

options = optimoptions(@ga);

This returns options with the default values for its fields. GA uses these default values if you do not pass in options as an input argument.

The value of each option is stored in a field of options, such as options.PopulationSize. You can display any of these values by entering options followed by a period and the name of the field. For example, to display the size of the population for the genetic algorithm, enter

options.PopulationSize ans = \'50 when numberOfVariables <= 5, else 200\'

To create options with a field value that is different from the default

For example, to set PopulationSize to 100 instead of its default value 50 — enter

options = optimoptions(\'ga\',\'PopulationSize\',100);

This creates options with all values set to their defaults except for PopulationSize, which is set to 100.

If you now enter,

ga(@fitnessfun,nvars,[],[],[],[],[],[],[],options)

GA runs the genetic algorithm with a population size of 100.

1. MATLAB Study and CODE

The following is the Stalagmite Function:

There are two functions that were developed in my study, one was to evaluate the function for different values of x and y and another was to call upon to get minimised/maximised.

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.

% This function will evaluate the value of stalagmite function for the
% given range of inputs.

function [s] =  evaluatestalagmite(p,q)
         
   for i = 1 : length(p)
       for j = 1 : length(q)
           
           m = p(i);
           n = q(j);
           
           A = ((sin((5.1*pi*m)+ 0.5))^6);                    % F1,x
         
           B = exp(((-4)*log(2)*((m - 0.0667)^2))/0.64);      % F2,x
         
           C = ((sin((5.1*pi*n)+ 0.5))^6);                    % F1,y
         
           D = exp(((-4)*log(2)*((n - 0.0667)^2))/0.64);      % F2,y
         
           s(i,j) = -(A*B*C*D);                               % F(x,y)
           
       end
   end
end

% This is the modified stalagmite function (-ve sign) that has to be minimised using GA
% and the point at which it is minima will mean that the stalagmite function is maximum
% at that point

function S = stalagmite(matrix)                             
       
         m = matrix(1);
         n = matrix(2);
         
         A = ((sin((5.1*pi*m)+ 0.5))^6);                    % F1,x
         
         B = exp(((-4)*log(2)*((m - 0.0667)^2))/0.64);      % F2,x
         
         C = ((sin((5.1*pi*n)+ 0.5))^6);                    % F1,y
         
         D = exp(((-4)*log(2)*((n - 0.0667)^2))/0.64);      % F2,y
         
         S = -(A*B*C*D);                                    % F(x,y)
         
end

The following is the code that is used to calculate the minima point for modified stalagmite function, using GA;

clear all
close all
clc

% Defining inputs

x = linspace(0,0.6,150);
y = linspace(0,0.6,150);
cycles = 50;

%We have defined a new function \"evaluatestalagmite\" to evaluate the value of
%stalagmite function at each value of x and y.

f = evaluatestalagmite(x,y);

% Study 1 - Unbounded Inputs

tic
for i = 1 : cycles
    
    [point1,fmin1(i)] = ga(@stalagmite,2);
             xmin1(i) = point1(1);
             ymin1(i) = point1(2);
                
end
Duration_Study1 = toc

figure(1)
subplot(2,1,1)
hold on
surfc(x,y,f)
shading interp
plot3(xmin1,ymin1,fmin1,\'marker\',\'o\',\'markersize\',5,\'markerfacecolor\',\'g\')
subplot(2,1,2)
plot(fmin1)
xlabel(\'Iterations\')
ylabel(\'Minimum Function Value (WRONG RESULTS)\')



% Study 2 - GA with bounded Inputs

tic
for i = 1 : cycles
    
    [point2,fmin2(i)] = ga(@stalagmite,2,[],[],[],[],[0 0],[0.6 0.6]);
             xmin2(i) = point2(1);
             ymin2(i) = point2(2);
             
end

Duration_Study2 = toc

figure(2)
subplot(2,1,1)
hold on
surfc(x,y,f)
shading interp
plot3(xmin2,ymin2,fmin2,\'marker\',\'o\',\'markersize\',5,\'markerfacecolor\',\'g\')
subplot(2,1,2)
plot(fmin2)
xlabel(\'Iterations\')
ylabel(\'Minimum Function Value\')


% Study 3 - This Study is conducted to find the actual minima point and
% this is achieved by increasing the population size from default 50 to
% higher number, so that more offsprings are produced in each more produced
% generation, and hence the results are more refined. We will increase it
% to 186 using optimoptions.

options = optimoptions(@ga);
options = optimoptions(\'ga\',\'PopulationSize\',186);

tic
for i = 1 : cycles
    
    [point3,fmin3(i)] = ga(@stalagmite,2,[],[],[],[],[0 0],[0.6 0.6],[],options);
             xmin3(i) = point3(1);
             ymin3(i) = point3(2);
                
end
Duration_Study3 = toc


figure(3)
subplot(2,1,1)
hold on
surfc(x,y,f)
shading interp
plot3(xmin3,ymin3,fmin3,\'marker\',\'o\',\'markersize\',5,\'markerfacecolor\',\'g\')
subplot(2,1,2)
plot(fmin3)
xlabel(\'Iterations\')
ylabel(\'Minimum Function Value\')


% We now find the max value of function and the point at which it occurs

maxvalue = (-(fmin3(3)));

display(\'The max value of Stalagmite function is :\')
display(maxvalue)

display(\'And the point at which this occurs is :\')
display(point3)

      
    
        
        

STUDY 1

In this Study we just do 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.

STUDY 2

In this Study, we set the 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.85 and -1. To get more refined answer we will increase the population size in the next study.

STUDY 3

Since, the default population size that GA considers is 50, we are getting absurd results as in above study. To get more proper and considerable result we will increase the population size of GA using optimoptions so that more generations and offsprings are produced that will lead the result towards convergence. We will set the population size to 186. Now, in the results we can see that the minimum value occurs between -0.9999999985 and -0.9999999995 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 (x,y).

Therefore, the max value of original stalagmite function (without the negatice sign) occurs at 0.0668, 0.0668. And the maximum value is 1.

Thus the function is optimised.

5. References

• https://in.mathworks.com/help/gads/what-is-the-genetic-algorithm.html
• https://www.analyticsvidhya.com/blog/2017/07/introduction-to-genetic-algorithm/