Optimization of stalagmite function

Aim

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

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.

Notion of Natural Selection

The process of natural selection starts with the selection of fittest individuals from a population. They produce offspring which inherit the characteristics of the parents and will be added to the next generation. If parents have better fitness, their offspring will be better than parents and have a better chance at surviving. This process keeps on iterating and at the end, a generation with the fittest individuals will be found.

This notion can be applied for a search problem. We consider a set of solutions for a problem and select the set of best ones out of them.

Five phases are considered in a genetic algorithm.

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

Initial Population

The process begins with a set of individuals which is called a Population. Each individual is a solution to the problem you want to solve.

An individual is characterized by a set of parameters (variables) known as Genes. Genes are joined into a string to form a Chromosome (solution).

In a genetic algorithm, the set of genes of an individual is represented using a string, in terms of an alphabet. Usually, binary values are used (string of 1s and 0s). We say that we encode the genes in a chromosome.

Fitness Function

The fitness function determines how fit an individual is (the ability of an individual to compete with other individuals). It gives a fitness score to each individual. The probability that an individual will be selected for reproduction is based on its fitness score.

Selection

The idea of selection phase is to select the fittest individuals and let them pass their genes to the next generation.

Two pairs of individuals (parents) are selected based on their fitness scores. Individuals with high fitness have more chance to be selected for reproduction.

Crossover

Crossover is the most significant phase in a genetic algorithm. For each pair of parents to be mated, a crossover point is chosen at random from within the genes.

For , consider the crossover point to be as shown below.

Offspring are created by exchanging the genes of parents among themselves until the crossover point is reached.

The new offspring are added to the population.

Mutation

In certain new offspring formed, some of their genes can be subjected to a mutation with a low random probability. This implies that some of the bits in the bit string can be flipped.

Termination

The algorithm terminates if the population has converged (does not produce offspring which are significantly different from the previous generation). Then it is said that the genetic algorithm has provided a set of solutions to our problem.

Comments

The population has a fixed size. As new generations are formed, individuals with least fitness die, providing space for new offspring.

The sequence of phases is repeated to produce individuals in each new generation which are better than the previous generation.

Syntax of GA for MATLAB

1. x= ga(fun,nvars) finds a local unconstrained minimum, x, to the objective function, fun. nvars is the dimension (number of design variables) of fun.

Note

Passing Extra Parameters explains how to pass extra parameters to the objective function and nonlinear constraint functions, if necessary.

2. x= ga(fun,nvars,A,b) 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) 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) 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) 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) minimizes with the default optimization parameters replaced by values in options. (Set nonlcon=[] if no nonlinear constraints exist.) Create options using optimoptions.
7. x= ga(fun,nvars,A,b,Aeq,beq,lb,ub,nonlcon,intcon) or x = ga(fun,nvars,A,b,Aeq,beq,lb,ub,nonlcon,intcon,options) requires that the variables listed in intcon take integer values.

Note

When there are integer constraints, ga does not accept nonlinear equality constraints, only nonlinear inequality constraints.

8. x= ga(problem) finds the minimum for problem, a structure described in problem.
9. [x,fval] = ga(___), for any previous input arguments, also returns fval, the value of the fitness function at x.

• [x,fval,exitflag,output] = ga(___) 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.
• [x,fval,exitflag,output,population,scores] = ga(___) also returns a matrix population, whose rows are the final population, and a vector scores, the scores of the final population.

Stalagmite function

MATLAB code for stalagmite function

%stalagmite function

function s = stalagmite(input_vector)       %definition of objective function
x = input_vector(1);
y = input_vector(2);
f1x = (sin((5.1*pi*x)+0.5))^6;
f1y = (sin((5.1*pi*y)+0.5))^6;
f2x = exp((-4*log(2))*(((x-0.0667)^2)/0.64));
f2y = exp((-4*log(2))*((y-0.0667)^2)/0.64);
s = -(f1x*f2x*f1y*f2y);                     %the minus sign gives the maximun value
end

• the input of the function is a row vector containing two columns
• a minus sign is used as GA evaluates the minimunm value and we are interested in the maximum value of the function

MATLAB code for GA

%genetic algorithm

clear all;
close all;
clc;

%input variables

x = linspace(0,0.6,150);
y = linspace(0,0.6,150);
population = input('Enter population size = ');


%2d array
%evalauting inputs for stalagmite functiom
[xx,yy] = meshgrid(x,y);                        %creating 2D grid from 1D array
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);         %calling stalgmite function and storing in f
    end
end



%study 1 - Unbounded stalagmite function
%calculates the study time 
tic
for i=1:population
[results,fopt1(i)]=ga(@stalagmite,2);      %calculating the minimum value of the function
xopt1(i)=results(1);                       %x coordinate of minimum
yopt1(i)=results(2);                       %y coordinate of minimum
end
study1_time = toc                          


% plot
figure(1)
subplot(2,1,1)
surfc(xx,yy,-f);        %plotting stalagmite function
shading interp
hold on

%plotting ga values
plot3(xopt1,yopt1,-fopt1,'-o','MarkerFaceColor','r');  
title('Stalagmite function - Unbonded');
xlabel('X axis');
ylabel('Y axis');
grid on
subplot(2,1,2)
plot(-fopt1);
xlabel('Iterations')
ylabel('Optimum value')


%study 2 - stalagmite function - bounded
%calculates the study time 
tic
for i=1:population
[results,fopt2(i)]=ga(@stalagmite,2,[],[],[],[],[0,0],[0.6,0.6]);      %calculating the optimized value of the function
xopt2(i)=results(1); %x coordinate of minimum
yopt2(i)=results(2);%y coordinate of minimum
end
study2_time = toc 


% plot
figure(2)
subplot(2,1,1)
surfc(xx,yy,-f);%plotting stalagmite function
title('Stalagmite function - Bounded')
hold on
shading interp

%plotting ga values
plot3(xopt2,yopt2,-fopt2,'-o','MarkerFaceColor','r');
xlabel('X axis');
ylabel('Y axis');
grid on
subplot(2,1,2)
plot(-fopt2);
xlabel('Iterations')
ylabel('Optimum Value')



%study 3 - increasing polpulation size
%calculates the study time 
tic


options = optimoptions('ga','PopulationSize',300);
for i=1:population
[results,fopt3(i)]=ga(@stalagmite,2,[],[],[],[],[0,0],[0.6,0.6],[],options);      %calculating the optimized value of the function
xopt3(i)=results(1); %x coordinate of minima
yopt3(i)=results(2);%y coordinate of minima
end
study3_time = toc 


% plot
figure(3)
subplot(2,1,1)
surfc(xx,yy,-f);%plotting stalagmite function
title('Stalagmite Function - increased population size')
hold on
shading interp

%plotting ga results
plot3(xopt3,yopt3,-fopt3,'-o','MarkerFaceColor','r');
xlabel('X axis');
ylabel('Y axis');
grid on
subplot(2,1,2)
plot(-fopt3);
xlabel('Iterations')
ylabel('Optimum Value')


%tabular values

X_coordiantes = [xopt1(1,50);xopt2(1,50);xopt3(1,50)];
Y_coordinates = [yopt1(1,50);yopt2(1,50);yopt3(1,50)];
Maxima = [-fopt1(1,50);-fopt2(1,50);-fopt3(1,50)];
Time = [study1_time;study2_time;study3_time];
row_names = {'Unbounded','Bounded','Population increase',};
t = table(X_coordiantes,Y_coordinates,Maxima,Time,'RowNames',row_names)

• user input is required to enter the size of the population
• meshgrid is used to create a 2D array for input of stalagmite function
• using for loop the input vector for the function is created
• the no of loop iterations are set to the size of the population defined at the start
• STUDY1 – Unbounded
  • The algorithm is predicting values in an unbounded search space.
  • This can be visualized through the plot3 command
  • The time required for executing the code is also calculated
• STUDY2 – Bounded
  • The search space is narrowed by providing the upper and lower bounds
• STUDY3 – Increasing population size
  • By increasing the population size the no of guesses at the initiation stage are increased which generates more mutations because of which the algorithm accuracy to predict the optimum increases.

Observation

Study 1

Study2

Study3

Table

The maxima is at (0.066832,0.066833) and the value is 1.