Matlab program to optimize the stalagmite function using the genetic algorithm, GA tool in matlab.
- Aim : To optimize the stalagmite function, finding the function maximum value using the genetic algorithm optimisation tool in matlab.
- Equations & Expressions : The following equations are used to express the Stalagmite function:
`f_(x1) = [sin(5.1pix+0.5)]^6`
`f_(y1) = [sin(5.1piy+0.5)]^6`
`f_(x2) = exp[-4ln(2)((x-0.0667)^2/0.64)]`
`f_(y2) = exp[-4ln(2)((y-0.0667)^2)/0.64]`
The stalagmite function is expressed as `f(x,y) = f_(x1)*f_(x2)*f_(y1)*f_(y2)`
The above equations provides us the single output(z-axis) for providing two inputs(x-axis and y-axis). The function we are optimising is the stalagmite function in which we will find the function's maximum value and which inputs gives us the maximum output.
Stalagmite Function : It is a mathematical model describing the shape of the stalactite(long thin mountains like structures hanging from ceiling of caves) developed by a group consisting professor, doctorate candidate, scientists from University of Arizona. They developed a general equation of motion for the growth of stalactites, they wanted to describe the universal shape and growth, characterstics of these stalactites found in nature, through a system(mathematical equation).
- Objectives & Workflow : In this project we will optimise the stalagmite function to find which inputs will give us the function's maximum value and what is the maximum value. The workflow of the project is as follow:
- First we describe some inputs for the function to provide the output then we calculate the output, after plotting the function we have a reference for our results from the optimisation.
- We use the GA function to do the optimisation of the stalagmite function, plotting the results from it with above plot in the same figure window. We do multiple iterations for getting verified results.
- Now we provide the search constraints to GA in lower and upper bounds to get better results and plot the results with function output plot. We plot the function's ouput with each iteration.
- Further on we provide a certain amount of starting population for the optimisation process giving the GA a better hand to work on and plot the results with same as above to get the results verified.
- Matlab Program
The programming language used is Matlab scripting language.
% Program for the optimisation of the stalagmite function
clc;
close all;
clear global;
%inputs for the stalgamite function
x= linspace(0,0.8,150);
y= linspace(0,0.8,150);
%creating the grid from the inputs for the surf function
[X, Y]=meshgrid(x,y);
%calculating the output of the stalagmite function
for i=1:length(X)
for j=1:length(Y) %using nested for loop to feed each respective x and y value
input(1)=X(i,j);
input(2)=Y(i,j);
f(i,j)=stalagmite(input); % calling the stalagmite function
end
end
%optismising using ga
%Case 1st using simple solver with no constraints or options etc.
for i=1:29 %using to loop to run ga multiple times
fit_fun=@stalagmite; %creating a function handle for the function
[result, fun_value(i)]=ga(fit_fun,2); %calling the ga function
x_result(i)=result(1); %storing the x value separately
y_result(i)=result(2); %storing the y value separately
end
%plotting the stalagmite function and x-y value and function value from ga
figure(1);
subplot(2,1,1);
surf(X,Y,f); %plotting the stalagmite function
shading interp;
title('progress 1 stalagmite optimisation');
hold on;
%plotting the x-y value and function value from ga solver using plot3
plot3(x_result,y_result,fun_value,'marker','o','markersize',5,'markerfacecolor','r');
xlabel('x-axis');
ylabel('y-axis'); %labelling the respective axis
zlabel('z-axis');
subplot(2,1,2);
plot(fun_value); %plotting the function value from each iteration
xlabel('iteration'); %labelling the respective axis
ylabel('optimised value')
%Case 2nd now giving upper and lower bounds for better search space
for i=1:29
fit_fun=@stalagmite; %creating a function handle for the function
%giving the ga solver lower and upper bounds for search restrictions
[result, fun_value(i)]=ga(fit_fun,2,[],[],[],[],[0 0],[0.8 0.8]);
x_result(i)=result(1); %storing the x value separately
y_result(i)=result(2); %storing the y value separately
end
%plotting the stalagmite function and x-y value and function value from ga
figure(2);
subplot(2,1,1);
surf(X,Y,f); %plotting the stalagmite function
shading interp;
title('progress 2 stalagmite optimisation');
hold on;
%plotting the x-y value and function value from ga solver
plot3(x_result,y_result,fun_value,'marker','o','markersize',5,'markerfacecolor','r');
xlabel('x-axis');
ylabel('y-axis'); %labelling the respective axis
zlabel('z-axis');
subplot(2,1,2);
plot(fun_value); %plotting the function value from each iteration
xlabel('iteration'); %labelling the respective axis
ylabel('optimised value')
%Case 3rd giving population size with more strict lower and upper bounds
for i=1:53
fit_fun=@stalagmite;
%giving an intial certain popluation size with lower and upper bounds
options = optimoptions('ga','populationsize',190);
[result, fun_value(i)]=ga(fit_fun,2,[],[],[],[],[0 0],[0.4 0.4],[],[],options);
x_result(i)=result(1);
y_result(i)=result(2);
end
%plotting the stalagmite function and x-y value and function value from ga
figure(3);
subplot(2,1,1);
surf(X,Y,f); %plotting the stalagmite function
shading interp;
title('stalagmite optimisation');
hold on;
%plotting the x-y value and function value from ga solver
plot3(x_result,y_result,fun_value,'marker','o','markersize',5,'markerfacecolor','r');
xlabel('x-axis');
ylabel('y-axis'); %labelling the respective axis
zlabel('z-axis');
subplot(2,1,2);
plot(fun_value); %plotting the function value from each iteration
xlabel('iteration'); %labelling the respective axis
ylabel('optimised value')
maxima=-fun_value(1)
- Explanation/Description :
- We use the linspace command to describe the inputs for the stalagmite function, creating 150 values for the x-axis and y-axis. Then we use the meshgrid command to create two matrix of the inputs so that we make a grid from that, helping us to plot a 3-d surface plot.
- Now calculating the output of the stalagmite function but to feed the inputs one by one we use the nested loop, providing us a way to feed each coordinate of the x-y plane to the function. We describe the x and y value to input and then calling the stalagmite function giving us the output and storing it in the variable f. The way we store the output in f is such that it is also a matrix like X and Y, we are storing the output in the matrix with respect to the input we took from the grid(x-y plane) for example for the first x and y value we are storing the output in first place in the matrix.
- We start the optimisation using the GA function, we create a function handle(fun_fit) for the evaluating function(stalagmite) according to the syntax of GA. Then we call the ga function giving fit_fun and number of variables, here the fit_fun is the fitness function according to which the GA function is going to assess the results. We are using a loop to run the GA mutiple times to see how every time random population gives us results. The function's maximum value is stored in the variable fun_value and the x and y values are in result.
- The x and y value in the result variable need to be put in different variables(x_result and y_result) as we need to plot the fun_value using the plot3 command.
- Now we start the plotting of the results taken till now, first we take a specify a figure window, then we create two subplots in the same figure window using the subplot command. In the first subplot we plot the output of the stalagmite function using the surf command creating 3-d surface. Then we use the shading interp command to remove the grid lines from the surface plot, also using the title command to give a title to window or plot.
- We use the hold on command so that all plots we wanna plot retains rather than the last one only. We plot the results and function value from GA operation in the same subplot as of stalagmite function so that we can easily verify our results. We use the plot3 command to plot the x-y and function value, giving marker size and color to analyse the results. The labelling of the axises been done.
- The function value we get from each execution of GA is plotted in the second subplot and the axises are labelled too.
- The results we get from the ga are not optimised so to improve the performance we need to provide the search constraints using the lower and upper bounds limits. The syntax of ga this time is a bit different, GA takes a lot of inputs to work at its best which has a defined syntax. So it takes the function, input variables, linear inequalities, linear equalities and inequalities, lower and upper bound, non linear equalities and inequalities, options defined in this manner. These are some extra parameters which one can define to pass extra constaints/parameter on the objective function to get results accordingly.
- If there is no linear inequalities and equalities constraints then we replace them with empty square brackets, so in this case we are giving lower and upper bounds to search for results in constrained space. Then we took the results from the ga after multiple executions and storing them in respective variables.
- Now we plot the outputs from the ga as we done before, open a new figure window using the figure command. Now create two subplot using the subplot command, in the first we plot the outputs from the stalagmite function and the ga operation using the surf and plot3 command respectively. Then we plot the objective function value we get from ga in the second subplot. The axises and titles are provided wherever needed accordingly.
- The results are still not optimum as the ga might also taking the local minimas. In the ga inputs there is 'options' to control the ga execution or how the ga do its various process like selection, mutation etc. so to get better results. These options need to be defined before calling the ga function using the options command. It takes inputs as first ga then option1, value, option2, value. In this case we are giving a certain size of population to have more randomness, giving a bigger population size than default.
- This time we are providing more restricted upper and lower bounds and after execution of ga we store the output in respective variables. The plotting of the outputs is done as before, this time after each execution we are getting the same result and it seem to be optimal.
- Stalagmite Function :
function [f]= stalagmite(input)
f1=(sin(5.1*pi*input(1)+0.5))^6;
f2=exp(-4*log(2)*((input(1)-0.0667)^2/0.64));
f3=(sin(5.1*pi*input(2)+0.5))^6;
f4=exp(-4*log(2)*((input(2)-0.0667)^2/0.64));
f=-(f1*f2*f3*f4);
end
- Explanation :
In this function we use the stalagmite mathematical model to get output. There are two inputs each coordinate from the x-y plane respectively the x and y values. We are taking the output in the variable f.
- Output :
https://drive.google.com/drive/folders/13qWN-ahM47afCBOjj_8dsRO8-Sud_h85?usp=sharing
- Outcomes/Learnings :
- Genetic Algorithm : We can formalize a definition as such it is a search algorithm for finding a optimised solution for your application/problem. It tends to find the input value for which you will get the maximum value of the objective function. It is based on the Charles Darwin theory of evolution involving the process like selection, crossover and mutation.
So lets start with a problem which has vast solution space, in this if we try evaluating each solution that will take us years to find the optimal solution. Here what genetic algorithm does magic, instead of trying each solution what if we guess a bunch of solution randomly over which a function can tell how good the solution is. This is the main theme on which GA works but to control this random process there are lot constraints to find the best optimal solution, also we need to take care the constraints/process(selection,crossover,mutation) otherwise it can also lead some varied results or take very long time to execute. In Matlab the GA tool finds the minimum of a function rather than maximum so one to made some tweeks in the function itself to get the maximum value.
The algorithm picks a certain number of solutions randomly called population. There is function called fitness function which evaluate the fitness of the solution. After evaluating we take two solutions with different combination of fitness and break each of them into two parts and then replace each others half creating a new solution that we hope is better. This way we can create a whole bunch of new solution called the new generation. There is one more process called mutation which make small changes in the solution to create more randomness in the population/generation. After a certain number of generations the GA will stop and gives us the results.
- Multiple Iterations & Stopping GA : We execute the GA in a loop, to run it multiple times as to get mutiple results, it may have happen that the ga find the local minima so to run it multiple times we get verified that the ga has done its best under the given constraints. There is options command which helps to give much more constraints on how the GA execute, in that we have a option in which we can set the number of generation the GA must run. There are some stopping conditions for the algorithm, after which it will stall, but some times we need to evaluate more to get better results. So for that we have a command 'MaxStallGenerations' which sets the GA to run for minimum set number of generations, to evaluate the function more. To optimise the stalagmite function we can also use another global optimisation algorithm such as differential evolution, tabu search as these are also heuristics algorithms like genetic algorithms.
- References :
https://natureofcode.com/book/chapter-9-the-evolution-of-code/
https://in.mathworks.com/help/gads/ga.html#mw_a07ddf7d-face-4243-a656-f1528b592de0
https://youtube.com/playlist?list=PLRqwX-V7Uu6bJM3VgzjNV5YxVxUwzALHV
https://in.mathworks.com/help/gads/genetic-algorithm-options.html
https://in.mathworks.com/help/gads/setting-the-maximum-number-of-generations.html
https://en.m.wikipedia.org/wiki/Global_optimization
https://www.sciencedaily.com/releases/2004/12/041206193416.htm