1. Genetic Algorithm

A Genetic Algorithm is an optimisation technique which is used to minimise a certain function. When used in optimisation, generally its purpose is to minimise the error of a function. It is often used in curve fitting. This algorithm is named 'Genetic', because of its structure. It is a ' fittest of the survival' algorithm.

How it works:

In Genetic algorithm (GA), 'n' values form the function/problem are randomly selected (n is the population size). The values are selected after being sorted from the best to the worst ones first ('best' condition is subjected to the need). The GA selects the best and rejects the other values. Now, each selected value is mutated (n-1) times, i.e., some randomness is added to the values. Then, the mutated values are again checked against the minimum criteria condition. This iteration is continued up to the number of iterations decided.

For example, if we need to check the least error condition for the fitness of a curve, we will use GA to optimise the fitness. First, we will furnish the algorithm (say) 100 values, and this becomes the population size. Then we will sort these values from the least to the most error. Now, we will select the first 10 values only. Now, these values are mutated (i.e., randomness is added). Each value gives 9 more mutations. So, again there are (old 10 + new 9*10) 100 values, which is the initial population size. Again, the values are sorted, the best ten values are selected then tested over and over again. The testing will continue until m generations. The value of 'm' is decided by the optimizer.

In the above example, we see that only the best values are selected. Each value is mutated. Then, the values are further tested. If the old values survive, they will pass on. If not, they will be replaced with the mutated values. Hence, this algorithm follows the idea of Survival of the Fittest.

GA in MatLab:

In MatLab, the function for a genetic algorithm is 'ga'. It always tries to minimise the given fitness function. So, GA can be used when a certain function needs to be minimised. For example, minimising some error.

The syntax for the 'ga' is as:

[input, fval]= ga (fitnessfcn, nvars, A, b, Aeq, beq, LB, UB, nonlcon, IntCon, options).

Where: 'input' returns the input values.

'fval' returns the value of the fitness function.

'fitnessfcn' is the function which we desire to minimise.

'nvars' are the number of variables associated with the function.

'A' and 'b' is the linear inequality constraint of the type A*x ≤ b, if any. (In MatLab, it's always 'Ax is less than b' form).

'Aeq' and 'beq' is the equality constraint of the type A*x = b, if any.

LB and UB set the lower and upper bound for the GA. These provide the GA the search-space. Otherwise, the GA will try to look at points beyond the function's domain.

'nonlcon' sets the non-linear constraint, if any.

'IntCon' sets the integer constraint, if any.

'options' are used whenever we need to change the way of the Matlab's default GA functioning. to the desired optimisation requirements. For example, MaxGenerations is used to change the default no. of generations that the GA uses before stopping the search. (generally 100*number of variables). Similarly, the 'populationSize' is used when we need to change the size of the population to be considered (default size = 50 for 5 or less no of input variables).

Note: When there are no constraints and say we need to change the upper and lower bounds, we need to provide [] for each entry that comes before the UB and LB in the syntax of GA function. For example in the given problem, we need to change the upper and lower bounds only. The syntax becomes:

[inputs, Fopt(i)] = ga (@stagfun, 2, [], [], [], [], [0;0], [0.6;0.6])

and not [inputs, Fopt(i)] = ga (@stagfun, 2, [0;0], [0.6;0.6]) XX

2. Code For Finding The Global Maxima of The Stalagmite Function

The code to find the global maxima of the stalagmite function in the x = (0,0.6) and y=(0,0.6) domain.

1. Code for the stalagmite function:

function[output]=stagfun(input)
%Function that creates a stalagmite function with base equations

x=input(1);
y=input(2);

%the equation for a stalagmite function depends upon two input variables, x
%and y. The output, "f" is calculated with the following formulae:
% f1=(sin(5.1*pi*x+0.5))^6;
% f2=(sin(5.1*pi*y+0.5))^6;
% f3=exp(-4*log(2)*((x-0.0667)^2)/0.64);
% f4=exp(-4*log(2)*((y-0.0667)^2)/0.64);

f1=(sin(5.1*pi*x+0.5))^6;
f2=(sin(5.1*pi*y+0.5))^6;
f3=exp(-4*log(2)*((x-0.0667)^2)/0.64);
f4=exp(-4*log(2)*((y-0.0667)^2)/0.64);

F = (f1*f2*f3*f4); This is the stalagmite function

% F will give a Stalgtite (i.e., un upward
% mounting peaks function)
% multiplying the product with (-1) will invert the values towards the
% negative side. This indeed is the Stalagtite function with peaks
% (troughs) facing downwards. This is required in order to find the minimum
% of the function as the ga (Genetic Algorithm) function  of MatLab only gives us the
% minimum values.



output = - F; % returns the output values for each x and y input to the function handle. This is the inverted stalagmite or stlagtite function. This is needed for the GA function which only find the minimum.

end

In the above function, the returned value is not the stalagmite function but, the inverted stalagmite (stalactite) function. This is needed for the GA to work. As the MatLab GA only finds minimum values.

2. Code for making calling of the Modified Stalagmite function:

function [t] = Stalagmite(x,y)
 
%A function to call the modified stalagmite function
    %This function takes in input values x and y and creates a meshgrid.
    %Each value of the meshgrid is fed to the "stagfun" function, which
    %returns the values and stores in the "g" matrix. 

[xx, yy] = meshgrid(x,y); %creates a meshgrid of input values x and y.

n = length(x); % defines the length of x,y, xx or yy (all must be of equal length).

g = zeros(n,n); 
%preallocating the size of matrix g in order to increase performance speed
 
for i = 1:n 
    for j=1: n
        input(1) = xx(i,j); % isolates each element of xx meshgrid with (i,j) index
        input(2) = yy(i,j); % isolates each element of xx meshgrid with (i,j) index

    %now each isolated value, i.e., xx(i,j) and yy(i,j) are stored in
    %the matrix "input"

        g(i,j)= stagfun(input); % creates a matrix g, of order n,
  % with each element calculated as the stalagmite function output point

    end
end

t = g; % t is the output from this function that will returned whenever this function is called.

end

  
  

3. Code For Determining The Global Maxima Using The Genetic Algorithm Function:

number of iterations = 100;

population size = 200;

%Determiinig the peak or global maxima of the stalagmite function

clear all  % clears all workspace variables 
close all  % closes all plot and figures
clc        % clears the command window screen

x=linspace(0,0.6,150);   % linespace for the X-axis, i.e., x values range
y=linspace(0,0.6,150);   % linespace for the Y-axis, i.e., y values range

%the linespace of x and y are not the search space for the ga. To set that, we 
%need to use the syntax inputs of the ga function.

f = Stalagmite(x,y);     % calling the Stalagmite function and storing the values in matrix "f"

num =100;  % number of iterations to be done. 
     
     %If we use ga only one time (or num = 1), we may not get the desired value
     %that we look for. Each iteration may give different results from the ga
     %and thus increasing the number of iterations will increase the changes of
     %the occurrence (converge to the actual) of the best solution.


Fopt= zeros(1,num);  

% The purpose of preallocating matrices with zeros (or ones etc.)
% is to create their size beforehand. When the
% iterations are done, MatLab only re-writes the values. This process takes
%less time than re-creating the sizes each time of the iteration.

xopt= zeros(1,num);
yopt= zeros(1,num);


options = optimoptions('ga');  
                              %this handle is used to input 'options' to the 
                              % ga function. First, we call on the options
                              % and provide the optimisation technique,
                              % here, the genetic algorithm (ga)
                              
options = reduce(options,'PopulationSize',200);

% providing the population size. More the population size, more number
%of data points will be checked against the ga, and better will be the
%results. However, with large population size, the algorithm will consume
%more time for the processing 
%notice that there are 150*150 data points in the f matrix. 

tic  %starts time counting. used for calculating the time it takes to
     % execute specific part of the script
       
for i = 1:num
    [inputs, Fopt(i)] = ga(@stagfun,2,[],[],[],[],[0;0],[0.6;0.6],[],options);

    % 'inputs' stores the x and y values that the ga algorithm returns
    % Fopt stores the optimum value that each iteration of the ga returns.
    % This stores the minimum z (function) values
    % The [] square brakets are necessary to provide in the function
    % handle even if they are empty. they are for linear equality and non 
    % equality constrains. the [0:0] and [0.6,0.6] are the lower and upper
    %bounds of the function. If this is not properly set, values beyond the
    %function domain will be returned by ga. In other words, these provide the
    % linespace for the ga function. The other two [] empty brackets are
    % for non linear inequality and interger inequality
    
    xopt(i) = inputs(1); %Stores the values of x returned by ga at each iteration 
    yopt(i) = inputs(2); %Stores the values of y returned by ga at each iteration 
    
  clc % clears the message "Optimization terminated: average change in the
      % fitness value less than options.FunctionTolerance."

end


Time = toc % stops time counting and provides the elapsed time value

figure(1)   % name of the figure

subplot(2,1,1);
%subplot command divides the current figure into a 2-by-1 grid and creates
% axes for a subplot in the position 1. Here, we will plot the modified
% Stalagmite function

hold on   %holds on the previous figure

surf(x,y,f); %plotting the surface of modified stalagmite function. 
              

shading interp; %creates a shading pattern instead of grid mesh. For better
                % visualisation.
                
xlabel('x-Values'); % labels the X-axis
ylabel('y-Values'); % labels the Y-axis
zlabel('Z/Function Values') %labels the Z-axis

plot3(xopt,yopt,Fopt,'marker','o','markersize',6,'markerfacecolor','r')
               % plot3 plots a line in 3-space through the points xopt, yopt and Fopt
               % plotting the point on the same figure with xopt, yopt and
               % zopt coordinates with a red marker i.e., the minimum value returned by 
               % the ga optimisation.
              
               
title('Modified Stalagmite Function') % a title for the figure
legend('Function','Minima') %provides visualisation on the figure

subplot(2,1,2) % creates the subplot at location 2.

plot(Fopt,'r') %plotting the Fopt value returned bu ga for each iteration

xlabel('Iterations') 
ylabel('Function Minimum')
title('Minimum Value At Different Iterations')
%It is seen that with the iterations, the GA predicts the same minimum value
% most of the times. Only at few iterations, does GA predict some different
% results. However, when the iteration are increased, we see that the most
% predicted value falls at (x,y) ~ (0.0668,0.0668) and the function's value
% is (-1)

%By looking at the subplot(2,1,2) [the iteration vs Fopt value] we can see
%that for most of the iterations the GA predicts a minima of (-)1. Only for
%few iterations does the GA predict some higher (less negative) values.
%This helps in deciding that the maxima must be at that point where most 
%iterations of the GA predicts.

%Now, to determine the maxima of the function, the xopt and yopt
%will remain the same. However, the value of Fopt needs to be put in
%positive value. Thus, multiplying Fopt values with (-1), will give the maximas,

Maximum = - (min(Fopt)); % to extract a single minimum value from the Fopt
                         % and to convert it into the maximum value


%To display the x and y  values at which maxima occurs. (Here, the inputs
%are of the last iteration. 
M = ['The global maxima of the said Stalagmite funtion is: z = '  , num2str(Maximum),' at (x,y) = ', num2str(inputs)];
disp(M)

N = ['The study took ',num2str(Time),' Seconds']
disp(N)

The function plot and the variation of Fopt against iterations plot:

Viewing the plot in 3D

We see from the above results that the Fopt returned by GA is always (-1) except for one iteration where the value is about (-0.85). Clearly, for most iterations (99% in this run), -1 is the global minima of the modified stalagmite function. Hence, it is obvious that we will select the most predicted value as the global minima.

Now for the given condition, we need to find the global maxima and not the minimum. To do this, simply multiply the minima with (-1), and the global maxima of the Stalagmite function is given. Also, the values for x and y (xopt, yopt) for the global maxima will remain the same.

In my study, the global maximum is at (x,y) = (0.0668, 0.0668). The value is 1.

% this code is to be run with the previous main code. It will plot the Stalagmite function and the global maxima.

Figure(2)
surf(x,y,-f)
shading interp
hold on;
plot3(inputs(1),inputs(2),1,'marker','o','markersize',5,'markerfacecolor','r') 

The stalagmite function and the maxima: