GENETIC ALGORITHM TO OPTIMISE THE STALAGMITE FUNCTION

OBJECTIVE:

To optimise the stalagmite function and find the global maxima of the function.

INTRODUCTION:

Optimization is the process of making something better. It refers to finding the values of inputs in such a way that the 'best' output values are obtained. The definition of 'best' varies from problem to problem, but in mathematical terms, it refers to maximizing or minimizing one or more objective functions, by varying the input parameters.

GENETIC ALGORITHMS:

Nature has always been an inspiration to all of humanity. Genetic Algorithms (GAs) are search based algorithms based on the concepts of natural selection and genetics. GAs are a subset of a much larger branch of computation known as Evolutionary Computation.

In GAs, we have a pool or population of potential solutions to the problem at hand. These solutions then undergo recombination and mutation like in natural genetics, resulting in new children, and the process is repeated over different generations. Each individual or candidate solution is assigned a fitness value based on its objective function value and the best fit people have a greater chance of mating and producing more fit people. This is in line with the Darwinian Theory of “Survival of the Fittest”.

In this way, we continue to evolve better people or solutions over generations, until we achieve a stop criterion.

Genetic algorithms are random enough in nature, but they perform far better than random local search, as they exploit historical information.

Advantages of GAs

GAs have various advantages which have made them immensely popular. These include −

• Does not need any derived information that may not be available for many real-world problems.

• Is faster and more effective in comparison with traditional methods.

• Has very good parallel capabilities.

• Optimizes both continuous and discrete functions and also multi-objective problems.

• Provides a list of good solutions and not just a single solution.

• Always gets a response to the problem, which improves over time.

• Useful when the search space is very large and there are a large number of parameters involved.

Limitations of GAs:

Like any technique, GAs also suffer from a few limitations. These include −

• GAs are not suited for all problems, especially problems which are simple and for which derivative information is available.

• The fitness value is calculated repeatedly, which can be expensive to calculate for certain problems.

• Being stochastic, there are no guarantees on the optimality or the quality of the solution.

• If not implemented properly, the GA may not converge towards the optimum solution.

Applications of GAs:

• Optimization 

• Economics

• Neural Networks

• Parallelization 

• Image Processing 

• Vehicle routing problems 

• Scheduling applications 

• Machine Learning 

• Robot Trajectory Generation

• Parametric Design of Aircraft 

• DNA Analysis 

• Multimodal Optimization

• Traveling salesman problem and its applications

BASIC TERMINOLOGY:

• Population − It is a subset of all the possible (encoded) solutions to the given problem. The population for a GA is analogous to the population for human beings except that instead of human beings, we have Candidate Solutions representing human beings.

• Chromosomes − A chromosome is one such solution to the given problem.

• Gene − A gene is one element position of a chromosome.

• Allele − It is the value a gene takes for a particular chromosome.

• Genotype − Genotype is the population in the computation space. The solutions are expressed in the computation space in a way that can be easily interpreted and manipulated by a computing device.

• Phenotype − The population in the actual real world solution space in which solutions are represented in the way they are represented in real world situations is referred to as the phenotype.

• Decoding and Encoding − For simple problems, the phenotype and genotype spaces are the same. However, in most of the cases, the phenotype and genotype spaces are different. Decoding is the process of converting a solution from the genotype space to the phenotype space, whereas encoding is the process of converting from the phenotype space to the genotype space. Decoding should be quick because it is repeated in a GA during the fitness value calculation.

• Fitness Function − A fitness function simply defined is a function which takes the solution as input and produces the suitability of the solution as the output. In some cases, the fitness function and the objective function may be the same, while in others it might be different based on the problem.

• Genetic Operators − These alter the genetic composition of the offspring. These include crossover, mutation, selection, etc.

CROSSOVER:

In crossover, more than one parent is selected and one or more off-springs are produced using the genetic material of the parents. Some of most popularly used crossover operators are discussed below

One Point Crossover

In this one-point crossover, a random crossover point is selected and the tails of its two parents are swapped to get new off-springs.

Multi Point Crossover:

Multi point crossover is a generalization of the one-point crossover wherein alternating segments are swapped to get new off-springs.

Uniform Crossover:

In a uniform crossover, the chromosome is not divided into segments, but rather each gene is treated separately. In this case, we essentially flip a coin for each chromosome to see if it will be included in the offspring. We can also tilt the coin in favour of one parent in order to have more genetic material from that parent in the child.

MUTATION:

A mutation is defined as a small random change in the chromosome that results in a new solution. It is used to keep and add diversity to the genetic population. It has been discovered that mutation is required for the GA to converge, whereas crossover is not.

Some of most popularly used crossover operators are discussed below

Bit Flip Mutation

In this bit flip mutation, we select one or more random bits and flip them. This is used for binary encoded GAs.

Random Resetting

Random Resetting is an extension of the bit flip for the integer representation. In this, a random value from the set of permissible values is assigned to a randomly chosen gene.

Swap Mutation

In swap mutation, we select two positions on the chromosome at random, and interchange the values. This is common in permutation based encodings.

Scramble Mutation

Scramble mutation is also popular with permutation representations. In this, from the entire chromosome, a subset of genes is chosen and their values are scrambled or shuffled randomly.

Inversion Mutation

In inversion mutation, we select a subset of genes like in scramble mutation, but instead of shuffling the subset, we merely invert the entire string in the subset.

SELECTION:

The Survivor Selection Policy determines who will be kicked out and who will be kept in the next generation. It is critical because it ensures that the fitter individuals are not pushed out of the population while also maintaining population diversity.

Elitism is used by some GAs. Simply put, it means that the current fittest member of the population is always passed down to the next generation. As a result, the fittest member of the current population can never be replaced.

The simplest policy is to kick out random members of the population, but this approach frequently leads to convergence issues, so the following strategies are commonly used.

Age based Selection:

Notion of a fitness is not used in Age based selection. It is based on the premise that each individual is allowed in the population for a finite generation during which it is allowed to reproduce, after which it is kicked out regardless of how fit it is.

Fitness Based Selection:

Children tend to replace the least fit individuals in the population in this fitness-based selection. 

TERMINATION CONDITIONS OF GAs:

A Genetic Algorithm's termination condition is critical in determining when a GA run will end. It has been observed that the GA progresses very quickly at first, with better solutions arriving every few iterations, but this tends to saturate in the later stages, where the improvements are very small. Typically, we want a termination condition that ensures our solution is close to optimal at the end of the run.

Termination conditions used for GAs are given below−

• Iterations: When there has been no improvement in the population for X iterations.
• Generations: When we reach an absolute number of generations.
• Fitness Value: When the objective function value has reached a certain pre-defined value.

The termination condition, like other GA parameters, is highly problem specific, and the GA designer should experiment with different options to see what works best for his specific problem.

BASIC STRUCTURE OF GAs:

STALAGMITE FUNCTION:

A stalagmite is a type of rock formation that rises from the cave floor as a result of material deposited on the floor by ceiling drippings.

Generally stalagmite function `f(x,y)` is given as 

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

where,

`f_(1,x) = [sin(5.1pix + 0.5)]^6`

`f_(1,y) = [sin(5.1piy + 0.5)]^6`

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

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

SYNTAX FOR GA IN MATLAB:

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.

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').

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').

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

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

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

x = ga(fun,nvars,A,b,[],[],lb,ub,nonlcon,IntCon) or x = ga(fun,nvars,A,b,[],[],lb,ub,nonlcon,IntCon,options) requires that the variables listed in IntCon take integer values.

x = ga(problem) finds the minimum for problem, a structure described in problem.

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

Syntax for 'options' using 'optimoptions':

options = optimoptions(SolverName) returns a set of default options for the SolverName solver.

options = optimoptions(SolverName,Name,Value) returns options with specified parameters set using one or more name-value pair arguments.

options = optimoptions(oldoptions,Name,Value) returns a copy of oldoptions with the named parameters altered with the specified values.

options = optimoptions(SolverName,oldoptions) returns default options for the SolverName solver, and copies the applicable options in oldoptions to options.

options = optimoptions(prob) returns a set of default options for the prob optimization problem or equation problem.

options = optimoptions(prob,Name,Value) returns options with specified parameters set using one or more name-value pair arguments.

%%%%%%%% STALAGMITE FUNCTION %%%%%%%%
function [fx_y] = stalagmite(input)
    x = input(1);
    y = input(2);
    f1_x = (sin((5.1*pi*x) + 0.5))^6;
    f1_y = (sin((5.1*pi*y) + 0.5))^6;
    f2_x = exp(-4*log(2)*((x - 0.0667)^2)/0.64);
    f2_y = exp(-4*log(2)*((y - 0.0667)^2)/0.64);
    % 'ga' function in MATLAB gives minimum value of the function.
    % In order to obtain maximum value, -ve of the function is given as input
    % -ve sign is used here to find out maxima of function using 'ga' function. 
    fx_y = -(f1_x * f1_y * f2_x * f2_y);
end

% GENERIC ALGORITHM  FOR STALAGMITE FUNCTION
clear all
close all 
clc

% Number of Iterations
limit = 100;

% Defining our search engine
x = linspace(0.0,0.6,250);
y = linspace(0.0,0.6,250);

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%% CASE 1 : STATISTICAL BEHAVIOUR WITH UNBOUNDED CONSTRAINTS %%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tic                 % Sarting of timer for case 1
for i = 1:limit
    % [x, fval] = ga(fun,nvars)
    % x - unconstrained minimum
    % fval - fitness value of function at x
    % nvars is the number of design variables of function stalagmite. Here x,y = 2.
    [inputs, fval(i)] = ga(@stalagmite, 2);
    x_out1(i) = inputs(1);
    y_out1(i) = inputs(2);
end
fval_case1 = -fval;
case1_time = toc    % End of timer for calculating case 1

% 'ga' function in MATLAB gives minimum value of the function.
% In order to obtain maximum value, -ve of the function is given as input
figure(1)
subplot(4,1,1:3)
% Plotting the Stalagmite surface
% Since 'ga' function in MATLAB gives minimum value for input function 'f',
% The function 'f' is given as -ve of that fuction for input to 'ga' function.
% Hence, while plotting function 'f', -ve of that function is plotted to get original function. 
surfc(xx,yy,-f)
% VShading the color in each line segment and face by interpolating the colormap index across the line or face
shading interp
xlabel('x');
ylabel('y');
title('CASE 1: STATISTICAL BEHAVIOUR WITH UNBOUNDED CONSTRAINTS');
hold on
% plotting the location of Maximum value in each iteration 
% Since 'fval' from 'ga' is minimum of the -ve function, -ve of that value is the maximum of +ve function
plot3(x_out1, y_out1, fval_case1, 'marker', 'o', 'markersize', 6, 'markerfacecolor', 'r');
xlabel('x');
ylabel('y');
subplot(4,1,4)
% Plotting Maximum value in each iteration
plot(-fval,'linewidth', 3);
axis([0 100 -0.25 1.25]);
xlabel('No of Iterations');
ylabel('Max F value');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% CASE 2 : STATISTICAL BEHAVIOUR WITH BOUNDED CONSTRAINTS WITH DEFAULT POPULATION SIZE OF 50 %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tic                 % Sarting of timer for case 2
for i = 1:limit
    % [x, fval] = ga(fun,nvars,A,b,Aeq,beq,lb,ub)
    % 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 - unconstrained minimum
    % fval - fitness value of function at x
    % nvars is the number of design variables of function stalagmite. Here x,y = 2.
    [inputs, fval(i)] = ga(@stalagmite,2,[],[],[],[],[0;0],[0.6;0.6]);
    x_out2(i) = inputs(1);
    y_out2(i) = inputs(2);
end
fval_case2 = -fval;
case2_time = toc    % End of timer for calculating case 2

% 'ga' function in MATLAB gives minimum value of the function.
% In order to obtain maximum value, -ve of the function is given as input
figure(2)
subplot(4,1,1:3)
% Plotting the Stalagmite surface
% Since 'ga' function in MATLAB gives minimum value for input function 'f',
% The function 'f' is given as -ve of that fuction for input to 'ga' function.
% Hence, while plotting function 'f', -ve of that function is plotted to get original function. 
surfc(xx,yy,-f)
% VShading the color in each line segment and face by interpolating the colormap index across the line or face
shading interp
xlabel('x');
ylabel('y');
title({['CASE 2: STATISTICAL BEHAVIOUR WITH BOUNDED CONSTRAINTS'];['Population Size = 50']});
hold on
% plotting the location of Maximum value in each iteration 
% Since 'fval' from 'ga' is minimum of the -ve function, -ve of that value is the maximum of +ve function
plot3(x_out2, y_out2, fval_case2, 'marker', 'o', 'markersize', 6, 'markerfacecolor', 'r');
xlabel('x');
ylabel('y');
subplot(4,1,4)
% Plotting Maximum value in each iteration
plot(-fval,'linewidth', 3);
axis([0 100 0.2 1.2]);
xlabel('No of Iterations');
ylabel('Max F value');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%% CASE 3 : STATISTICAL BEHAVIOUR WITH BOUNDED CONSTRAINTS WITH POPULATION SIZE OF 100 %%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
options = optimoptions('ga');
options = optimoptions(options,'populationsize',100);

tic                 % Sarting of timer for case 3
for i = 1:limit
    % [x, fval] = ga(fun,nvars,A,b,Aeq,beq,lb,ub,nonlcon,options)
    % Linear equalities Aeq*x = beq and A*x ? b. 
    % Set A=[] and b=[] if no linear inequalities exist.
    % Set nonlcon = [] if no nonlinear inequalities exist.
    % ga evaluates the matrix product Aeq*x as if x is transposed (Aeq*x').
    % x - unconstrained minimum
    % fval - fitness value of function at x
    % nvars is the number of design variables of function stalagmite. Here x,y = 2.
    [inputs, fval(i)] = ga(@stalagmite,2,[],[],[],[],[0;0],[0.6;0.6],[],options);
    x_out3(i) = inputs(1);
    y_out3(i) = inputs(2);
end
fval_case3 = -fval;
case3_time = toc    % End of timer for calculating case 3

% 'ga' function in MATLAB gives minimum value of the function.
% In order to obtain maximum value, -ve of the function is given as input
figure(3)
subplot(4,1,1:3)
% Plotting the Stalagmite surface
% Since 'ga' function in MATLAB gives minimum value for input function 'f',
% The function 'f' is given as -ve of that fuction for input to 'ga' function.
% Hence, while plotting function 'f', -ve of that function is plotted to get original function. 
surfc(xx,yy,-f)
% VShading the color in each line segment and face by interpolating the colormap index across the line or face
shading interp
xlabel('x');
ylabel('y');
title({['CASE 3: STATISTICAL BEHAVIOUR WITH BOUNDED CONSTRAINTS'];['Population Size = 100']});
hold on
% plotting the location of Maximum value in each iteration 
% Since 'fval' from 'ga' is minimum of the -ve function, -ve of that value is the maximum of +ve function
plot3(x_out3, y_out3, fval_case3, 'marker', 'o', 'markersize', 6, 'markerfacecolor', 'r');
xlabel('x');
ylabel('y');
subplot(4,1,4)
% Plotting Maximum value in each iteration
plot(-fval,'linewidth', 3);
axis([0 100 0.5 1.5]);
xlabel('No of Iterations');
ylabel('Max F value');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%% CASE 4 : STATISTICAL BEHAVIOUR WITH BOUNDED CONSTRAINTS WITH POPULATION SIZE OF 1000 %%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
options = optimoptions('ga');
options = optimoptions(options,'populationsize',1000);

tic                 % Sarting of timer for case 4
for i = 1:limit
    % [x, fval] = ga(fun,nvars,A,b,Aeq,beq,lb,ub,nonlcon,options)
    % Linear equalities Aeq*x = beq and A*x ? b. 
    % Set A=[] and b=[] if no linear inequalities exist.
    % Set nonlcon = [] if no nonlinear inequalities exist.
    % ga evaluates the matrix product Aeq*x as if x is transposed (Aeq*x').
    % x - unconstrained minimum
    % fval - fitness value of function at x
    % nvars is the number of design variables of function stalagmite. Here x,y = 2.
    [inputs, fval(i)] = ga(@stalagmite,2,[],[],[],[],[0;0],[0.6;0.6],[],options);
    x_out4(i) = inputs(1);
    y_out4(i) = inputs(2);
end
fval_case4 = -fval;
case4_time = toc    % End of timer for calculating case 4

% 'ga' function in MATLAB gives minimum value of the function.
% In order to obtain maximum value, -ve of the function is given as input
figure(4)
subplot(4,1,1:3)
% Plotting the Stalagmite surface
% Since 'ga' function in MATLAB gives minimum value for input function 'f',
% The function 'f' is given as -ve of that fuction for input to 'ga' function.
% Hence, while plotting function 'f', -ve of that function is plotted to get original function. 
surfc(xx,yy,-f)
% VShading the color in each line segment and face by interpolating the colormap index across the line or face
shading interp
xlabel('x');
ylabel('y');
title({['CASE 4: STATISTICAL BEHAVIOUR WITH BOUNDED CONSTRAINTS'];['Population Size = 1000']});
hold on
% plotting the location of Maximum value in each iteration 
% Since 'fval' from 'ga' is minimum of the -ve function, -ve of that value is the maximum of +ve function
plot3(x_out4, y_out4, fval_case4, 'marker', 'o', 'markersize', 6, 'markerfacecolor', 'r');
xlabel('x');
ylabel('y');
subplot(4,1,4)
% Plotting Maximum value in each iteration
plot(-fval,'linewidth', 3);
axis([0 100 0.5 1.5]);
xlabel('No of Iterations');
ylabel('Max F value');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% STALAGMITE FUNCTION
function [fx_y] = stalagmite(input)
    x = input(1);
    y = input(2);
    f1_x = (sin((5.1*pi*x) + 0.5))^6;
    f1_y = (sin((5.1*pi*y) + 0.5))^6;
    f2_x = exp(-4*log(2)*((x - 0.0667)^2)/0.64);
    f2_y = exp(-4*log(2)*((y - 0.0667)^2)/0.64);
    % 'ga' function in MATLAB gives minimum value of the function.
    % In order to obtain maximum value, -ve of the function is given as input
    % -ve sign is used here to find out maxima of function using 'ga' function. 
    fx_y = -(f1_x * f1_y * f2_x * f2_y);
end