Using genetic algorithm in MATLAB to optimize a mathematical model

A MATLAB program has been written to optimize a Mathematical function, called the stalagmite function, and find it\'s global maxima.
This function helps us generate a number of local minima and maxima. The optimization in this program is done using the Genetic Algorithm inbuilt in MATLAB.

What is Genetic Algorithm?

Genetic Algorithms(GAs) are adaptive heuristic search algorithms that belong to the larger part of evolutionary algorithms. Genetic algorithms are based on the ideas of natural selection and genetics. These are intelligent exploitation of random search provided with historical data to direct the search into the region of better performance in solution space. They are commonly used to generate high-quality solutions for optimization problems and search problems,

Genetic algorithms simulate the process of natural selection which means those species who can adapt to changes in their environment are able to survive and reproduce and go to next generation. In simple words, they simulate “survival of the fittest” among individual of consecutive generation for solving a problem. Each generation consist of a population of individuals and each individual represents a point in search space and possible solution. Each individual is represented as a string of character/integer/float/bits. This string is analogous to the Chromosome.

Search space

The population of individuals are maintained within search space. Each individual represent a solution in search space for given problem. Each individual is coded as a finite length vector (analogous to chromosome) of components. These variable components are analogous to Genes. Thus a chromosome (individual) is composed of several genes (variable components).

Fitness Score

A Fitness Score is given to each individual which shows the ability of an individual to “compete”. The individual having optimal fitness score (or near optimal) are sought.

Flow Chart

Steps Involved:

Step 1: Initialize the population. A population of n individuals (chromosomes) are selected, either randomly or based on the boundary conditions. This defines the search space. (In MATLAB, by default, the initial population size is 50 in the ga function)

Step 2: A fitness criteria is developed and individuals (chromosomes) with better fitness scores are selected. If these chromosomes satisfy the stop condition, the operation is stopped and an Output is recieved. If not, the chromosomes are selected at random in Roulette wheel fashion, with chromosomes with a higher fitness score having a greater chance of being selected i.e: Greater sector area on the Roulette wheel would be given to chromosomes with better fitness. The idea is to give preference to the individuals with good fitness scores and allow them to pass there genes to the successive generations.
After the selection the chromosomes are placed in an order.

Step 3: Now from these selected chromosomes, a crossover operation is performed between 2 chromosomes, in a single or multiple point style, between any of selected chromosomes. This operation represents mating between individuals. Two individuals are selected using the selection operator and crossover sites are chosen randomly. Then the genes at these crossover sites are exchanged thus creating a completely new individual (offspring). Example:

Step 4: From these newly formed offspring chromosomes, mutation is performed at random genes and in random chromosomes. The key idea is to insert random genes in offspring to maintain the diversity in population to avoid the premature convergence i.e: to ensure that the search algorithm does not necessarily get stuck at the local maxima. Example of Mutation:

Step 5: After Mutation, the fitness of the resultant chromosomes evaluated and if they meet the convergence criterion, they are given ad the output or else the cycle repeats.

To see solved examples Genetic Algorithm, search for Knapsack or MAXONE problems.

Syntax

1.) x = ga(@function,Number of variables, A, b, Aeq, Beq, lb, ub, noncolon, options)

• 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)
  Finds a local minimum x to fun, subject to the linear inequalities.
  
  
• 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 or gaoptimset. We can change parameters like population size and other using this command.
  
  
• 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.
  
  
• Further modifications or specifications can be set using the Optimization Toolbox.

In this program, we use the ga function in a situation where [inputs,fopt(n)] = ga(@stalagmite,2,[],[],[],[],[0;0],[1;1],[],[],options) , where the function to be minimized, the number of vaariables, specified as positive integers, and the upper and lower bound coordinates, along with an increase in population size, using the options command are used.

2.) subplot(m,n,p)

• subplot(m,n,p) divides the current figure into an m-by-n grid and creates axes in the position specified by p. MATLAB numbers subplot positions by row. The first subplot is the first column of the first row, the second subplot is the second column of the first row, and so on. If axes exist in the specified position, then this command makes the axes the current axes.

3.) plot3(X1,Y1,Z1)

• The plot3 function displays a three-dimensional plot of a set of data points.
• plot3(X1,Y1,Z1), where X1, Y1, Z1 are vectors or matrices, plots one or more lines in three-dimensional space through the points whose coordinates are the elements of X1, Y1, and Z1

A Total of 3 statistical studies are carried out and plots are made for the maximum value of the function vs no. of iterations

The program code is given below:

clear all
close all
clc

%Defining our search space
x = linspace(0,0.6,150);
y = linspace(0,0.6,150);
num_cases = 50; %Number of Iterations

% Creating a 2D Mesh matrix
[xx, yy] = meshgrid(x,y);

% Evaluating Stalagmite Function at all values of x and y
for k = 1:length(xx)
    for j = 1:length(yy)
        
        input_vector(1) = xx(k,j);
        input_vector(2) = yy(k,j);
        
        f(k,j) = stalagmite(input_vector);
    end
end

tic
% Study 1 - Statistical Behaviour without boundary conditions
for i = 1:num_cases
    [inputs, fopt(i)] = ga(@stalagmite,2);
    
    xopt(i) = inputs(1);
    yopt(i) = inputs(2); 
    
end

study1_time = toc;

figure(1)
subplot(2,1,1)
hold on
surfc(x,y,f)
xlabel(\'x value\')
ylabel(\'y value\')
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 minimum\')

%in this study, the problem is that the ga takes random values and does not
%stick to the boundaries of x and y values, we will fix this in study 2

% Study 2 - Statistical Behaviour - with upper and lower bounds
tic
for m = 1:num_cases
    [inputs,fopt(m)] = ga(@stalagmite,2,[],[],[],[],[0;0],[1;1]); 
    % Bounded from 0 to 1 for x and y
    xopt(m) = inputs(1);
    yopt(m) = inputs(2);
end
study2_time = toc;

figure(2)
subplot(2,1,1)
hold on
surfc(x,y,f)
xlabel(\'x value\')
ylabel(\'y value\')
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 minimum\')

% Study 3 - Increasing GA iterations by increasing intial population size
% to more than 50

options = gaoptimset(\'ga\')
options = gaoptimset(options,\'PopulationSize\',170);

tic
for n = 1:num_cases
    [inputs,fopt(n)] = ga(@stalagmite,2,[],[],[],[],[0;0],[1;1],[],[],options);
    xopt(n) = inputs(1);
    yopt(n) = inputs(2);
end

study3_time = toc;

figure(3)
subplot(2,1,1)
hold on
surfc(x,y,f)
xlabel(\'x value\')
ylabel(\'y value\')
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 minimum\')

figure(4)
surfc(xx,yy,f)
xlabel(\'x value\')
ylabel(\'y value\')

shading interp

[inputs,fval] = ga(@stalagmite,2,[],[],[],[],[0;0],[1;1],[],[],options)

study1_time,study2_time,study3_time

• x and y co-ordinates define our searchspace.
  
  
• num_cases represents the number of iterations we want to run the algorithm for.
  
  
• Two for loops that are nested are used to access the coordinates of the 2-dimensional grid. Inner for loop traces column elements and outer for loop traces the row elements of the coordinate matrices(xx and yy). The traced x and y coordinates are stored as elements of the \'input_vector\'.
  
  
• The Stalagmite function takes the \'input_vector\' as input and calculates the corresponding z value.
  
  
• We obtain the contour using \'surfc\' command. \'shading interp\' removes the superimposed black mesh lines.
  
  
• For the first study, Study 1, we run the genetic algorithm for stalagmite function without using any boundary conditions and make plots for the result. These plots represent the change in the function maximum value along with the iteration number.
  
  
• For the second study, Study 2, we run the genetic algorithm for stalagmite function specifying the lower boundary as 0 in a matrix of doubles, and the upper boundary as 1 in a matrix of doubles. The result plots represent the change in the function maximum value along with the iteration number
  
  
• For the second study, Study 3, we run the genetic algorithm for stalagmite function specifying the lower boundary as 0 in a matrix of doubles, and the upper boundary as 1 in a matrix of doubles along with increasing the intial population size from the default 50 to 170 using the gaoptiset function.
  
  
• The plots show the sequence in which the genetic algorithm moved through elements to find the minimum of the function. However the stalagmite function is coded in such a way that the minimum value found by the \'ga\' function gives the actual maximum value of the stalagmite function.
  
  
• The output is time taken by each study and the x and y coordinates at which function value is maximum.
  
  
• The tic and toc commands are used to measure the time taken by a particular operation

The Stalagmite function:

The eqn for the function is given by:

The two variables to be input are x and y. The code for the function in MATLAB is given below:

function [f] = stalagmite(input_vector)

x = input_vector(1);
y = input_vector(2);
f1_x = (sin(5.1*pi*x + 0.5)).^6;
f1_y = (sin(5.1*pi*y + 0.5)).^6;
f2_x = exp(-4*(log(2))*(((x-0.06670).^2)./0.64));
f2_y = exp(-4*(log(2))*(((y-0.06670).^2)./0.64));   

g = f1_x*f2_x*f1_y*f2_y;

f = 1/(1+g); % Both are correct but this one gives more accurate results and faster computation
%f = -g;
end

The input vector is a two element array from the main code which provides the x and y input values for the function. Since the MATLAB \'ga\' function find the minimum value of a function, i have taken the inverse of the actual function g represented as f = 1/1+g, so that the minimum value of function f gives the maximum value of function g.

Error Faced:

The function optimoptions() does not work in my version of MATLAB. Instead, i have used the function gaoptimset().

Results:

A screenshot of the result is shown below:

inputs gives the x and y coordinates at which the function value is maximum. The function maximum value at the x and y coordinates given above is:

Plots:

Surface plot with x,y and the function value as z. This plot shows the local maxima and the global maxima. This plot has also been shown in the 3 studies given below so as to compare the progression of the solution generated by the genetic algorithm in MATLAB to the actual answer required.

Study 1:

Since the initial population in study 1 is unbounded, the results of this function minimization of f, are extremely varied and are between 0.5 and 1 (Function f = 1/1+g, where g is original function)

Study 2:

Since the initial population in study 2 is bounded between 0 and 1, the results are not as varied as study 1. The ga function results move between 0.5 and 0.7, sometimes getting stuck at local maxima.(Function f = 1/1+g, where g is original function)

Study 3:

Since the initial population of study 3 is increased from 50 to 170, it does not get at local maxima and maintains a constant resulting solution between 0.5000 and 0.5003.

Even though study 3 gives the most precise and accurate results, it takes the most amount of time to compute.

Comparsion between time taken to run each of the studies depending on how the inverse for the stalagmite function is taken

Function Inverse Type	Study 1	Study 2	Study 3
f = 1/(1+g)	11.6960 s	14.3320 s	33.6266 s
11.7309 s	14.3680 s	34.3845 s
11.9081 s	14.5991 s	32.5586 s
f = - g	13.4479 s	14.9242 s	39.145 s
12.5500 s	14.6420 s	43.1763 s
11.0490 s	14.4825 s	39.4547 s

The above table shows that, when the initial population size increases, the inverese function type f = 1/(1+g), where g is the original stalagmite function, computes the minimum value of the function faster than when f = -g.