Global optimization of Stalagmite function using Genetic Algorithm

Objective:

To write a code for optimization of the stalagmite function and finding the global maxima value using Genetic Algrithm.

Description:

Stalactite and Stalagmite, elongated forms of various minerals deposited from solution by slowly dripping water. A stalactite hangs like an icicle from the ceiling or sides of a cavern. A stalagmite appears like an inverted stalactite, rising from the floor of a cavern.

The Stalagmite function is written as,

For a Stalagmite function to run from a given set of inputs, we need to define the input parameters which are 2 dimensional arrays. Here, we have defined x and y which are sets of 1 dimensional array and having values between 0 to 0.6, each having 150 variables which are equally spaced. This means x and y are row vector of the form 1X150.

In order to form a stalagmite formation we require 2 dimensional arrays, and they are formed from the given input using the meshgrid command, and in the process we obtain xx and yy which are the matrices of the form 150x150 (2 dimensional array) from the given input x and y. Implementing the value of x and y in the above equations, we can obtain the function value for a given x and y coordinate.

Now, we can plot the surface for the stalagmite function with the values of xx, yy and Function value (F). From the plot, we can observe the surface in the form of stalagmite formed which consist of various local maxima and a global maxima. In order to plot the location and calculating the coordinate of the function value at the global maxima we need to incorporate the study of 'Genetic Algorithm'. But the main property of the genetic algorithm is to minimize a particular function and the value obtained will be the value of the minimum point on the surface. So, for calculating value of function at global maxima we need to make some changes in the function value of the Stalagmite function. If we make the main function value negative then we can obtain the value of the global maxima and its coordinate.

What is Genetic Algorithm?

Genetic algorithm is a search heuristic which is inspired by the process of natural selection based on the Darwin's Theory of Evolution. The algorithm is used to generate high quality solutions or population to optimization or search problem. The algorithm uses biologically inspired operators like crossover, mutation and selection for obtaining the quality solutions. This operators apply random changes to the solution available to the function in order to increase the fitness value of the solution.

In genetic algorithm, a population is referred to a group of solutions and each individual solution has a set of properties which combined together are known as chromosome. Each property is known as gene, this set of genes are contained in the chromosome. This gene are represented as bianry values of '0' and '1' which are joined together in a form of binary string. Other type of representation system are also possible and they vary according to the problem to be optimized.

In genetic algorithm, the evolution starts from a population of randomly generated individual and it's an iterative process, where population in each iteration is kown as generation. The population size depends upon the nature of problem, but usually contains hundreds and thousands of solutions. The fitness of each of the individual in the population is evaluated with help of the fitness function associated with the optimization problem. Because of the fitness function, solutions with high fitness value are selected for mating to produce an offspring. The fitness of offspring generated by mating of two high quality parent is also high since the offspring will carry the properties of both the parents. This transfer of properties into the offspring occurs due to crossover and this crossover of properties is dependent on the cossover point. The crossover point vary from individual to individual and in addition to crossover, mutation also takes place. Mutation is necessary to avoid occurance of offspring with properties same as that of their parent which ensures that the flaws of the parent will not be continued in the next generation.

The solutions obtained from the first generation are forwarded to next generation so that high quality individual are available for further mating and fitness of the population is improved from one generation to another. The algorithm terminates the process either when maximum number of generations are produced or when the fitness level of the population has reched to a desired point.

We have various syntaxes for GA in matlab, those are as follows,

• 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.) Create options using optimoptions.
• 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(fun,nvars,A,b,[],[],lb,ub,nonlcon,IntCon,options) minimizes with the default optimization parameters replaced by values in options. (Set nonlcon=[] if no nonlinear constraints exist.) Create options using optimoptions.
• x = ga(problem) finds the minimum for problem, where problem is a structure.
• [x,fval] = ga(___), for any previous input arguments, also returns fval, the value of the fitness function at x.
• [x,fval,exitflag,output] = ga(___) = 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(___) = ga(___) also returns a matrix population, whose rows are the final population, and a vector scores, the scores of the final population.

The program for optimizing the function and calculating the global maxima an be written as,

% Program for Global optimization problem using Stalagmite function and genetic algorithm

clear all
close all
clc

% Preparing input database

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


% Creating 2 Dimensional array
[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

surfc (xx,yy,F)

shading interp

[result,Fval] = ga(@stalagmite,2);

hold on 

plot3(result(1),result(2),Fval,'marker','o','markersize',3,'markerfacecolor','r')

The surface plot for the Stalagmite function can be observed as,

The value of function after running the Genetic algorithm and its respective coordinate are,

From above images, we can establish that the value of the Function obtained is not that of the global maxima, and every time we run the program we obtain different value for the same function.

So, if we run the program for n, number of times we will obtain n different values of the function ,one of which might be or might not be the value of the global maxima. The value obtained are not linear and scattered over the space. To ensure that we get the same value of the global maxima for n number of times the program is run, then we need to make some addition to the original program. Let us study in three steps in which we can overcome the above mentioned issue.

The program for obtaining the same value of function for the global maxima consistently can be written as,

% Program for Global optimization of Stalagmite function using genetic algorithm

clear all
close all
clc

% Preparing input database

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


% Creating 2 Dimensional array
[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

num_cases = 100;

% Case 1 - Stalagmite function without any bounds

for m = 1:num_cases
    
    [results, Fval(m)] = ga(@stalagmite,2);
    x_val(m) = results(1);
    y_val(m) = results(2);
    
end

% Plotting of function value without any bounds

figure(1)
subplot (2,1,1)
hold on
surfc (xx,yy,F)
shading interp
plot3(x_val,y_val,Fval,'marker','o','markersize',3,'markerfacecolor','r')
title('Stalagmite Function without bounds')
subplot (2,1,2)
plot(Fval,'color','b')
axis([-0.1 100 -1.2 0])
xlabel ('Iterations')
ylabel ('Function value')

% Case 2 - Stalagmite function with lower and upper bounds

tic

for h = 1:num_cases
    [results1,Fopt(h)] = ga(@stalagmite,2,[],[],[],[],[0;0],[1;1]);
    xopt(h) = results1(1);
    yopt(h) = results1(2);
end
 
study_time = toc
 
figure(2)
subplot (2,1,1)
hold on
surfc (xx,yy,F)
shading interp
plot3(xopt,yopt,Fopt,'marker','o','markersize',3,'markerfacecolor','r')
title ('Stalagmite Function with bounds')
subplot (2,1,2)
plot(Fopt)
axis([-0.1 50 -1.2 0])
xlabel ('Iterations')
ylabel ('Function value')

% Case 3 - Stalagmite function with lower and upper bounds and defined population size
 
tic
 
options = optimoptions('ga','PopulationSize',300)
 
for k = 1:num_cases
    [results2, Fopt1(k)] = ga(@stalagmite,2,[],[],[],[],[0;0],[1;1],[],[],options);
    xopt1(k) = results2(1);
    yopt1(k) = results2(2);
end
 
study_time1 = toc
 
figure(3)
subplot (2,1,1)
hold on
surfc (xx,yy,F)
shading interp
plot3(xopt1,yopt1,Fopt1,'marker','o','markersize',3,'color','r')
title ('Function value Plot with Initial Population of 300')
subplot (2,1,2)
plot(Fopt1)
axis([-0.1 100 -1.2 0])
xlabel ('Iterations')
ylabel ('Function value')

We have written the program in three steps so we can observe the improvement in the functioning of the program and the consistency in the outputs obtained for the same.

Case 1:

In the first case, we have ran the program for 100 times and obtained the results without applying any modification in the existing program and plotted the results. First plot depicts the stalagmite function surface plot within which we have then plotted the function value for each case. The subplot in the same figure illustrates the plot of function for the number of iterations. The output for Case 1 can be observed as follows,

We can see that the function value is scattered across the space for each of the number of times the program is run. This is not the result that anyone desires. Our main aim is to calculate the value of function at the global maxima. To achieve this, we need to make some changes to the existing program. Now, we move on to Case 2.

Case 2:

In the second case, we have introduced upper and lower bound in the genetic algorithm. The lower bound has been defined as [0;0] and upper bound has been defined as [1;1], since the variation in the output varies in the range 0 to 1. Along with upper and lower bound other variables such as A, b, Aeqand beq are also the input whose default value has been applied since the solutions are not linear. The main effect of applying the bounds to the algorithm is that, the algorithm will now search the populations containing the values between the bounds. Because of the bounds applied, the solutions or population of each generation are also of good quality or fitness hence improving the overall performance of the algorithm. The results obtained by applying the boundary condition are contained and not scattered as in case 1. The plot for Case 2 are illustrated as follows,

From above results we can establish that, after implementation of the boundary condition the results now are refined than that of case 1 but they are still not that consistent as desired. So, to achieve consistent result we need to make another addition to the program in the form of population size. Lets move to third case.

Case 3:

In this case, we have made one more addition to the existing program i.e using the options command which acts as an additional input to the problem and refines the solution or output. Options utilizes optimoptions as the command. Options and optimoptions are explained as below,

Options are a way of combining a set of name-value pairs. They are useful because they allow you to: