Optimization using evolutionary theory

Aim: Optimization of stalagmite function using Genetic Algorithm

Objective:

1) Maximizing Stalagmite function using GA for x,y ranging from 0 to 0.6

2) Checking for repeatability

3) Plotting the function in 3D along with optimized values

Genetic Algorithm:

GA is a search process inspired by Charles Darwin's Natural Evolutionary Theory. The process is based on natural selection where the fitest individuals are allowed to produce offsprings, the weak and old are removed from population. The fitness in this search algorithm is determined by the objective function (function that has to be optimized), the search space becomes the population and each value in the search space becomes an individual or a chromosome.[1]

The process begins with an initial population which is a subset of the population (for simplicity this program uses an equally spaced array generated within the specified limits as initial population). Using the fitness function, individuals are selected based on pc (cross over probability) for reproduction. These are called parents. The parents then undergo crossover to generate offsprings which will be added to the population. Few individuals are discarded and few are mutated. Crossover, discarding and mutation lead to a new generation of population and the cycle starts again. The cycle is terminated upon achieving certain condition. Each step of the cycle is discussed in detail below.

Selection:[2]

There are quite a few techniques that are being extensively used in parent selection, here we use a Ranking scheme for parent selection, where the rank of an individual is based on the fitness and higher ranked individuals are given more priority during selection.

In the following function the population is sorted based on their fitness. Now their row number becomes the rank. The population is then divided into three sets having equal number of individuals, where set1 holds the highest ranked individuals, set3 holds the least and set2 has the in-betweens. Two-thirds of the first set, a third of the second and a sixth of the third set get selected for crossover.

• To account for a possibility of no selection happening due to small population a null vector was first assigned to the output array to avoid errors.
• The lengths of the arrays were floored as the probabilities being assigned may result in real numbers that aren't integers that can't act as limiters.
• Conditional statements 26,43 were used to avoid errors that may occur in case no assignments were done in preceding lines.
• Separate variable j was used when calling the output variable to avoid overwriting any previous assignments.

function selection = Selection(population)
%ps = 0.5,0.33,0.16
% 0.5 * 1st 1/3 population will become parents
% 0.33 2nd 1/3 population will become parents
% 0.16 3rd 1/3 population will become  parents

%initializing array
selection = [];

%first 1/3
for i = 1:floor(length(population(:,3))/3)
    a(i,:) = population(i,:);
end
a = rearrange(a);
for i = 1:floor(2*length(a(:,1))/3)
    selection(i,:) = a(i,:);
end

%second 1/3
j = 1;
for i = floor(length(population(:,3))/3)+1:floor(2*length(population(:,3))/3)
    a(j,:) = population(i,:);
    j = j+1;
end
a = rearrange(a);
if size(selection) ~= [0,0]
    j = length(selection(:,1))+1;
else
    j = 1;
end
for i = 1:floor(length(a(:,1))/3)
    selection(j,:) = a(i,:);
    j = j+1;
end

%third 1/3
j = 1;
for i = floor(2*length(population(:,3))/3)+1:length(population(:,3))
    a(j,:) = population(i,:);
    j = j+1;
end
a = rearrange(a);
if size(selection) ~= [0,0]
    j = length(selection(:,1))+1;
else
    j = 1;
end
for i = 1:floor(length(a(:,1))/6)
    selection(j,:) = a(i,:);
    j = j+1;
end
end​

Crossover:[3]

Crossover function is the techniques used to merge parent chromosomes and generate off springs. The number of off springs generated, number of parents involved in crossover to generate one off spring and the technique used for crossover vary based on the application. Here we are using a Weighted Arithmetic Recombination where the fitter parents is given a higher weight age (In this case 0.6).

`Child = α.x + (1-α).y`

A randomly arranged array of selected parents is used to generate off springs. One off spring is generated per pair of parents 1,2 and 2,3 and so on till n-1,n and n,1. Entries in first and second columns of the output array are the chromosomes, entry in the third column is the corresponding fitness function and entry four is age which is reset to 1 when off springs are introduced to the population.

function child = crossover(parent,optimizing_function)

%initializing child array
child = rand(length(parent(:,1)),4);

%Weighted algorithmic crossover for parents 1-2, 2-3, 3-4 ..... (n-1)-n
%resulting in one offspring per pair
for i = 1:length(parent(:,1))-1
    if parent(i,3) < parent(i+1,3)%comparing fitness
        x = parent(i,1:2);
        y = parent(i+1,1:2);
    else
        x = parent(i+1,1:2);
        y = parent(i,1:2);
    end
    child(i,1) = 0.6*(x(1))+0.3*(y(1));%twice the weightage to fitter parent
    child(i,2) = 0.6*(x(2))+0.3*(y(2));%twice the weightage to fitter parent
    child(i,3) = optimizing_function(child(i,1:2));
    child(i,4) = 1;
end

%Weighted algorithmic crossover for parents n,1
%resulting in one offspring
if parent(length(parent(:,1)),3) < parent(1,3)%comparing fitness
    x = parent(length(parent(:,1)),1:2);
    y = parent(1,1:2);
else
    x = parent(1,1:2);
    y = parent(length(parent(:,1)),1:2);
end
child(length(parent(:,1)),1) = 0.6*(x(1))+0.3*(y(1));%twice the weightage to fitter parent
child(length(parent(:,1)),2) = 0.6*(x(2))+0.3*(y(2));%twice the weightage to fitter parent
child(length(parent(:,1)),3) = optimizing_function(child(length(parent(:,1)),1:2));
child(length(parent(:,1)),4) = 1;
end

New Generation:[4]

New Generation is created by replacing few members of the existing population with the off springs. This can be done by replacing either old individuals (ones that have been a part of the population for more number of generations) or weak individuals (low fitness).

This project uses a combination where two tenth of the off springs replace oldest individuals and the remaining replace weak individuals. This is achieved by first sorting the population based on fitness and replacing the individual with 8/10 of the off springs starting from the bottom. Further the revised population is sorted based on age and remaining off springs are substituted into the array starting from the bottom.

function newgeneration = newgen(population,children)

%sorting based on fitness
sorted = foldsort(population);

%replacing week fits with offsprings
for i = 1:floor(length(children)*0.8)
    sorted(length(sorted)-i+1,:) = children(i,:);
end

%sorting based on age
Sortforage = sortforage(sorted);

%replacing older population with offsprings
j = 1;
for i = floor(length(children)*0.8)+1:length(children)
    Sortforage(length(Sortforage)-j+1,:) = children(i,:);
    j = j+1;
end

%output
newgeneration = Sortforage;
end

Mutation:[5]

Mutation is a crucial step to evolution which checks on pre mature convergence. It stops convergence to a local optimization by adding random individuals outside the population.

Following function replaces a random individual with random chromosomes(variables) within the search space and assigns corresponding fitness function while also resetting the age. The function does not deal with the number of mutations and hence the genetic algorithm has to call the function multiple times based on the population and mutation probability(pm) [for this case of total population of 200 and a pm of 0.2, the function is called 40 times].

function mutation_set = mutation(selected_parents,optim_func)

a = length(selected_parents(:,1));
a = randi([1,a]);%random integer between 1 and length of the parent array

%assigning values of the random element from population to random element in parent
I = [0.6*rand,0.6*rand];
f = optim_func(I);
selected_parents(a,1:2) = I;
selected_parents(a,3) = f;
selected_parents(a,4) = 1;

%output
mutation_set = selected_parents;
end

Termination:[6]

Once a new generation is obtained after mutation, the population then becomes initial population for the next cycle. A genetic algorithm needs a termination condition to check when the cycle has to stop. This can be achieved by comparing the population in successive generations and cycle can be terminated when there is no significant difference in the population. Another way to do this is by terminating after a predefined number of cycles or when the population has achieved a certain level of fitness. As in this case we do not know the exact fitness threshold nor acceptable difference in the population we will be terminating the function after 100 iterations/generations.

The GA function:

Function takes in the fitness function(optim_function) and the limiting values(I,J) of the search space to minimize the fitness function and return the minimized value(population(1,3)) and corresponding state points where the optimization has occurred(population(1,1:2)).

The heuristic starts with an initial population which is a random selection within the search space. To achieve the same a 400 data point array was generated using a mesh grid of equally spaced data points and half of them were chosen in random. Although is not a perfect definition of random, allows for data points to be spread across the search space rather than a truly random set that could concentrate the initial population locally.

A population array is defined where the first two columns are the variables of the function to be optimized, the third represents the fitness corresponding to the variables and the fourth is the age of the individual (number of generations the individual has been in the population).

Age of the entire population is then incremented before starting the obtaining a new generation. The population is then subjected to selection, crossover, substitution and mutation as already discussed to arrive at a new generation. The process is run for 100 generations and the fittest in the final population is returned to main function

%ps = 0.66,0.33,0.16
%pc = 0.33*generation population [0.38 = 0.66/3 + 0.33/3 + 0.16/3]
%pm = 0.2*generation population
%pold = 0.2*children
%pweek = 0.8*children

%procedure:
%initial population 300
%fitness
%sorting based on fitness
%selection
%crossover
%mutation
%new generation
%go back to line9 

% optim_function = @stalagmite;
% I = [0,0.6];
% J = I;

function output = geneticalgorithm(optim_function,I,J)

%input
ll1 = I(1);
ul1 = I(2);
ll2 = J(1);
ul2 = J(2);

%initial data set
x = linspace(ll1,ul1,20);
y = linspace(ll2,ul2,20);
[X,Y] = meshgrid(x,y);

k = 1;
for i = 1:length(X(:,1))
    for j = 1:length(X(1,:))
        %column1-variable,column2-variable,column3-fitness,column4-age:
        original_population(k,:) = [X(i,j),Y(i,j),0,1];
        k = k+1;
    end
end

%rearranging
original_population = rearrange(original_population);

%population initialisation
for i = 300
    population(i,:) = original_population(i,:);
end

%fitness function
for i = 1:length(population(:,1))
    population(i,3) = optim_function(population(i,:));
end

%sort
sorted = foldsort(population);
population = sorted;

for count = 1:100
    %Aging before going to new generation
    population(:,4) = population(:,4)+1;


    %selection by rank
    selection = Selection(population);
    parents = rearrange(selection);

    %crossover by weighted arithmatic recombination
    children = crossover(parents, optim_function);

    %new generation
    Newgen = newgen(population,children);
    population = Newgen;

    for count2 = 1:60
        %mutation by random resetting
        population = mutation(population, optim_function);
    end

    population = rearrange(population);

    %sort
    sorted = foldsort(population);
    population = sorted;

end

%output
output = population(1,1:3);

end

Sorting and Rearranging:

1)foldsort function:

This function is used to sort the array based on fitness with the fittest on top to the least in bottom

function A = foldsort(A)
i = 1;
j = 1;
for i = 1:length(A(:,1))
    for j = 1:length(A(:,1))
        if A(i,3) < A(j,3);
            %swapping
            b = A(i,:);
            A(i,:) = A(j,:);
            A(j,:) = b;
        end
    end
end

2)sortforage function:

This function is used to sort the array based on age with the youngest on top to the oldest in bottom

function A = sortforage(A)
i = 1;
j = 1;
for i = 1:length(A(:,1))
    for j = 1:length(A(:,1))
        if A(i,4) < A(j,4);
            %swapping
            b = A(i,:);
            A(i,:) = A(j,:);
            A(j,:) = b;
        end
    end
end     

3)rearrange function:

Arranges the array using a random permutation, hence shuffling the data randomly.

function Data = rearrange(Data1)
if size(Data1) == [0,0]
    Data = Data1;
elseif length(Data1(:,1)) == 1
    Data = Data1;
else
a = randperm(length(Data1(:,1)));
for i = 1:length(Data1(:,1))
    Data(i,:) = Data1(a(i),:);
end
end
end

Stalagmite Function:

`f1,x = sin^6(0.5+5.1*pi*x)`

`f1,y = sin^6(0.5+5.1*pi*y)`

`f2,x = e^(-4*log(2)*((x-0.0667)^2)/0.64)`

`f2,y = e^(-4*log(2)*((y-0.0667)^2)/0.64)`

`F(x,y) = f1,x * f1,y * f2,x * f2,y`

Input I = [x,y]

function f = stalagmite_func(I)
A = (sin(0.5+(5.1*pi*I))).^6;
B = exp(-4*log(2)*((I-0.0667).^2)/0.64);
f = -A(1)*B(1)*A(2)*B(2);
end

Main:

Arrays for x and y are generated with 150 equally spaced values between 0 and 0.6. A mesh grid is generated and the corresponding stalagmite values are assigned to f and plotted in a surface plot of x vs y vs f.

Genetic algorithm is called by sending the stalagmite function and the search space limits. Since the algorithm minimizes the function, stalagmite function has been defined to return the -ve values, therefore the final values have a negative multiplier to neutralize the change before plotting. The GA function is repeated to check for repeatability.

Optimum values obtained from each iteration are plotted in the surface plot along with a plot for fitness values of each set to visualize the change.

clear all
close all
clc

%Input
x = linspace(0,0.6,150);
y = linspace(0,0.6,150);
[X,Y] = meshgrid(x,y);
optimization_func = @stalagmite_func;

for i = 1:length(X)
    for j = 1:length(Y)
      in(1) = X(i,j);
      in(2) = Y(i,j);
        f(i,j) = -optimization_func(in);
    end
end


for i = 1:50
    optim_result(i,:) = geneticalgorithm(optimization_func,[0,0.6],[0,0.6]);
end

subplot(2,1,1)
surfc(X,Y,f)
xlabel('X')
ylabel('Y')
zlabel('Stalagmite')
hold on
shading interp
plot3(optim_result(:,1),optim_result(:,2),-optim_result(:,3),'*','color',[0 0 0])
subplot(2,1,2)
plot(-optim_result(:,3))
ylabel('Stalagmite optimization results')
xlabel('Iteration')
axis([1 length(optim_result(:,3)) 0.5 1.1])

Result & Discussion:

The program was studied for various Initial population used in the GA function along with changing mutation probability. Results as follows:

figure 1: Initial population: 100 | No of Iterations: 100 | No of Mutations: 40

figure 2: Initial population: 200 | No of Iterations: 100 | No of Mutations: 40

figure 3: Initial population: 300 | No of Iterations: 100 | No of Mutations: 40

figure 4: Initial population: 300 | No of Iterations: 100 | No of Mutations: 30

figure 5: Initial population: 300 | No of Iterations: 100 | No of Mutations: 50

figure 6: Initial population: 300 | No of Iterations: 100 | No of Mutations: 60

As the plots show increasing the population size and no of mutations improve the quality of fit obtained. We should note that this conclusion is only for the studied range and cannot be used to hypothesis that a larger population or pm always indicate a better fit. One must also consider that the computation time is largely increased for a small change in these parameters.

Scope:

1) Characteristics of the GA function can be studied based on the standard deviation obtained from results in each iteration of the function with respect to increasing population, termination condition and mutation permutation to understand the function better.

2) More over the effect of fixed quantities used in this project like a) the priorities given to individuals based on ranking to get selected for crossover, b) the weight age given in the crossover algorithm and c)the ratio of off springs taking over the old and the weak can be studied further to make the function better.

3) Also the effects of above parameter have to be studied on various functions to study if the GA can produce similar results for any function or the parameter optimization is different for different functions.

References:

1)https://towardsdatascience.com/introduction-to-genetic-algorithms-including-example-code-e396e98d8bf3

2)https://www.tutorialspoint.com/genetic_algorithms/genetic_algorithms_parent_selection.htm

3)https://www.tutorialspoint.com/genetic_algorithms/genetic_algorithms_crossover.htm

4)https://www.tutorialspoint.com/genetic_algorithms/genetic_algorithms_survivor_selection.htm

5)https://www.tutorialspoint.com/genetic_algorithms/genetic_algorithms_mutation.htm

6)https://www.tutorialspoint.com/genetic_algorithms/genetic_algorithms_termination_condition.htm