Optimization of a function using Genetic Algorithm in MATLAB

The project aims at optimization of a stalagmite function and finding the global maxima.

For this, we will be using the concept of genetic algorithm.

Following are the steps used to optimise the function:

  1. Code the function given in a separate function code. This is the function to be optimized.
  2. Define the search space space and create the mesh grid (2D Array).
  3. Evaluate the function for individual values of the mesh grid.
  4. Run the ga command for the fitness function (function to be optimized) to return x and y values, and maxima values of the function in a loop from 1 to total number of executions previously defined.
  5. Plot the same, first as a surface plot using the meshgrid values and function values, and make a second plot as a 3D plot of the maximas of the function returned by the ga funcion.
  6. Repeat the same two more times, once with bounds, and last with population defined.

Stalagmite Function:

function [f] = stalagmite(inputvector)
  
  ivx = inputvector(1);
  ivy = inputvector(2);
  
  f1x = (sin((5.1*pi*ivx)+0.5))^6;
  f1y = (sin((5.1*pi*ivy)+0.5))^6;
  f2x = exp(-4*log(2)*(((ivx-0.0667)^2)/0.64));
  f2y = exp(-4*log(2)*(((ivy-0.0667)^2)/0.64));
  
  f = -(f1x*f1y*f2x*f2y);
  %f = 1/(f+1);

end

 Primary Code:

clc
clear all
close all

% Search Space
x = linspace(0, 0.6, 150);
y = linspace(0, 0.6, 150);
numcases = 100

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

% Evaluate stalagmite function
for i = 1:length(xx)
  for j = 1:length(yy)
    inputvector(1) = xx(i, j);
    inputvector(2) = yy(i, j);
    f(i, j) = stalagmite(inputvector);
  end
end

%-------------------------------------------------------------------------------------
% Study 1, with no bounds or population
tic
for i = 1:numcases
    [inputs1, fval1(i)] = ga(@stalagmite, 2);
    xval1(i) = inputs1(1);
    yval1(i) = inputs1(2);
end

executiontime1 = toc

figure(1)
subplot(2,1,1)
hold on
surfc(xx, yy, -f)
shading interp %shading instead of grid lines on plot
plot3(xval1, yval1, -fval1, 'marker', 'o', 'markersize', 5, 'markerfacecolor', 'y')
xlabel('x-value');
ylabel('y-value');
zlabel('z-value');
title('Unbounded inputs/ no defined population')
subplot(2,1,2)
plot(-fval1)
xlabel('Iterations');
ylabel('Function Maximum');

%-------------------------------------------------------------------------------------
% Study 2, with bounds, no population
tic
for i = 1:numcases
    [inputs2, fval2(i)] = ga(@stalagmite, 2, [], [], [], [], [0 0], [1 1]);
    xval2(i) = inputs2(1);
    yval2(i) = inputs2(2);
end

executiontime2 = toc

figure(2)
subplot(2,1,1)
hold on
surfc(xx, yy, -f)
shading interp %shading instead of grid lines on plot
plot3(xval2, yval2, -fval2, 'marker', 'o', 'markersize', 5, 'markerfacecolor', 'y')
xlabel('x-value');
ylabel('y-value');
zlabel('z-value');
title('Bounded inputs/ no defined population')
subplot(2,1,2)
plot(-fval2)
xlabel('Iterations');
ylabel('Function Maximum');

%-------------------------------------------------------------------------------------
% Study 3, with bounds, and population defined

opts = optimoptions('ga');
opts.PopulationSize = 500;
% opts = optimoptions (opts, 'populationSize', 500);

tic
for i = 1:numcases
    [inputs3, fval3(i)] = ga(@stalagmite, 2, [], [], [], [], [0 0], [1 1], [], [], opts);
    xval3(i) = inputs3(1);
    yval3(i) = inputs3(2);
end

executiontime3 = toc

figure(3)
subplot(2,1,1)
hold on
surfc(xx, yy, -f)
shading interp %shading instead of grid lines on plot
plot3(xval3, yval3, -fval3, 'marker', 'o', 'markersize', 5, 'markerfacecolor', 'y')
xlabel('x-value');
ylabel('y-value');
zlabel('z-value');
title('Bounded inputs with defined population')
subplot(2,1,2)
plot(-fval3)
xlabel('Iterations');
ylabel('Function Maximum');

About Genetic Algorithm in MATLAB:

A genetic algorithm (GA) is a method for solving both constrained and unconstrained optimization problems based on a natural selection process that mimics biological evolution, influenced by Darwin's Theory of Evolution.

Sometimes the function is required to be optimised and maxima or minima found. Genetic Algorithm uses a random set of inputs for the calculation and uses the result as inputs for the next iteration, thereby reducing the error in the inputs. Over successive generations, the population evolves towards an optimal solution.

Syntax:

x = ga(fun,nvars)

Here, function is the fitness function, which is the function to be optimized. "nvars" is the number of variables in the fitness function.

x = ga(fun,nvars,A,b,[],[],lb,ub,nonlcon,IntCon,options)

Above is the complete syntax, where

  1. A and b are used such that the function calculates the minimum for x subject to linear inequalities A*x  <= b.
  2. [], [] are in reality Aeq and beq, which find x subject to linear equalities Aeq*x = beq. "[]" at any instance indicate no value.
  3. lb and ub ae lower and upper bounds of the design variables, so x is found in the range lb <= x < ub.
  4. "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.
  5. "IntCon" is such that if variables are listed in that, they take integer values.
  6. options:
    1. options = optimoptions('ga','Param1', value1, 'Param2', value2, ...);‚Äč

      Details the other options such as "population", which are a set of points in a design space. Initial population is generated randomly by default. Specifying the population increases accuracy and gives consistent results.

GA finding Minima:

GA by default finds the minima. However, here, we are required to find the maxima. That can be done in 2 ways:

  1. Using a minus sign inside the function in front of the final output statement
  2. f = 1/(function output + 1): this computes mean, and optimisation output shows the nean values.

Here, negative sign was used as it provides the real maxima.

f = -(f1x*f1y*f2x*f2y)

This however resulted in a flipped graph, so edits had to be made at the the time of plotting to flip it back.

...
surfc(xx, yy, -f)
...
plot3(xval1, yval1, -fval1, 'marker', 'o', 'markersize', 5, 'markerfacecolor', 'y')
...
...
plot(-fval1)
...

This way, maxima was plotted with respect to iterations.

Outputs:

  1. No bounds or population:

  2. Bounded but without population

  3. Bounded, and with defined population

As we can see, the more defined our input constraints get, like specified bounds and poulation, the more accurate our result of the maxima is, showing fewer troughs and hence, few variations, essentially indicating a more precise results.

Having bounds defined and population as well also helps in repetitive results over multiple executions.

As mentioned by Mathworks, The default population size used by ga is 50 when the number of decision variables is less than 5 and 200 otherwise. This size can be a poor one for some problems; a smaller population size can be sufficient for smaller problems. We have selected 500 here for accuracy.

Errors:

  1. First attempt at getting maxima:

    Took the minus and put it in function output

  2. Matrix operation error: Dots misplaced in formula

    Corrected the dots.

  3. Study 3 error: syntax error

    Corrected syntax based on what's on Mathworks site for A, B, Aeq, Beq, IntCon and nonlcon.

References:

  1. https://in.mathworks.com/matlabcentral/answers/52240-how-to-find-max-fuction-with-genetic-algorithm
  2. Global vs. Local Minima Using ga - MATLAB & Simulink - MathWorks India
  3. Genetic Algorithm - MATLAB & Simulink
  4. Genetic Algorithm Options - MATLAB & Simulink Example - MathWorks India
  5. Genetic Algorithm Options - MATLAB & Simulink - MathWorks India
  6. Find minimum of function using genetic algorithm - MATLAB ga - MathWorks India

Projects by Arjun Bhat

This project explores the simulation of the motion of a pendulum for the equation of damped vibration of a pendulum. Below is the equation (taken from "Challenges" page): The above equation is solved using ODE function, and then odeint is used to get an ar Read more

Project involves plotting an Otto Cycle graph in MATLAB The goal is to plot all strokes for a petrol operating by the otto cycle from intake, compression, power to exhaust, including both adiabatic curves. Description: The Otto Cycle is the idealised cycle that descr Read more

Objective: The project involves simulating a 2-joint robotic arm.   Code and Description: import math import matplotlib.pyplot as plt import numpy as nm # Lengths l1 = 1 l2 = 0.5 # Angles ang1 = nm.radians(nm.linspace(0, 90, 19)) ang2 = nm.radians Read more

Flow over bicycle in Python
Arjun Bhat · 2019-11-18 16:09:30

Objective: To calculate the drag force for various drag coefficient and velocity values for a certain geometry.   Introduction: Drag force is the force exerted by the air as the vehicle moves through it. The force exerted is directly proportional to the square Read more

Aim: To parse the NASA Thermodynamc Properties data file, calculate thermodynamic properties and plot the characteristics. Introduction: NASA's thermodynamic properties file consists of multiple species of gases, their respective operating temperature ranges, and Read more

Curve Fitting using MATLAB
Arjun Bhat · 2019-08-19 20:20:53

Project aims at finding the best fit curve for a given data and finding the error between the original and said curve. The data being used is that of Specific heat (Cp) vs Temeperature. Plot will be made to show variation of Cp with Temp first.  Code: % Data cp Read more

This project explores the simulation of the motion of a pendulum for the equation of damped vibration of a pendulum. Below is the equation (taken from "Challenges" page): An attempt is made to solve the above equation, and then call it using ode45 to get an array of Read more

Project involves plotting an Otto Cycle graph in MATLAB The goal is to plot all strokes for a petrol operating by the otto cycle from intake, compression, power to exhaust, including both adiabatic curves. Description: The Otto Cycle is the idealised cycle that descr Read more

Code and Description: The project involves simulating a 2-joint robotic arm. % Lengths l1 = 1; l2 = 0.5; % Angles ang1 = linspace (0, 90, 20); ang2 = linspace (0, 90, 20); count = 1; % loop for i = 1:length(ang1) ANG1 = ang1(i); for j = 1:length(ang2 Read more


Loading...

The End