Global optimization of stalagmite function using Genetic algorithm
AIM:To optimise the stalagmite function and find the global maxima of the function.
THEORY:
Global Optimization:Global optimization is a branch of applied mathematics and numerical analysis that attempts to find the global minima or maxima of a function or a set of functions on a given set.
Some variables used in global optimization are,
- LocalMaximum: A (local) maximum ˆxl ∈ X of one (objective) function
f : X 7→ R is an input element with f(ˆxl) ≥ f(x) for all x neighboring ˆxl. - Local Minimum:A (local) minimum ˇxl ∈ X of one (objective) function
f : X 7→ R is an input element with f(ˇxl) ≤ f(x) for all x neighboring ˇxl. - Local Optimum: A (local) optimum x⋆l ∈ X of one (objective) function
f : X 7→ R is either a local maximum or a local minimum. - GlobalMaximum: A global maximum ˆx ∈ x of one (objective) function
f : X 7→ R is an input element with f(ˆx) ≥ f(x) ∀x ∈ X. - Global Minimum: A global minimum ˇx ∈ X of one (objective) func-
tion f : X 7→ R is an input element with f(ˇx) ≤ f(x) ∀x ∈ X.
Figure below illustrates the variables ,
Genetic algorithm: Genetic algorithm developed by Goldberg was inspired by Darwin's theory of evolution which states that the survival of an organism is affected by rule "the strongest species that survives". Darwin also stated that the survival of an organism can be maintained through the process of reproduction, crossover and mutation.
In the genetic algorithm process is as follows :
Step 1. Determine the number of chromosomes, generation, and mutation rate and crossover rate value
Step 2. Generate chromosome-chromosome number of the population, and the initialization value of the genes chromosome-chromosome with a random value
Step 3. Process steps 4-7 until the number of generations is met
Step 4. Evaluation of fitness value of chromosomes by calculating objective function
Step 5. Chromosomes selection
Step 6. Crossover
Step 7. Mutation
Step 8. Solution (Best Chromosomes)
The flowchart of algorithm can be seen in Figure
In MATLAB,
The syntax for Genaral algorithm is,
x = ga(fun,nvars,A,b)
x = ga(fun,nvars,A,b,Aeq, beq)
x = ga(fun,nvars,A,b,Aeq,beq,lb ,ub)
x = ga(fun,nvars,A,b,Aeq,beq,lb ,ub,nonlcon)
x = ga(fun,nvars,A,b,Aeq,beq,lb,ub ,nonlcon,options)
x = ga(fun,nvars,A,b,[],[],lb,ub ,nonlcon,IntCon)
x = ga(fun,nvars,A,b,[],[],lb,ub ,nonlcon,IntCon,options)
where,
(fun):Objective function, specified as a function handle or function name.
ex:fun=@x(x.^2);
nvars:Number of variables, specified as a positive integer. The solver passes row vectors of length nvars to fun.
A:Linear inequality constraints, specified as a real matrix. A is an M-by-nvars matrix, where M is the number of inequalities.
A encodes the M linear inequalities
A*x <= b,
where x is the column vector of nvars variables x(:), and b is a column vector with M elements.
b:Linear inequality constraints, specified as a real vector. b is an M-element vector related to the A matrix. If you pass b as a row vector, solvers internally convert bto the column vector b(:).
b encodes the M linear inequalities
A*x <= b,
where x is the column vector of N variables x(:), and A is a matrix of size M-by-N.
Aeq:Linear equality constraints, specified as a real matrix. Aeq is an Me-by-nvars matrix, where Me is the number of equalities.
Aeq encodes the Me linear equalities
Aeq*x = beq,
where x is the column vector of N variables x(:), and beq is a column vector with Me elements.
beq:Linear equality constraints, specified as a real vector. beq is an Me-element vector related to the Aeq matrix. If you pass beq as a row vector, solvers internally convert beq to the column vector beq(:).
beq encodes the Me linear equalities
Aeq*x = beq,
where x is the column vector of N variables x(:), and Aeq is a matrix of size Meq-by-N.
lb:Lower bounds, specified as a real vector or array of doubles. lb represents the lower bounds element-wise in lb ≤ x ≤ ub.
Internally, ga converts an array lb to the vector lb(:).
ub:Upper bounds, specified as a real vector or array of doubles. ub represents the upper bounds element-wise in lb ≤ x ≤ ub.
Internally, ga converts an array ub to the vector ub(:).
noncolon:Nonlinear constraints, specified as a function handle or function name. nonlcon is a function that accepts a vector or array x and returns two arrays, c(x) andceq(x).
Optimization options, specified as the output of optimoptions or a structure.
The following outline summarizes how the genetic algorithm works in MATLAB:
-
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:
-
Scores each member of the current population by computing its fitness value. These values are called the raw fitness scores.
-
Scales the raw fitness scores to convert them into a more usable range of values. These scaled values are called expectation values.
-
Selects members, called parents, based on their expectation.
-
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.
-
Produces children from the 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.
-
Replaces the current population with the children to form the next generation.
-
-
The algorithm stops when one of the stopping criteria is met.
The genetic algorithm uses the following conditions to determine when to stop:
- Generations — The algorithm stops when the number of generations reaches the value of Generations.
- Time limit — The algorithm stops after running for an amount of time in seconds equal to Time limit.
- Fitness limit — The algorithm stops when the value of the fitness function for the best point in the current population is less than or equal to Fitness limit.
- Stall generations — The algorithm stops when the average relative change in the fitness function value over Stall generations is less than Function tolerance.
- Stall time limit — The algorithm stops if there is no improvement in the objective function during an interval of time in seconds equal to Stall time limit.
- Stall test — The stall condition is either average change or geometric weighted. For geometric weighted, the weighting function is 1/2n, where n is the number of generations prior to the current. Both stall conditions apply to the relative change in the fitness function over Stall generations.
- Function tolerance — The algorithm runs until the average relative change in the fitness function value over Stall generations is less than Function tolerance.
- Constraint tolerance — The Constraint tolerance is not used as stopping criterion. It is used to determine the feasibility with respect to nonlinear constraints. Also, max(sqrt(eps),Const aintTolerance) determines feasibility with respect to linear constraints.
Stalagmite Function:
Code explanation:
We initialize two random vaiables x&y through which we define a matrix.
x=linspace(0,0.6,150);
y=linspace(0,0.6,150);
[xx,yy]=meshgrid(x,y);To get the stalagmite functions, we use nested for loop and we compute the value f,
then we compute the local maximas to get get optimam values,xopt,yopt,fopt
%Get the 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
endthen we study the cases using different conditions to get the global maxima value at unbounded,bounded& bounded with population size respectively
% case study 1 with unbounded conditions
tic
for i=1:no_cases
[inputs,fopt(i)]=ga(@stalagmite,2);
xopt(i)=inputs(1);
yopt(i)=inputs(2);
end
study_time1=toc;
figure(1)
subplot(2,1,1)
surfc(x,y,-f)
shading interp
hold on
plot3(xopt,yopt,-fopt,'marker','o','markerSize',5)
title('Unbounded inputs')
xlabel('x value')
ylabel('y value')
zlabel('f value')
subplot(2,1,2)
plot(-fopt)
xlabel('iterations')
ylabel('function minimum')%Case 2 study with bounded conditions
tic
for i=1:no_cases
[inputs,fopt(i)]=ga(@stalagmite,2,[],[],[],[],[0;0],[0.6;0.6]);
xopt(i)=inputs(1);
yopt(i)=inputs(2);
end
study_time2=toc;
figure(2)
subplot(2,1,1)
surfc(x,y,-f)
shading interp
hold on
plot3(xopt,yopt,-fopt,'marker','o','markerSize',5)
title('bounded inputs')
xlabel('x value')
ylabel('y value')
zlabel('f value')
subplot(2,1,2)
plot(-fopt)
xlabel('iterations')
ylabel('function minimum')
%Case 3 study with bounded conditions
tic
options=optimoptions('ga');
options=optimoptions(options,'populationSize',300);
for i=1:no_cases
[inputs,fopt(i)]=ga(@stalagmite,2,[],[],[],[],[0;0],[0.6;0.6],[],[],options);
xopt(i)=inputs(1);
yopt(i)=inputs(2);
end
study_time3=toc;
figure(3)
subplot(2,1,1)
surfc(x,y,-f)
shading interp
hold on
plot3(xopt,yopt,-fopt,'marker','o','markerSize',5,'markerfacecolor','r')
title('Bounded inputs with population size')
xlabel('x value')
ylabel('y value')
zlabel('f value')
subplot(2,1,2)
plot(-fopt)
xlabel('iterations')
ylabel('function minimum')
fprintf('Study time for Study one is:%f \n',study_time1);
fprintf('Study time for Study two is:%f \n' ,study_time2);
fprintf('Study time for Study three is:%f \n',study_time3);
FULL CODE:
%Global optimization of a stalagmite function using genetic algorithm
clear all
close all
clc
% initialization
x=linspace(0,0.6,150);
y=linspace(0,0.6,150);
% No of cases to be studied
no_cases=50;
[xx,yy]=meshgrid(x,y);
%Get the 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
% case study 1 with unbounded conditions
tic
for i=1:no_cases
[inputs,fopt(i)]=ga(@stalagmite,2);
xopt(i)=inputs(1);
yopt(i)=inputs(2);
end
study_time1=toc;
figure(1)
subplot(2,1,1)
surfc(x,y,-f)
shading interp
hold on
plot3(xopt,yopt,-fopt,'marker','o','markerSize',5)
title('Unbounded inputs')
xlabel('x value')
ylabel('y value')
zlabel('f value')
subplot(2,1,2)
plot(-fopt)
xlabel('iterations')
ylabel('function minimum')
%Case 2 study with bounded conditions
tic
for i=1:no_cases
[inputs,fopt(i)]=ga(@stalagmite,2,[],[],[],[],[0;0],[0.6;0.6]);
xopt(i)=inputs(1);
yopt(i)=inputs(2);
end
study_time2=toc;
figure(2)
subplot(2,1,1)
surfc(x,y,-f)
shading interp
hold on
plot3(xopt,yopt,-fopt,'marker','o','markerSize',5)
title('bounded inputs')
xlabel('x value')
ylabel('y value')
zlabel('f value')
subplot(2,1,2)
plot(-fopt)
xlabel('iterations')
ylabel('function minimum')
%Case 3 study with bounded conditions
tic
options=optimoptions('ga');
options=optimoptions(options,'populationSize',300);
for i=1:no_cases
[inputs,fopt(i)]=ga(@stalagmite,2,[],[],[],[],[0;0],[0.6;0.6],[],[],options);
xopt(i)=inputs(1);
yopt(i)=inputs(2);
end
study_time3=toc;
figure(3)
subplot(2,1,1)
surfc(x,y,-f)
shading interp
hold on
plot3(xopt,yopt,-fopt,'marker','o','markerSize',5,'markerfacecolor','r')
title('Bounded inputs with population size')
xlabel('x value')
ylabel('y value')
zlabel('f value')
subplot(2,1,2)
plot(-fopt)
xlabel('iterations')
ylabel('function minimum')
fprintf('Study time for Study one is:%f \n',study_time1);
fprintf('Study time for Study two is:%f \n' ,study_time2);
fprintf('Study time for Study three is:%f \n',study_time3);
function:
function [f]=stalagmite(input_vector)
x=input_vector(1);
y=input_vector(2);
f_1x=(sin(5.1*pi*x+0.5))^6;
f_1y=(sin(5.1*pi*y+0.5))^6;
f_2x=exp(-4*log(2)*((x-0.0667)^2)/0.64);
f_2y=exp(-4*log(2)*((y-0.0667)^2)/0.64);
f=-(f_1x*f_1y*f_2x*f_2y);
endThese questions may occur when the whole program gets completed,
1.Why to use the loop for 50 times for ga?
We use number of cases to run ga,to get the maximum number of optimum values.In order to get the maximum number of optimum values ,we help matlab to search the whole searchspace.In other words, we're defining the number of generations that we are opting out for case study.
2.Why to use/increase the population size?
As we use the population size command in matlab,matlab will search the random number of population or points in searchspace per generation.With the increase in population size,we get more accurate results or the single global maxima.
CONCLUSION:
We get the Global maxima when we use bounded conditions with the defined population size.
Under Unbounded conditions,
for bounded conditions,
For bounded with population size condition,
Time taken to study the search space is ,