Particle swarm optimisation in matlab code

Particle swarm optimization is one of the most used algorithms for optimisation problems. This algorithm is inspired from nature, where flocks of birds move in the air, altering their positions frequently to reach the destination. 

Generally, in a flock of birds, there will be a lead, whom other birds follow. That lead will guide other birds in order to reach the target. This lead generally loses more energy when compared to other birds because of several factors of retardation forces that oppose the lead more than it’s members in the flock’s movement. Hence, it is not advisable for only one bird to guide the flock to the destination. Hence each bird in the flock leads when it’s turn come, so that energy loss is optimised.

Code:
clc;
clear;
close all;
nVar = 2;        % Number of Unknown (Decision) Variables
VarSize = [1 nVar];         % Matrix Size of Decision Variables
VarMin = -10;	% Lower Bound of Decision Variables
VarMax = 10;    % Upper Bound of Decision Variables
%% Parameters of PSO
MaxIt = 50;   % Maximum Number of Iterations
nPop = 4;     % Population Size (Swarm Size)
w = 1;           % Intertia Coefficient
wdamp = 0.99;   % Damping Ratio of Inertia Coefficient
c1 = 2;         % Personal Acceleration Coefficient
c2 = 2;         % Social Acceleration Coefficient
% The Flag for Showing Iteration Information
MaxVelocity = 0.2*(VarMax-VarMin);
MinVelocity = -MaxVelocity;
labels = {'P1','P2','P3','P4'};
%% Initialization
% The Particle Template
empty_particle.Position = [];
empty_particle.Velocity = [];
empty_particle.Cost = [];
empty_particle.Best.Position = [];
empty_particle.Best.Cost = [];
empty_particle.Energy = [];
empty_particle.Ch = [];
empty_particle.Alive = [];
% Create Population Array
particle = repmat(empty_particle, nPop, 1);
% Initialize Global Best
GlobalBest.Cost = inf;
c_i = 0;
% Initialize Population Members
for i=1:nPop
% Generate Random Solution
particle(i).Position = unifrnd(VarMin, VarMax, VarSize);
particle(i).Alive = 1;
% Initialize Velocity
particle(i).Velocity = zeros(VarSize);
particle(i).Energy = 20;
particle(i).Ch = 0;
end
for i = 1:nPop
% Evaluation
particle(i).Cost = Cost_function(i,particle,nPop);
% Update the Personal Best
particle(i).Best.Position = particle(i).Position;
particle(i).Best.Cost = particle(i).Cost;
% Update Global Best
if particle(i).Best.Cost < GlobalBest.Cost
GlobalBest = particle(i).Best;
c_i = i;
end
end
%Update cluster head
particle(c_i).ch = 1; 
% Array to Hold Best Cost Value on Each Iteration
BestCosts = zeros(MaxIt, 1);
C_Energy = zeros(MaxIt, 1);
y = [];
x = [];
for it=1:MaxIt
c_new = c_i;
for i=1:nPop
%Energy dissipated by non cluster head nodes to send the sensor values and their positions 
if particle(i).Ch == 0 && particle(i).Alive ==1
distance = (particle(i).Position(1)-particle(c_i).Position(1)).*(particle(i).Position(1)-particle(c_i).Position(1));
distance = distance + (particle(i).Position(2)-particle(c_i).Position(2)).*(particle(i).Position(2)-particle(c_i).Position(2));
particle(i).Energy = particle(i).Energy - 0.01.*distance;
if particle(i).Energy <=0
particle(i).Alive = 0;
end
end
%Energy dissipated by cluster head to base station
if particle(i).Ch == 1 && particle(i).Alive ==1
distance = (particle(i).Position(1)).*(particle(i).Position(1));
distance = distance + (particle(i).Position(2)).*(particle(i).Position(2));
particle(i).Energy = particle(i).Energy - 0.01.*distance;
if particle(i).Energy <=0
particle(i).Alive = 0;
end
end
end
for i=1:nPop
if particle(i).Alive ==1
% Update Velocity
particle(i).Velocity = w*particle(i).Velocity ...
+ c1*rand(VarSize).*(particle(i).Best.Position - particle(i).Position) ...
+ c2*rand(VarSize).*(GlobalBest.Position - particle(i).Position);
% Apply Velocity Limits
particle(i).Velocity = max(particle(i).Velocity, MinVelocity);
particle(i).Velocity = min(particle(i).Velocity, MaxVelocity);
% Update Position
particle(i).Position = particle(i).Position + particle(i).Velocity;
% Apply Lower and Upper Bound Limits
%particle(i).Position = max(particle(i).Position, VarMin);
%particle(i).Position = min(particle(i).Position, VarMax);
% Evaluation
particle(i).Cost = Cost_function(i,particle,nPop);
% Update Personal Best
if particle(i).Cost < particle(i).Best.Cost
particle(i).Best.Position = particle(i).Position;
particle(i).Best.Cost = particle(i).Cost;
% Update Global Best
if particle(i).Best.Cost < GlobalBest.Cost
GlobalBest = particle(i).Best;
c_new = i;
end            
end
y(i) = particle(i).Position(1);
x(i) = particle(i).Position(2);
end
end
%Energy dissipated by cluster head to all non cluster head nodes
for j = 1:nPop
if particle(j).Alive ==1
distance = (particle(c_i).Position(1)-particle(j).Position(1)).*(particle(c_i).Position(1)-particle(j).Position(1));
distance = distance + (particle(c_i).Position(2)-particle(j).Position(2)).*(particle(c_i).Position(2)-particle(j).Position(2));
particle(c_i).Energy = particle(c_i).Energy - 0.01.*distance;
end
if particle(c_i).Energy <=0
particle(c_i).Alive = 0;
end
end
cls = zeros(1,nPop);
cls(c_i)=1;
c=zeros(4,3);
c(cls==1,:)=[1 0 0]; % 1 is red
% ...
scatter(x,y,[],c);
hold on
scatter(20,20,'x')
hold off
text(x,y,labels,'VerticalAlignment','top','HorizontalAlignment','left')
xlim([-50,50]);
ylim([-50,50]);
pause(1.5);
if c_i~=c_new
%Energy dissipated by new cluster head
for j = 1:nPop
if particle(j).Alive ==1
distance = (particle(c_new).Position(1)-particle(j).Position(1)).*(particle(c_new).Position(1)-particle(j).Position(1));
distance = distance + (particle(c_new).Position(2)-particle(j).Position(2)).*(particle(c_new).Position(2)-particle(j).Position(2));
particle(c_new).Energy = particle(c_new).Energy - 0.01.*distance;
end
end
if particle(c_new).Energy <=0
particle(c_new).Alive = 0;
end
particle(c_i).Ch=0;
particle(c_new).Ch=1;
c_i = c_new;
end
% Store the Best Cost Value
BestCosts(it) = GlobalBest.Cost;
C_Energy(it) =0;
for j=1:nPop
C_Energy(it) = C_Energy(it) + particle(i).Energy;
end
% Display Iteration Informatio
disp(['Iteration ' num2str(it) ': Best Cost = ' num2str(BestCosts(it)) ': Particle =' num2str(c_i) ': Energy =' num2str(particle(c_i).Energy)]);
% Damping Inertia Coefficient
w = w * wdamp;
end
plot([1:MaxIt],C_Energy);
xlabel('No. of Iterations');
ylabel('Total Energy of Cluster');
Cost_function:
function z = Cost_function(i,particle,npop)
alpha = 0.4;
beta = 0.4;
gamma = 0.19;
avg_distance = 0;
avg_energy = 0;
for j = 1:npop
x =  (particle(j).Position(1) - particle(i).Position(1)).*(particle(j).Position(1) - particle(i).Position(1));
y = (particle(j).Position(2) - particle(i).Position(2)).*(particle(j).Position(2) - particle(i).Position(2));
avg_distance = avg_distance + sqrt(x + y);
end
avg_distance = avg_distance/npop;
for j = 1:npop
avg_energy = avg_energy + particle(j).Energy;
end
avg_energy = avg_energy/particle(i).Energy;
x =  (particle(i).Position(1) - 20).*(particle(i).Position(1) - 20);
y = (particle(i).Position(2) - 20).*(particle(i).Position(2) - 20);
avg_distance1 = sqrt(x + y);
z = alpha.*avg_distance + beta.*avg_energy+gamma.*avg_distance1;
end

In this code, 4 particles are defined initially with a fixed energy. A destination is fixed, to which the flock moves, and all 4 birds try to move to the destination from the initial place, by changing lead. The lead is adjusted on various factors, like energy each bird has, distance between destination and bird, etc. They do this for every iteration. The birds stop if the energy is completed or they reach the destination.

Cost function :

Cost function is the function which depends upon the factors which opposes the flocks movement. It depends upon 3 factors:

1. The average distance between lead and other members,
2. The average distance between destination, and
3. The average energy which is consumed for a particular bird to send information to other birds.

The cost function taken is : 

Cost function =  alpha.*avg_distance_between_lead_and_other_birds + beta.*avg_energy + gamma.*avg_distance_between_lead_and_destination.

Where alpha, beta, and gamma are constants which varies upon the situation considered.