Optimization of a Stalagmite Function using Genetic Algorithm in MATLAB

GENETIC ALGORITHM:

Introduction:

Genetic Algorithm is based upon Charles Darwin's principle of Natural selection which explains the natural law of survival of the fittest. It is based on a heurisitic approch which means, it learns and takes decision itself during the processing time span.

Application:

It is generally used in engineering, statistics and economics for optimising a function to its minimum or near-minimum. Genetic Algorithm (GA) is a subset of much larger concept of Evolutionary Computation.

Detailed concept:

In GA, there is an initial population randomly selected/ defined by the algorithm itself, just as in nature. It then calculates the fitness of the genes of population and ranks the fittest genes in decreasing order. Further, just as in nature, it allocates a pool of fittest population and then cross overs the parents with each other. A cross over point is chosen within each parent and the respective genes are crossed over. It does not discard the weaker population entirely so that it does not lose the good genes from it. After the cross over, a new set of population is created and is called 1st generation. It then repeats the above process of calcualating the fitness, creating mating pool, crossing over the genes and creating new generation until the fitness value remains unchanged. This proves we have reached towards the best population.

Correlating with a minimisation problem, the ga first defines an initial set of solutions. It calculates the value of fitness function using the population. The fitness ranking is based upon the minimum value of fitness function obtained. The new population is then selected near these coordinates and again fitness function is evaluated. The process is repeated until the minimum value is found.

Going into further detail,

Population selection strategies:

1. Elite: The parents with highest fitness are selected. The king in a particular population is left untouched in crossover and is sent to next population. If a higher fitness parent is found in new population, it is then crossed over with this king.

2. Roulette: The parents are arranged so as to occupy the area of a circle i.e. a roulette wheel. The area occupied by a parent is directly proportional to its fitness score. Selection is done from this wheel and hence the parent with higher fitness has a higher chance to get selected. Duplication is avoided.

3. Rank: Similar to roulette but the proportion of area is proportional to the rank based on ascending order of fitness scores. This allows more diversity and avoids local convergence.

4. Tournament: Random population is selected for mating and tournaments (equal to required size of mating population) are conducted between them. The parent with higher fitness wins the particular tournament. 'Selection pressure' determines the randomness of selection by selecting the fittest individuals first and is directly proportional to the time for convergence. However, high pressure may be lead to less diversity and may result in local minima rather than global minima.

Mutation: To maintain diversity, the genes of parents (chromosomes) are sometimes mutated within itself and new children are formed.

GA solves 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. However this process of genesis leads to huge computation time. In order to reduce this, techniques of hybrid computation is used wherein, GA is used to find the approximate location of global minima and traditional numerical method technique is used to converge upon the solution in the located area since it consumes very less time.

PARENT CODE:

clear all
close all
clc

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

%Creating 2D mesh:
[xx yy]=meshgrid(x,y);

%Evaluating 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);
       val(i,j)=stalagmite(input_vector);
    
   end
   
end

%Plotting the stalagmite function
figure('Name','Stalagmite function','NumberTitle','off')
    surfc(xx,yy,-val);
    shading interp;

%Running the genetic algorithm
[coord, htval]=ga(@stalagmite,2);


%Case 1: Statistical behaviour

%Defining no. of iterations:
num_cases=50;

tic %Setting up timer to know the time required (tic-toc command)

for k=1:num_cases %Running the genetic algorithm multiple times

    [inputs, fopt(k)]= ga(@stalagmite,2);
    xopt(k)=inputs(1);
    yopt(k)=inputs(2);
     
end

fopt=fopt*(-1); %Since we want the maxima of stalagmite function

case_1_time=toc

%Creating new figure window and plotting the data
figure('Name','Case 1: Statistical Behaviour','NumberTitle','off') 

subplot(2,1,1);
    hold on
    srf_pt= surfc(x,y,-val);
    shading interp;
    fopt_plot(3)= plot3(xopt,yopt,fopt,'marker','o','markersize',5,'markerfacecolor','r');
    title('Unbounded inputs');

subplot(2,1,2);
    twodplot(4)= plot(fopt);
    axis([0 50 -1 2 -1 1]);
    xlabel('Iterations');
    ylabel('Function minimum');


%Case 2: Statistical behaviour with upper and lower bounds

tic
for l=1:num_cases
    
    %The seventh and eighth position of 'ga' syntax is for setting limits
    [inputs, fopt_2(l)]= ga(@stalagmite,2,[],[],[],[],[0;0],[1;1]);
    xopt_2(l)=inputs(1);
    yopt_2(l)=inputs(2);
     
end
case_2_time=toc

fopt_2=fopt_2*(-1);

figure('Name','Case 2: Bounded upper and lower limits of search space','NumberTitle','off')

    subplot(2,1,1);
        hold on
        srf_pt=surfc(x,y,-val);
        shading interp;
        fopt_plot(5)=plot3(xopt_2,yopt_2,fopt_2,'marker','o','markersize',5,'markerfacecolor','r');
        title('Bounded inputs');

    subplot(2,1,2)
        twodplot(6)= plot(fopt_2);
        axis([0 50 -1 2 -1 1]);
        xlabel('Iterations');
        ylabel('Function minimum');


%Case 3: Increasing the no. of GA iterations

%Setting population size to 200
%optimoptions allows to change the options in an optimization algorithm
options=optimoptions('ga','PopulationSize',200); 

tic

for m=1:num_cases
    
  [inputs, fopt_3(m)]= ga(@stalagmite,2,[],[],[],[],[0;0],[1;1],[],[],options);
  xopt_3(m)=inputs(1);
  yopt_3(m)=inputs(2); 
        
end

case_3_time=toc

fopt_3=fopt_3*(-1);

figure('Name','Case 3: Increasing the Population Size','NumberTitle','off')

    subplot(2,1,1);
        hold on
        srf_pt= surfc(x,y,-val);
        shading interp;
        fopt_plot(7)= plot3(xopt_3,yopt_3,fopt_3,'marker','o','markersize',5,'markerfacecolor','r');
        title('Bounded inputs with high population size');

    subplot(2,1,2);
        twodplot(8)= plot(fopt_3);
        axis([0 50 -1 2 -1 1]);
        xlabel('Iterations');
        ylabel('Function minimum');


%Subplots for comparison

figure('Name','Comparison plots','NumberTitle','off')

    subplot(4,2,[1,2]);
        surfc(xx,yy,-val);
        shading interp;

ctr=1;

for si=3:2:8

      
    sbp(si)= subplot(4,2,si);
    hold on
    copyobj(srf_pt,sbp(si));
    shading interp;
    copyobj(fopt_plot(si),sbp(si));
        
        if ctr==1
            title('Unbounded inputs');
        end
        
        if ctr==2
            title('Bounded inputs');
        end
        
        if ctr==3
            title('Bounded inputs with high population size');
        end
      ctr=ctr+1;  
        
end

    
      
for pi=4:2:8

    sbp(pi)= subplot(4,2,pi);
    copyobj(twodplot(pi),sbp(pi));
    axis([0 50 -1 2 -1 1]);
    xlabel('Iterations');
    ylabel('Function minimum');
    
end

CODE for stalagmite.m :

function output=stalagmite(input_vector)

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

val1 = (sin(5.1*pi*x + 0.5))^6;

val2 = exp(-4*log(2)*((x-0.0667)^2)/0.64);

val3 = (sin(5.1*pi*y + 0.5))^6;

val4 = exp(-4*log(2)*((y-0.0667)^2)/0.64);

output = val1*val2*val3*val4;

output=output.*(-1);

end

Explanation of code:

1. Parent code

The code begins with clearing previous memory and command window.

Then it defines the coordinates wherein the fitness function is defined. These coordinate values are then used to create a meshgrid. Further, two 'for' loops are nested within each other to provide input vector for 'stalagmite' function and simultaneously provide the output values. A new figure window with the name 'Stalagmite function' is created and the surface plot for the function is plotted with interpolated shading so that the plot color pattern looks smooth. It is to be noted that the 'z' values are plotted in negative so that it resembles the original stalagmite function. Clarity will be obtained during the code explanantion of 'stalagmite.m'

The genetic algorithm is then made to run to evaluate the optimum position and values for the fitness function through 'coord' and 'htval' respectively.

ga syntax:

[X, FVAL] = ga(FITNESSFCN,NVARS,A,b,Aeq,beq,lb,ub,NONLCON,options)

where,

X : coordinates of optimization solution

FVAL : value of optimization solution

FITNESSFCN : fitness function handle (here '@stalagmite')

NVARS : No. of variables (here 2)

A,b,Aeq,beq : used when linear inequalities are used

lb : defines lower bound of search space

ub : defines upper bound of search space

NONLCON : used for non linear inequalities

options : used to change parameters in the 'ga'

It is observed that coord and htval keeps on changing every time we run the code. This does not provide a reliable output and needs to be changed.

Case 1:

Here we define multiple iterations so that the ga runs multiple times. Our aim is to understand how ga algorithm searches for the optimum value.

The 'num_cases' variable defines the no. of iterations and the 'for' loop runs the ga for the same iterations. The value is stored in 'inputs'.

Then we create a new figure window with subplots. In the first subplot, we plot the surface for stalagmite function and overlap plot showing the optimum value of ga as per iterations. In the second subplot, we plot the optimum value of plot against the no. of iterations.

The tic-toc command displays the time required for operation to be carried out.

The result shows how the 'ga' algorithm searches at many different locations but does not provide the final, consistent result.

Case 2:

Here we define the upper and lower bounds for search space to allow 'ga' to concentrate on a smaller area. Plotting is done just as in case 1. The result shown is much better and we realise that the scatter is much lower and the solution is more concentrated near optimum solution. Very few values show deviation.

Case 3:

Here we have extracted the Population Size from 'ga' parameter using 'options' command and changed the value from 50(default) to 200. This method differs from case 2 in a way that since the orginal population is higher, higher no. of crossovers will take place and convergence towards optimum solution will be faster and more accurate right from first crossover. Plotting is to be followed just as in case 1.

For ease in representation and comparison, a new figure window is created and all the above subplots are placed in it for better comparison. Here two different for loops are used and the axes from first figure (alongwith graphcs) are called in the respective subplot using 'copyobj' command.

2. Code for stalagmite.m

The stalagmite function takes a row vector named 'input_vector' as an input and gives an output value named 'output'. The equations are defined in the function and is finally multiplied by '-1' to invert the function about X axis. This is done since 'ga' algorithm always minimises the fitness function, but we needed maximisation as the output. Minimum value searched by 'ga' in the inverted function will be the actual maximum value of original function.

Output:

Stalagmite function:

Case 1:

Case 2:

Case 3:

Comparison plots:

Output plot:

Using inbuilt optimization app for Genetic algorithm.

MATLAB provides an inbuilt optimization app with various solvers. Here we have selected 'ga-Genetic algorithm' solver.

1. The fitness function and no. of variables are given as input. All the remaining parameters are kept as default. It is observed that the result is inaccurate.

2. Now we set the upper and lower bounds to 0,0 and 1,1 and increase population size to 150. Better result than previous is obtained but the answer is still not accurate. Also the coordinate shown is 0.262,0.067.

3. Keeping the bounds same, we change the population size to 300. Accurate result is obtained and algorithm stops at approx 80 generations. Coordinate shown is 0.067, 0.067

4.Now we set the upper and lower bounds to 0,0 and 0.5,0.5 and change population size to 150. Accurate result is obtained and algorithm stops at around 75 generations

5. Keeping the bounds same, we change the population size to 300. Accurate result is obtained and algorithm stops at around 70 generations. This is the best case where optimization happens in the least no. of generations.

For a more complex problem, other parameters like Crossover fraction, Mutation function can be changed to obtain the best result in minimum generations.