Stalagmite Function Optimization using Genetic Algorithm

• Selection rules select the individuals, called parents, that contribute to the population at the next generation. For Eg: parent1 = [1010011000], parent2 = [1001001010] are two points/solutions which are good fits to the objective fucntion and thus are selected for the next generation.

• Crossover rules combine two parents to form children for the next generation.For Eg: parent1 = [1010011000], parent2 = [1001001010], child = [1000011010]

• Mutation rules apply random changes to individual parents to form children. For Eg: parent = [1010011000], child = [0101010001]

In each generation, the fitness of every individual candidate solution in the population is evaluated. Fitness usually is decided on the basis of the value of the objective function for a candidate solution.

The more fit solutions are stochastically (randomly) selected from the current populations and each solution property is modified (recombined/cross-over or mutated) to form a new generation. The new generation of candidate solutions is then used in the next iteration of GA.

The termination point of GA occurs at two instances, either when the maximum number of generations have been obtained or when satisfactory levels of fitness have been reached.

METHODOLOGY:

Ist Iteration:

Selection:

IInd Iteration (generation formation by either mutation or recombination)

Termination/Convergence:

Syntax for GA in MATLAB:

x = ga(fun,nvars,A,b) finds the design input variables x corresponding to a local minimum fun, subject to the linear inequalities A*x ≤ b. ga evaluates the matrix product A*x as if x is transposed (A*x').

where x = `[[x(1)],[x(2)],[x(3)],[x(4)],[.],[.]]` where x(1),x(2),x(3) etc. are the design variables for the objective function. (take A = [] and b = [] if there are no linear inequality constraints)

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 Aeq=[] and beq=[] if no linear equality constraints 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 for the range lb ≤ x ≤ ub. (Set lb=[] and ub=[] if no bounds exist.)

x = ga(fun,nvars,A,b,Aeq,beq,lb,ub,nonlcon) subjects the minimization to the constraints defined in nonlcon (non-linear constraints). 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 nonlcon=[] if no nonlinear constraints exist.)

x = ga(fun,nvars,A,b,Aeq,beq,lb,ub,nonlcon,options) minimizes with the default optimization parameters replaced by values in options i.e. the minimization is done with the modified paramters from options.

x = ga(fun,nvars,A,b,Aeq,beq,lb,ub,nonlcon,Intcon,options) requires that the variables listed in IntCon take integer values. (Set IntCon = [] if there exist no interger constraints)

OPTIMIZING THE STALAGMITE FUNCTION IN ORDER TO FIND GLOBAL MAXIMA:

The stalagmite function is given by:

Since, we are interested in the global maxima of the function but the GA algorithm in MATLAB only calculates the local/global minima, we have to modify the function. In order to get the global maxima, our stalagmite function will now look like:

`f(x,y) = -underset(1,x)(f) . underset(2,x)(f) . underset(1,y)(f) . underset(2,y)(f)`

The 3D plot for the modified stalagmite function looks like:

The code for the stalagmite function is given by:

function [func] = stalagmite(variables)

variables = [variables(1) variables(2)];
term1 = (sin(5.1*pi*variables(1) + 0.5))^6;
term2 = (sin(5.1*pi*variables(2) + 0.5))^6;
term3 = exp(-4*log(2)*((variables(1)-0.0667)^2)/0.64);
term4 = exp(-4*log(2)*((variables(2)-0.0667)^2)/0.64);
func = -(term1*term2*term3*term4);%minus sign added in order to find out the global maxima of the staligmite function with the help of the geetic algo which calculates the minima only
end

Now, since the function has been inverted, the minimum values found by the GA will lie on the peaks, and the modulus of those minimum values will give the maximas.

We will be writing the code in 3 different ways:

1). The GA opertated in an unbound manner

2). The design variables are bounded

3). The design variables are bounded and the population size for GA has been increased.

We will first define a search space for the GA to operate on. Note that this search space is not automatically detected by the GA. This search space also acts as the design variable space for the objective fucntion i.e. the stalagmite function.

%sample space for the GA (but not recognised by it automatically)
%defining 2 row vectors
x = linspace(0,0.6,150);
y = linspace(0,0.6,150);

%creating 2-D arrays
[xx yy] = meshgrid(x,y);

%giving different values to the 2 input variables in stalagmite function
for i=1:length(xx)
    for j=1:length(yy)
        input_var(1) = xx(i,j);
        input_var(2) = yy(i,j);
        f(i,j) = stalagmite(input_var);
    end
end

Now, we will be applying the GA on this function for a fixed number of times. Our aim is to ensure that the GA gives the same optimal soution every time it is applied.

num_cases = 50; %the number of times GA will be run (number of iterations)
%inbuilt MATLAB genetic algo syntax for finding minima of a function
for i=1:num_cases
    %objective function is the stalagmite function
    %number of design variables for the objective function is 2
    %A = [], b = [] as there are no constraints on the objective function
    %with inequalities
    %Aeq = [], beq = [] as there are no constraints on the objective function
    %with inequalities
    %no lower bound for both the input variables (xx(i,j) and yy(i,j)) lb = []
    %no upper bound for both the input variables (xx(i,j) and yy(i,j)) ub ] = []
    %nonlcon = [] for no nonlinear constraints
    %Intcon = [] for no interger constraints
    [opt_inputs,f_val(i)] = ga(@stalagmite,2);%gives the optimum values of inputs for 
     which function is min and min value of function
     opt_f(i) = -f_val(i); %maximum value of the stalagmite function
     %optimal design variables for which the stalagmite function minimises
     xxopt(i) = opt_inputs(1);
     yyopt(i) = opt_inputs(2);
end

In the above code, we can see that for each iteration in the for loop, the GA is giving an optimal maxima (opt_f 1x50 matrix) and 2 optimal design variables (xxopt and yyopt, both 1x50 matrices).

If we plot the 3D surface plot for stalagmite function and the optimal values obtained by the GA function along with the change in optimal value with changing iterations, we will get the following plots:

As seen above, since the GA is unbounded, it gives optimal solutions which vary with each iteration and are way off the downward peaks.

The code for plotting is given below:

figure(2)
%plotting to visualise the optimal value obtained from GA with respect to
%the 3D surface plot of the stalagmite function
subplot(2,1,1)
hold on
surfc(xx,yy,f)
shading interp
plot3(xxopt,yyopt,f_val,'-o','markers',5,'markerface','r')
xlabel('xx input variable')
ylabel('yy input variable')
title('Unbounded Inputs')

%plotting to understand how the value of optimal solution for the
%stalagmite function changes with each use of ga algorithm (iteration)
subplot(2,1,2)
plot(opt_f)
xlabel('Iterations')
ylabel('Maximum Value of Function')
      

2). If we bind the design variables with a lower bound of 0 and an upper bound of 0.6, we will be able to increase the accuracy of GA.

for i=1:num_cases
    %objective function is the stalagmite function
    %number of design variables for the objective function is 2
    %A = [], b = [] as there are no constraints on the objective function
    %with inequalities
    %Aeq = [], beq = [] as there are no constraints on the objective function
    %with inequalities
    %lower bound for both the input variables (xx(i,j) and yy(i,j)) is 0
    %upper bound for both the input variables (xx(i,j) and yy(i,j)) is 0.6
    %nonlcon = [] for no nonlinear constraints
    %Intcon = [] for no interger constraints
    [opt_inputs,f_val(i)] = ga(@stalagmite,2,[],[],[],[],[0;0],[0.6;0.6]);%gives the optimum values of inputs for which function is min and min value of function
     opt_f(i) = -f_val(i); %maximum value of the stalagmite function
     %optimal design variables for which the stalagmite function minimises
     xxopt(i) = opt_inputs(1);
     yyopt(i) = opt_inputs(2);
end

In the above code, the lower and upper bounds for the design variables has been set as 0 and 0.6 respectively. If we plot the optimal outputs and the 3D plots now, we get much accurate plots:

As seen in the above plot, even though the accuracy has increased, the GA stills gives values of local maxima for certain iterations instead of global maxima.

3). In order to further increase the accuracy, we will increase the population size for the GA function using the gaoptimset function. Though this function is not as accurate as its later versions (optionoptim) but still we can achieve a global maxima using gaoptimset.

num_cases = 50; %the number of times GA will be run (number of iterations)
%inbuilt MATLAB genetic algo syntax for finding minima of a function
%increasing the population size of the search space for the ga algorithm
options = gaoptimset('PopulationSize',500)%population size should be less than the maximum size of the objective function (150x150 in this case)
for i=1:num_cases
    %objective function is the stalagmite function
    %number of design variables for the objective function is 2
    %A = [], b = [] as there are no constraints on the objective function
    %with inequalities
    %Aeq = [], beq = [] as there are no constraints on the objective function
    %with inequalities
    %lower bound for both the input variables (xx(i,j) and yy(i,j)) is 0
    %upper bound for both the input variables (xx(i,j) and yy(i,j)) is 0.6
    %nonlcon = [] for no nonlinear constraints
    %Intcon = [] for no interger constraints
    [opt_inputs,f_val(i)] = ga(@stalagmite,2,[],[],[],[],[0;0],[0.6;0.6],[],[],options);%gives the optimum values of inputs for which function is min and min value of function
     opt_f(i) = -f_val(i); %maximum value of the stalagmite function
     %optimal design variables for which the stalagmite function minimises
     xxopt(i) = opt_inputs(1);
     yyopt(i) = opt_inputs(2);
end

NOTE: The population size for the GA, at all times should be less than the maximum size of range of the objective function, 150x150 in this case.

The plots obtained are:

If we look closely in the above plots, the maximum optimal value of the stalagmite function remains 1 (global maxima) for all its iterations.

The optimal outputs achieved are given below:

Thus now, the GA is working accurately so as to find the global maxima of the stalagmite function.

The full code with the appropriate comments is given below:

clear all
close all
clc

%sample space for the GA (but not recognised by it automatically)
%defining 2 row vectors
x = linspace(0,0.6,150);
y = linspace(0,0.6,150);

%creating 2-D arrays
[xx yy] = meshgrid(x,y);

%giving different values to the 2 input variables in stalagmite function
for i=1:length(xx)
    for j=1:length(yy)
        input_var(1) = xx(i,j);
        input_var(2) = yy(i,j);
        f(i,j) = stalagmite(input_var);
    end
end

%plotting 3-D plots of the function
surf(xx,yy,f)
shading interp %added to shade the contours
xlabel('X-AXIS')
ylabel('Y-AXIS')
colorbar

num_cases = 50; %the number of times GA will be run (number of iterations)
%inbuilt MATLAB genetic algo syntax for finding minima of a function
%increasing the population size of the search space for the ga algorithm
options = gaoptimset('PopulationSize',500)%population size should be less than the maximum size of the objective function (150x150 in this case)
for i=1:num_cases
    %objective function is the stalagmite function
    %number of design variables for the objective function is 2
    %A = [], b = [] as there are no constraints on the objective function
    %with inequalities
    %Aeq = [], beq = [] as there are no constraints on the objective function
    %with inequalities
    %lower bound for both the input variables (xx(i,j) and yy(i,j)) is 0
    %upper bound for both the input variables (xx(i,j) and yy(i,j)) is 0.6
    %nonlcon = [] for no nonlinear constraints
    %Intcon = [] for no interger constraints
    [opt_inputs,f_val(i)] = ga(@stalagmite,2,[],[],[],[],[0;0],[0.6;0.6],[],[],options);%gives the optimum values of inputs for which function is min and min value of function
     opt_f(i) = -f_val(i) %maximum value of the stalagmite function
     %optimal design variables for which the stalagmite function minimises
     xxopt(i) = opt_inputs(1);
     yyopt(i) = opt_inputs(2);
end


figure(2)
%plotting to visualise the optimal value obtained from GA with respect to
%the 3D surface plot of the stalagmite function
subplot(2,1,1)
hold on
surfc(xx,yy,f)
shading interp
plot3(xxopt,yyopt,f_val,'-o','markers',5,'markerface','r')
xlabel('xx input variable')
ylabel('yy input variable')
title('Bounded Inputs')

%plotting to understand how the value of optimal solution for the
%stalagmite function changes with each use of ga algorithm (iteration)
subplot(2,1,2)
plot(opt_f)
xlabel('Iterations')
ylabel('Maximum Value of Function')
      
    

      

REFERENCES:

https://in.mathworks.com/help/gads/what-is-the-genetic-algorithm.html#targetText=The%20genetic%20algorithm%20is%20a,a%20population%20of%20individual%20solutions.

https://en.wikipedia.org/wiki/Genetic_algorithm

https://in.mathworks.com/help/gads/ga.html#mw_5430547e-4111-4753-a80d-727f85b8f63e

https://in.mathworks.com/help/optim/ug/optim.problemdef.optimizationproblem.optimoptions.html