Challenge

Finding Global Maxima of the Stalagmite function by optimizing it using Genetic Algorithm in MATLAB. (Matlab/Challenge-Week4.1-Genetic Algorithm)

SSachin Kandpal22 Feb 2021

AIM : To find Global Maxima of the Stalagmite function by optimizing it using Genetic Algorithm in MATLAB.

OBJECTIVES :

Write a code in MATLAB to optimise the stalagmite function and find the global maxima of the function.
Clearly explain the concept of genetic algorithm in your own words and also explain the syntax for ga in MATLAB in your report.

THEORY :

Genetic Algorithm (GA) is a search-based optimization technique based on the principles of Genetics and Natural Selection. It is frequently used to find optimal or near-optimal solutions to difficult problems which otherwise would take a lifetime to solve. It is frequently used to solve optimization problems, in research, and in machine learning.

The genetic algorithm is a method for solving both constrained and unconstrained optimization problems that are based on natural selection, the process that drives biological evolution. The genetic algorithm repeatedly modifies a population of individual solutions. At each step, the genetic algorithm selects individuals at random from the current population to be parents and uses them to produce the children for the next generation. Over successive generations, the population "evolves" toward an optimal solution. You can apply the genetic algorithm to solve a variety of optimization problems that are not well suited for standard optimization algorithms, including problems in which the objective function is discontinuous, nondifferentiable, stochastic, or highly nonlinear. The genetic algorithm can address problems of mixed-integer programming, where some components are restricted to be integer-valued.

The genetic algorithm uses three main types of rules at each step to create the next generation from the current population:

Selection rules select the individuals, called parents, that contribute to the population at the next generation.
Crossover rules combine two parents to form children for the next generation.
Mutation rules apply random changes to individual parents to form children.

Genetic Algorithm Terminology -

Fitness Functions
The fitness function is the function you want to optimize. For standard optimization algorithms, this is known as the objective function. The toolbox software tries to find the minimum of the fitness function.
Write the fitness function as a file or anonymous function, and pass it as a function handle input argument to the main genetic algorithm function.
Populations and Generations
A population is an array of individuals. For example, if the size of the population is 100 and the number of variables in the fitness function is 3, you represent the population by a 100-by-3 matrix. The same individual can appear more than once in the population. For example, the individual (2, -3, 1) can appear in more than one row of the array. At each iteration, the genetic algorithm performs a series of computations on the current population to produce a new population. Each successive population is called a new generation.
Diversity
Diversity refers to the average distance between individuals in a population. A population has high diversity if the average distance is large; otherwise, it has low diversity. In the following figure, the population on the left has high diversity, while the population on the right has low diversity.
Diversity is essential to the genetic algorithm because it enables the algorithm to search a larger region of the space.
Individuals
An individual is any point to which you can apply the fitness function. The value of the fitness function for an individual is its score.
Fitness Values and Best Fitness Values
The fitness value of an individual is the value of the fitness function for that individual. Because the toolbox software finds the minimum of the fitness function, the best fitness value for a population is the smallest fitness value for any individual in the population.
Parents and Children
To create the next generation, the genetic algorithm selects certain individuals in the current population, called parents, and uses them to create individuals in the next generation, called children. Typically, the algorithm is more likely to select parents that have better fitness values.

Working Outline of the Genetic Algorithm :

The following outline summarizes how the genetic algorithm works:

The algorithm begins by creating a random initial population.
The algorithm then creates a sequence of new populations. At each step, the algorithm uses the individuals in the current generation to create the next population. To create the new population, the algorithm performs the following steps:
1. Scores each member of the current population by computing its fitness value. These values are called raw fitness scores.
2. Scales the raw fitness scores to convert them into a more usable range of values. These scaled values are called expectation values.
3. Selects members, called parents, based on their expectation.
4. Some of the individuals in the current population that have lower fitness are chosen as elite. These elite individuals are passed to the next population.
5. Produces children from their parents. Children are produced either by making random changes to a single parent—mutation—or by combining the vector entries of a pair of parents—crossover.
6. Replaces the current population with the children to form the next generation.
The algorithm stops when one of the stopping criteria is met.

Genetic Algorithm Syntax in MATLAB

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. When there are integer constraints, ga does not accept linear or nonlinear equality constraints, only inequality constraints.
x = ga(problem) finds the minimum for problem, a structure described in problem.
[x,fval] = ga(___), for any previous input arguments, also returns fval, the value of the fitness function at x.
[x,fval,exitflag,output] = 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(___) also returns a matrix population, whose rows are the final population, and a vector scores, the scores of the final population.

Stopping criteria of Genetic Algorithm -

The termination condition of a Genetic Algorithm is important in determining when a GA run will end. It has been observed that initially, the GA progresses very fast with better solutions coming in every few iterations, but this tends to saturate in the later stages where the improvements are very small. We usually want a termination condition such that our solution is close to the optimal, at the end of the run.

Usually, we keep one of the following termination conditions −

When there has been no improvement in the population for X iterations.
When we reach an absolute number of generations.
When the objective function value has reached a certain pre-defined value.

For example, in a genetic algorithm, we keep a counter which keeps track of the generations for which there has been no improvement in the population. Initially, we set this counter to zero. Each time we don’t generate off-springs which are better than the individuals in the population, we increment the counter.

However, if the fitness of any of the off-springs is better, then we reset the counter to zero. The algorithm terminates when the counter reaches a predetermined value.

Stalagmite and Stalactite

A stalagmite is a type of rock formation that rises from the floor of a cave due to the accumulation of material deposited on the floor from ceiling drippings.

A stalactite is a mineral formation that hangs from the ceiling of caves, hot springs, or man-made structures such as bridges and mines. Any material that is soluble and that can be deposited as a colloid, or is in suspension, or is capable of being melted, may form a stalactite.

Stalagmite Function

A stalagmite function generally refers to the function with which we are generating values over a stalagmite kind of surface for finding out the local and global maxima and minima of that particular given function.

`f(x,y) = f_(1x)*f_(2x)*f_(1y)*f_(2y)`

Where,

`f_(1x) = [sin(5.1*pi*x+0.5)]^6`

`f_(1y) = [sin(5.1*pi*y+0.5)]^6`

`f_(2x) = e^[-4*ln(2)*(x-0.0667)^2/0.64]`

`f_(2y) = e^[-4*ln(2)*(y-0.0667)^2/0.64]`

PROGRAM :

1. Main Function -

clear all
close all
clc

% defining search space
x = linspace(0,0.6,150);
y = linspace(0,0.6,150);

num_cases = 50;
f = get_stalagmite(x,y);

%...............study of statistical behaviour...............

% STUDY 1  - In this the inputs of GA are unbounded i.e. GA can take input as any sampling from global range(has no range limit) 

tic;
for i = 1:num_cases
    [inputs, f_opt(i)] = ga(@stalagmite, 2);         % fopt is the the optimized value after optimization of the genetic algorithm.  
    x_opt(i) = inputs(1);                            % ga minimize a function.. gives function minima... here Fopt is the minimum value           
    y_opt(i) = inputs(2);                            % and xopt, yopt are the inputs at which we get the function minima. 
end                                                  % ga checks for minima by several inputs.. and returns optimiun value and its inputs.
study1_time = toc

figure(1);
    subplot(2,1,1);
    hold on;
        surfc(x,y,-f);                               % Since we used negative sign in stalagamite function to maximize using GA.. Now plotting results a  
        shading interp;
        xlabel('X-Axis');     ylabel('Y-Axis');    zlabel('Z-Axis');
        title('STUDY 1  - Inputs of GA are unbounded');
        plot3(x_opt, y_opt, -f_opt,'marker','o','markersize',5,'markerfacecolor','k');
    subplot(2,1,2);   
        plot(-f_opt);
        xlabel('Iterations');     ylabel('Function Maxima');
        title('Iteration vs Function Maxima');

% STUDY 2 -In this the inputs of GA are bounded i.e. GA take inputs or samplings from local user-specified range

tic;
for i = 1:num_cases
    [inputs, f_opt(i)] = ga(@stalagmite, 2,[],[],[],[],[0;0],[1;1]);
    x_opt(i) = inputs(1);
    y_opt(i) = inputs(2);
end
study2_time = toc

figure(2);
    subplot(2,1,1);
    hold on;
        surfc(x,y,-f);   
        shading interp;
        xlabel('X-Axis');     ylabel('Y-Axis');    zlabel('Z-Axis');
        title('STUDY 2  - Inputs of GA are Bounded');
        plot3(x_opt, y_opt, -f_opt,'marker','o','markersize',5,'markerfacecolor','k');
    subplot(2,1,2);   
        plot(-f_opt);
        xlabel('Iterations');     ylabel('Function Maxima');
        title('Iteration vs Function Maxima');
        
% STUDY 3 - here in this the inputs of GA are bounded same as study 2, but here 
%           we are raising the population size(GA as default have population size of 50) inorder to increase optimization accuracy, 
%           as increasing population leads to more cross-over and variety for selection of the optimum value.

options = optimoptions('ga');
options = optimoptions(options,'PopulationSize',250);

tic;
for i = 1:num_cases
    [inputs, f_opt(i)] = ga(@stalagmite, 2,[],[],[],[], [0;0], [1;1],[],[], options);
    x_opt(i) = inputs(1);
    y_opt(i) = inputs(2);
end
study3_time = toc

figure(3);
    subplot(2,1,1);
    hold on;
        surfc(x,y,-f);   
        shading interp;
        xlabel('X-Axis');     ylabel('Y-Axis');    zlabel('Z-Axis');
        title('STUDY 3 - Bounded Inputs with Increased Population');
        plot3(x_opt, y_opt, -f_opt,'marker','o','markersize',5,'markerfacecolor','k');
    subplot(2,1,2);   
        plot(-f_opt);
        xlabel('Iterations');     ylabel('Function Maxima');
        title('Iteration vs Function Maxima');
        
% STUDY 4 - its same as study 3, but here were reducing bounded region as we know that optimium is in lesser range

options = optimoptions('ga');
options = optimoptions(options,'PopulationSize',250);

tic;
for i = 1:num_cases
    [inputs, f_opt(i)] = ga(@stalagmite, 2,[],[],[],[], [0;0], [0.5;0.5],[],[], options);
    x_opt(i) = inputs(1);
    y_opt(i) = inputs(2);
end
study4_time = toc

figure(4);
    subplot(2,1,1);
    hold on;
        surfc(x,y,-f);   
        shading interp;
        xlabel('X-Axis');     ylabel('Y-Axis');    zlabel('Z-Axis');
        title('STUDY 4 - Reducing Input bounded region Further');
        plot3(x_opt, y_opt, -f_opt,'marker','o','markersize',5,'markerfacecolor','k');
    subplot(2,1,2);   
        plot(-f_opt);
        xlabel('Iterations');     ylabel('Function Maxima');
        title('Iteration vs Function Maxima');

2. get_stalagmite Funtion -

function f = get_stalagamite(x,y)
   
    % creating 2 Dimensional mesh array
    [xx, yy] = meshgrid(x,y);

    % Evaluating 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
    
end

3. Stalagmite fitness function which is gonna be used for optimization from Genetic Algorithm -

function [f] = stalagmite(i_vec)

    x = i_vec(1);
    y = i_vec(2);

    pi = 22/7;

    f1x = (sin((5.1*pi*x)+0.5))^6;
    f1y = (sin((5.1*pi*y)+0.5))^6;

    f2x = exp( -4*log(2)*((x-0.0667)^2 / 0.64) );
    f2y = exp( -4*log(2)*((y-0.0667)^2 / 0.64) );

    f = f1x * f1y * f2x * f2y;

    f = f*(-1);                 % Since GA in matlab only minimizes a function therefore, 
                                % inorder to find maxima we change the to negative form
end

OUTPUT :

1) Study 1 -