Using Genetic Algorithm for Optimising Stalagmite Function

Aim:- To implement Genetic Algorithm for finding global maxima of stalagmite function for a given range of input, and perform a statistical study of the ga() function of MATLAB for 50 iterations.

Theory:-

What is Stalagmite function?

Stalagmite function is a 3-dimensional mathematical function which takes two inputs x and y and produce output along the z-axis. It is given as:-

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

`f_(1,x)= sin^6⁡(5.1πx+0.5)` `f_(1,y)= sin⁡^6(5.1πx+0.5)`

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

The function gets the name due to its resemblance to the stalactites and stalagmites found in the caves.

What is Genetic Algorithm?

Genetic Algorithm is a search algorithm which implements the process of natural selection where the fittest individuals (variables with their values) are selected for reproduction in order to produce offspring of the next generation.

It is based on Charles Darwin’s theory of natural evolution.

A 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.

Classical Algorithm	Genetic Algorithm
Generates a single point at each iteration. The sequence of points approaches an optimal solution.	Generates a population of points at each iteration. The best point in the population approaches an optimal solution.
Selects the next point in the sequence by a deterministic computation.	Selects the next population by computation which uses random number generators.

Five phases are considered in a genetic algorithm.

1. Initial population
2. Fitness function
3. Selection
4. Crossover
5. Mutation

Syntax of genetic algorithm function in MATLAB.

[x , fval , exitflag , output , population , scores] = ga(fun , nvar , A,b,Aeq,beq,lb,ub,nonlcon,Intcon,options)

First, let’s talk about the input parameters of function ga().

1. fun – It is the function handle of the function to be optimised. It’s given as ‘@’+function name. In our case it is @Stalagmites
2. nvar – This integer represents the number of input required by the function passed. Stalagmites function takes two input as discussed earlier.
3. A , b – These parameters are used to specify the Linear inequalities among inputs.
   1. A*x <= b
   2. A is M-by-nvars matrix, where M is number of inequalities
   3. b is 1-by-M column matrix
4. Aeq , beq – Similar to previous case, Aeq and beq are used to specify Linear Equalities
   1. Aeq*x <= beq
   2. Aeq is M-by-nvars matrix, where M is number of inequalities
   3. beq is 1-by-M column matrix
5. lb , ub – These row or column vectors specify lower and upper bounds of inputs respectively.
6. nonlcon – This again is a function handle which is used to mention non-linear (such as quadratic, cubic, exponential,etc) constraints.
7. options – This , a sort of vector, consists of other parameters of genetic algorithm such as initial PopulationSize, CrossoverFraction, MaxGenerations etc. These parameters are enclosed within braces “{…}”.

Coming to the Output side:-

1. x – It is a 1-by-nvar row matrix consisting of values of inputs for which ga() calculated optimum value of passed function.
2. fval – It is the optimum value of the function. Generally fval = fun(x).
3. exitflag – exit lag is an integer, which states the reason for which ga() has stopped. A list of integers and their meaning can be found in the documentation of ga()
4. output - Information about the optimization process, returned as a structure with different fields such as problemtype, rngstate, etc.
5. population – This population-by-nvars matrix contains the final population.

Code:-

For the purpose of studying various parameters of the genetic algorithm function of MATLAB, the code is divided into three parts. All these three cases uses the same set of inputs and function output for testing genetic algorithm.

1. Part –I : Running ga() function on stalagmite function for 50 times (num_cases) and observing time required and optimum output and its corresponding inputs at each iteration.

   1. % Programme for statistical behaviour of ga function
      clear all;
      close all;
      clc;
      
      % Creating input vectors
      x = linspace(0,0.6,150);
      y = linspace(0,0.6,150);
      % Calling stalgmite function to get 2-D Output in f
      f = stalagmite(x,y);
      num_cases = 50;     % Defining number of iterations
      
      tic;
      for i=1:num_cases
          [inp , fopt(i)] = ga(@Stalagmites,2);   % Calling Genetic Algorithm function
          xopt(i) = inp(1);                       % Splitting inp vector into x and y coordinate
          yopt(i) = inp(2);
      end
      % Recording time taken by ga for 50 iterations
      study1_time = toc;     
      figure(1);
      subplot(2,1,1);
      hold on;
      % Plot - Surface Plot for Stalagmite function
      surfc(x,y,f);
      % Line plot for fopt, xopt and yopt
      plot3(xopt,yopt,fopt,'marker','o','markersize',5,'markerfacecolor','r');
      grid on;
      shading interp;
      title('Unbounded Inputs');
      xlabel('X-Axis');
      ylabel('Y-Axis');
      zlabel('Z-Axis Stalagmite Function');
      
      subplot(2,1,2);
      plot(fopt);
      xlabel('Iteration');
      ylabel('Function Minimum');
      grid on;
      
      % Calling User-Defined function for getting useful data
      useful_ga_data(fopt,xopt,yopt,study1_time);

   2. The codes starts with basic three commands and setting up workspace by clearing variables to avoid clashes, if any and closing previous graphs.
   3. In next two lines we use linspace() function to specify ranges for both inputs , x and y from 0 to 0.6 with 150 points in between.
   4. After this user-defined function stalagmite() function takes these inputs and returns a 2 – dimensional array consisting of corresponding values of stalagmite function in the given range to f variable.
   5. Then we declare number of iterations as 50 into num_cases
   6. Now we run for loop for num_cases times to implement genetic algorithm.
   7. In this our first study we just mention bare minimum info to the function, that is the function and number of variables.
   8. The optimum (minimum) value at each iteration and its corresponding input is stored in a row vector into inp and fopt The inp variable is further broken down into xopt and yopt variables.
   9. MATLAB’s tic-toc command is used to measure and record time taken by ga() function to finish these 50 iterations into study1_time
   10. subplot() function divides the figure 1 into 2 rows.
   11. The first row contains the inverted stalagmite plot created by in-built surf() function by passing x, y and corresponding stalagmite output value f.
   12. The 3D plot of different optimum values as predicted by the ga() function over 50 iterations. The 3D plot of optimum value and corresponding inputs is achieved by plot3()
   13. The second row contains a simple 2D plot depicting optimum value over different iterations. In-built plot() function plots array element along y-axis against their index along x-axis.
   14. Once the study is over, now we display useful data onto console by calling user-defined function named usefuk_ga_data() function and passing variables fopt, xopt, yopt and study1_time as parameters.
2. Part –II– In the second study, in addition to bare minimum information, we also specify the upper and lower limits of inputs x and y as per syntax mentioned earlier. Lower and Upper limit of both x and y are set as 0 and 1 respectively.

   1. tic;
      for i=1:num_cases
          [inp , fopt(i)] = ga(@Stalagmites,2,[],[],[],[],[0,0],[1,1]);
          xopt(i) = inp(1);
          yopt(i) = inp(2);
      end
      study2_time = toc;
       
      figure(2);
      subplot(2,1,1);
      hold on;
      surfc(x,y,f);
      plot3(xopt,yopt,fopt,'marker','o','markersize',5,'markerfacecolor','r');
      grid on;
      shading interp;
      title('Bounded Inputs');
      xlabel('X-Axis');
      ylabel('Y-Axis');
      zlabel('Z-Axis Stalagmite Function');
      zlabel('Z-Stalagmite  Function');
       
      subplot(2,1,2);
      plot(fopt);
      xlabel('Iteration');
      ylabel('Fuunction Minimum');
      grid on;

   2. This section also uses the same input and function output.
   3. The same tic-toc command is again used to monitor CPU run time and store recorded value into study2_time
   4. Here also ga() is running for 50 iterations, but with specified limits, restricting it to take both inputs x and y between 0 and 1. These are specified in a row-vector format as [0,0] and [1,1] as the lower and upper limit respectively.
   5. Similar to the first case optimal value and corresponding inputs are stored into inp and fopt variables.
   6. These are again plotted into respectively 3-D and 2-D (fopt vs iteration) in row fashion. Like previous case useful data is also printed onto the console to facilitate the user.
3. Part - III - In this section we specify initial population

   1. % Study-3 - Including Upper and Lower Bound along with Initial Population
      options = optimoptions('ga');
      options = optimoptions(options,'PopulationSize',170);
      tic;
      for i=1:num_cases
          [inp , fopt(i)] = ga(@Stalagmites,2,[],[],[],[],[0,0],[1,1],[],[],options);
          xopt(i) = inp(1);
          yopt(i) = inp(2);
      end
      study3_time = toc;
       
      figure(3);
      subplot(2,1,1);
      hold on;
      surfc(x,y,f);
      plot3(xopt,yopt,fopt,'marker','o','markersize',5,'markerfacecolor','r');
      grid on;
      shading interp;
      title('Bounded Inputs');
      xlabel('X-Axis');
      ylabel('Y-Axis');
      zlabel('Z-Axis Stalagmite Function');
      zlabel('Z-Stalagmite  Function');
       
      subplot(2,1,2);
      plot(fopt);
      xlabel('Iteration');
      ylabel('Fuunction Minimum');
      grid on;

   2. In the first line of this section we create an option object for ga function, a part of optimisation toolbox by using optimoptions() function by mentioning ‘ga’ as a parameter. This create option instance specially for ga.
   3. In the next line we specify 170 as initial population by using ‘PopulationSize’ By default initial ‘PopulationSize’ is set to only 50.
   4. This option instance is then passed to ga()function for each iteration as per the syntax.
   5. Similar to the previous two cases, time is recorded into study3_time variable. Since initial population is increased, CPU takes more time.
   6. The two graphs 3-D graph showing stalagmite function, fopt and inp and 2-D graph are also plotted.
4. stalagmite() Function - This user-defined function is returns 2-D array “ f ” of outputs for a 3-D mathematical function.
   

   1. function [f] = stalagmite(x,y)
          [xx,yy] = meshgrid(x,y);
          for (i=1:length(xx))
              for (j=1:length(yy))
                  x_c = xx(i,j);     % x-coordinate
                  y_c = yy(i,j);     % y-coordinate
                  f(i,j) = Stalagmites([x_c,y_c]);
              end
          end
      end

   2. It takes two row-vectors as input or parameter.
   3. Then it creates 2-D array of x and y using meshgrid() function and returns it to xx and yy variables respectively.
   4. xx is a 2-D array, whose rows are a repetition of the specified range. Similarly, yy is a 2-D array whose columns are repeated with the specified range.
   5. Now we use two nested loops, the outer one for taking index of x and the inner one for taking index of y.
   6. These two indices are then used to extract x and y These two coordinates are further passed into Stalagmites() function to get respective function value, which is then stored in “ f “ (another 2-D array).
5. Stalagmites() Function - This is the main stalagmite function, which return single output for given pair of x and y inputs.

   1. function [z] = Stalagmites(inp_vect)
          x = inp_vect(1);
          y = inp_vect(2);
          f1x = (sin(5.1*pi*x+0.5))^06;
          f1y = (sin(5.1*pi*y+0.5))^06;
          
          f2x = exp(-4*log(2)*((x-0.0667)^2)/0.64);
          f2y = exp(-4*log(2)*((y-0.0667)^2)/0.64);
          
          z = -1*f1x*f1y*f2x*f2y;
      end

   2. The function takes a row vector inp_vect as input parameter. The first two elements of this vector, inp_vect(1) and inp_vect(2) are assigned as x and y
   3. After this, four terms f1x, f1y, f2x and f2y are calculated according to the formula discussed in the theory.
   4. These four terms are then multiplied to get the final answer “ z “. All these values are positive. We want to find optimum, that is the height of the tallest peak of the stalagmite function.
   5. Since genetic algorithm finds the minimum of a function, therefore we multiply final output/answer by -1 to get inverted stalagmite, and thus find most negative (minimum), output.
6. useful_ga_data() Function - Function to print useful data from statistical study of ga().

   1. % Function to display usegul data for statistical study of Genetic Algorithm over Stalagmite function.
      function [z] = useful_ga_data(fopt,xopt,yopt,time)
         clc;
         % Extracting out Useful Data from study
         f_optimum = min(fopt);          % Finding minimum output (most optimum function output)
         pos = find(fopt==min(fopt));    % Finding index/position of minimum output (most optimum function output)
         x_min = xopt(pos);              % Finding index of x coordinate corresponding to most minimum fopt
         y_min = yopt(pos);              % Finding index of y coordinate corresponding to most minimum fopt
         
         % Displaying extracted information with proper formatting of floating point number
         fprintf("\nOptimum Value's Position                     = %d" , pos);   
         fprintf("\nMost Optimum Function Value in 50 Iterations = %4.3f" , f_optimum);
         fprintf("\n(x - y) Coordinate at Optimum Values         = (%4.3f , %4.3f)" , x_min , y_min);
         fprintf("\nTime taken for 50 Iterations                 = %4.3f sec\n" , time);
      end
      
      

   2. This function is a non-returning type function, that is doesn’t give any value back to the main programme, from where it is being called. It just prints the formatted information onto the command window.
   3. f_optimum – This variable stores the most optimum function output produced in 50 iterations. We use min() function on fopt array for this purpose.
   4. pos – find() function is used to get index of f_optimum from fopt array and stores it in pos
   5. x_min and y_min – The above position is then used to extract x and y coordinates corresponding to f_optimum from xopt and yopt
   6. These information; index, most optimum output, its corresponding x-y inputs along with the time required by ga() for 50 iterations is then displayed onto console using fprintf() function with proper formatting of floating/double data types.
   7. Format specifier %4.3f tells MATLAB to print double data type with four fields, 3 of which should be after the decimal point

Console Outputs and Plots:-

Case - I Simple Run

Case - II With Upper and Lower Limits

Case - III With Limits and Increased Initial Population

Observation:-

We observe that ga() predicts the correct set of answer (optimum value and corresponding inputs) in the second and third case. It can be seen that time required for each case increases as we add more parameters and to get accurate results. By taking only limits in the second case, the frequency of ga() deviating from optimum value is more as compared to the third case. Increasing initial population size in the third case reduces this ga()’s frequency (two times)of deviating from optimum but increases computation time.

Result and Conclusion:-

For getting accurate results, one should give proper instruction of limits, initial conditions, even if it increases computation time.

References:-

• Mathworks Documentation
• Youtube - MATLAB
• Using Genetic Algorithm in MATLAB
• Genetic Algorithm for Travelling Salesman Problem
• Toweardsdatascience
• Genetic Algorithm - Wikipedia