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.

Categories: Matlab

0 Comments

Leave a Reply