Week 4.1 - Genetic Algorithm

% Program to find the global maxima of the function

x = linspace(0,0.6,150);
y = linspace(0,0.6,150);
 
% converting 1D array into 2D array
[xx,yy] = meshgrid(x,y);

num_cases = 50;

tic
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_h(input_vector);
    end
end

%surf(xx,yy,f)
figure(1)
surfc(xx,yy,f);
xlabel('x value')
ylabel('y value')
shading interp



% study 1 statistical behaviour

tic
for i = 1:num_cases
    [inputs,fopt(i)] = ga(@stalagmite_h,2);
    xopt(i) = inputs(1);
    yopt(i) = inputs(2);
end
study1_time = toc;
figure(1)
subplot(2,1,1)
hold on
surfc( xx,yy,f)
shading interp
plot3(xopt,yopt,fopt,'marker','o','markersize',5,'markerfacecolor','r')
title('Unbounded inputs')
subplot(2,1,2)
plot(fopt)
xlabel('Iterations')
ylabel('Function maximum')

% study 2 Bounded inputs, upper and lower bound
tic
for i = 1:num_cases
    [inputs,fopt(i)] = ga(@stalagmite_h,2,[],[],[],[],[0;0],[0.6;0.6]);
    xopt(i) = inputs(1);
    yopt(i) = inputs(2);
end
study2_time = toc;
figure(2)
subplot(2,1,1)
hold on
hold on
surfc(x,y,-f)
shading interp
plot3(xopt,yopt,-fopt,'marker','o','MarkerSize',5,'markerfacecolor','r')
title('Optimum solution')
subplot(2,1,2)
plot(-fopt)
xlabel('iterations')
ylabel('function maximum')


%study3 Bounded inputs with modified polpulation
options = optimoptions('ga');
options = optimoptions(options,'PopulationSize',300);
tic
for i = 1:num_cases
    [inputs,fopt(i)] = ga(@stalagmite_h,2,[],[],[],[],[0;0],[1;1],[],[],options);
    xopt(i) = inputs(1);
    yopt(i) = inputs(2);
end
study3_time = toc;
figure(3)
subplot(2,1,1)
hold on
surfc(x,y,-f)
shading interp
plot3(xopt,yopt,fopt,'marker','o','MarkerSize',5,'markerfacecolor','r')
title('Optimum solution')
subplot(2,1,2)
plot(-fopt)
xlabel('Iterations')
ylabel('Function maximum')
global_maxima = -fopt;




%stalagmite func

function [ f ] = stalagmite_h(input_vector)
x = input_vector(1);
y = input_vector(2);
f_x_1 = (sin(5.1*pi*x + 0.5)).^6;
f_y_1 = (sin(5.1*pi*y + 0.5)).^6;
f_x_2 = exp(-4.*log(2)*((x - 0.0667).^2)/(0.64));
f_y_2 =  exp(-4.*log(2)*((y - 0.0667).^2)/(0.64));
[f] = -(f_x_1.*f_y_1.*f_x_2.*f_y_2);