To maximize the stalagmite function using genetic algorithm in MATLAB

Aim:
To optimise the stalagmite function and obtain the global maxima for the given function

Theory:
Genetic algorithm is a sort of heuristic or a probabilistic search inspired by Charles Darwin’s theory of natural selection. The process of natural selection starts with the selection of fittest individuals from a population. 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. Using this principle, the genetic algorithm function selects the inputs randomly from the initial workspace and determines the local maximas. The process repeats itself again and again and at the end, a global-maxima is achieved.

Let us consider that there is a function that is unsolvable or too difficult to solve using analytical methods and thus its optimisation becomes very hard and time consuming. In such cases, genetic algorithm proves to be very helpful. It generates an initial population of random solutions. These solutions are then evaluated based or their fitness to the fitness function, or the function under consideration. The selected population then breeds within itself to produce a fitter population, that is, better solutions to the function. Thus, each subsequent set of solution would be better than the previous one and the process is iterated till the fittest population, in other words, the optimum value to the function is obtained.

The population (solutions) are evaluated based on the environment, to give better off-breeds (local minimas or maximas) which eventually provide the desired results. Thus, genetic algorithm is considered to be inspired for the theory of natural selection, or, the survival of the fittest.

Syntax of ga function in MATLAB
[inputs, f] = ga(fun, nvars)
fun = function to be optimised
nvars = number of associated variables/ inputs

[inputs, f] = ga(fun, nvars, A, b, Aeq, beq, lb, ub, nonIcons, IntCon, options)
Inputs is a matrix having the values of inputs to the function corresponding to the optimum value and f is the optimum value of the function
A, b = linear inequality constraints. These are the matrix of coefficients of the inputs to the functions, satisfying certain linear inequalities.
Aeq, beq = linear equality constraints. These are the matrix of coefficients of the inputs to the functions, satisfying certain linear equalities.
lb, ub = It refers to lower and upper bounds respectively. These are used to set limits to the values of input variables. For example, in the program given below, 0<=x<=0;6 and 0<=y<=0;6. Thus we provide these to ga function as [0;0], [0.6;0.6]
nonlcon: Use for nonlinear constraints
IntCon: to get solution in integer.
options: It replaces the default optimization parameters. It is created using optimoptions command.

Options = optimoptions(‘function’, ‘constraint’, ‘limit’)

Note: If any of the above values need not be provided to the ga function or is not applicable, it is replaced by [ ] in the function syntax.

Program

Initially, the inputs are assigned the values and a 3D array is generated from these values using meshgrid command. The xx and yy matrices so obtained are stored in column 1 and 2 respectively of the array, input_vector and is used as input in the stalagmite function.

clear all
close all
clc

%defining the inputs
x = linspace(0,0.6,150);
y = linspace(0,0.6,150);

% creating 3D array
[xx yy] = meshgrid(x,y);

%looping 

for i= 1:length(x)
   for j = 1:length(y)
   input_vector(1) = xx(i,j);
   input_vector(2) = yy(i,j);
   f(i,j) = stalagmite(input_vector);
   end
end

%plotting
surfc(x,y,-f)
xlabel('x value')
ylabel('y value')
shading interp

% study 1 unbounded inputs study
tic
num_cases = 100
for i = 1:num_cases
    [inputs, f_opt(i)] = ga(@stalagmite,2)
     x_opt(i) = inputs(1);
     y_opt(i) = inputs(2);
end
study_time1 = toc %to calculate compilation time

%plotting
figure(2)
subplot(2,1,1)
hold on
surfc(x,y,-f);
shading interp
plot3(x_opt,y_opt,-f_opt,'marker','o','MarkerSize',5,'MarkerFaceColor','r')
title('unbounded inputs')
subplot(2,1,2)
plot(-f_opt)
xlabel('iterations')
ylabel('function maximum')


%study 2 setting upper and lower bounds to the inputs
tic
for i = 1:num_cases
    [inputs, f_opt(i)] = ga(@stalagmite,2,[],[],[],[],[0;0],[0.6;0.6]);
     x_opt(i) = inputs(1);
     y_opt(i) = inputs(2);
end
study_time2 = toc

%plotting
figure(3)
subplot(2,1,1)
hold on
surfc(x,y,-f);
shading interp
plot3(x_opt,y_opt,-f_opt,'marker','o','MarkerSize',5,'MarkerFaceColor','r')
title('bounded inputs')
subplot(2,1,2)
plot(-f_opt) %plot optimum value obtained in each iteration
xlabel('iterations')
ylabel('function maximum')

% study3 defining the population size
options = optimoptions('ga');
options = optimoptions(options,'PopulationSize',200);
tic
for i = 1:num_cases
    [inputs, f_opt(i)] = ga(@stalagmite,2,[],[],[],[],[0;0],[0.6;0.6],[],[],options);
     x_opt(i) = inputs(1);
     y_opt(i) = inputs(2);
end
study_time3 = toc

%plotting
figure(4)
subplot(2,1,1)
hold on
surfc(x,y,-f)
shading interp
plot3(x_opt,y_opt,-f_opt,'marker','o','MarkerSize',5,'MarkerFaceColor','r')
title('bounded inputs with mod. population size')
subplot(2,1,2)
plot(-f_opt)
xlabel('iterations')
ylabel('function maximum')

% Getting the optimum value
f_max = -min(f_opt);
fprintf('maximum value of the function is')
disp(f_max)

In the above program, plot3 command is used to plot 3 variables, x_opt, y_opt and -f_opt in the same plot to determine the optimum value of the function (local as well as global) and corresponding values of the inputs x and y.

Stalagmite function

function [ fxy ] = stalagmite(input_vector)
x = input_vector(1);
y = input_vector(2);
fx = (sin(5.1 * pi * x + 0.5))^6;
fy = (sin(5.1 * pi * y + 0.5))^6;
gx = exp((-4) *log(2) * (((x - 0.0667)^2)/0.64));
gy = exp((-4) *log(2) * (((y - 0.0667)^2)/0.64));
fxy= -(fx * gx * fy * gy); 
end

Generic algorithm, ga function in MATLAB minimizes the function. But here we are required to maximise the function and get the global maxima. In order to do so, we use the ga function to minimize the negative of the function, which is same as maximizing the function. Hence, negative sign (-) is used in the coding of stalagmite function (-fx * fy * gx * gy).

For the same reason, negative sign is also used while plotting the results in the main program.

Output

Study 1

In this study, the ga function was executed without setting the lower and upper bounds to the inputs due to which, we get a scattered plot with no particular value of a global maxima. This means that the global maxima varied for each iteration.

Study 2

In this study, the lower and upper bounds were set to the input variables. This constrained the scattered plot obtained in the study 1 to a concise area, i.e. the range of the global maxima was more precise than that obtained in the previous study.

Study 3

In this study, the default optimization parameters were altered to get better results. The population size was increased using optimoptions command, which by-default, is set to 50 by MATLAB. As the population size increases, the possible combination for the input values also increases which helps in getting an accurate, or close to accurate value of the maxima of the function.

However, there is a problem associated with this approach, that is, it increases the execution time quite significantly.

Results

According to study 3, we get an accurate value of global maxima for the stalagmite function. The value of global maxima is equal to 1.

Errors

Error 1

optimoptions command requires a function as its input. This function is to be provided inside ''. This error was generated because while writing the optimoptions command, the ga fuunction was passed as it is, instead of 'ga'.

error 2

The above error was generated due to syntax error while passing the constraint to the optimoptions.

Conclusion

In order to maximize a function, we just need to minimize the negative of the function. Genetic algorithm is a kind of heuristic search for the optimum value; therefore, it is possible that it may lead to results that are good enough for the function but are not accurate. Also, the results depend upon the size of initial population. The precision of the result improves with greater population size or a greater number of iterations.

As the population size increases, number of combinations to get a fitter solution also increases which increases the accuracy of the results. However, it also increases the time required to get the results. Therefore, the population size should be as low as possible.