Shortest path using Particle Swarm Optimisation

Question

I want to solve the Shortest Path problem using PSO in MATLAB.I encoded the path using priority encoding [1], and I'm using constriction and velocity clamping [2].

The problem I'm facing is that the code is extremely slow compared to Dijkstra. I, first, test using Dijkstra to get a benchmark time, and then run PSO to find the lowest cost it can achieve in that time. The result of PSO is always much higher.

If I check how fast each iteration finishes, I find that it takes few seconds for a path with 1000+ nodes on an Intel Core i3-2120 processor.

In the code below, you need to run data.m first to initialise the cost matrix, and then run Dijkstra to get a time benchmark. After that, modify the allowedTime variable in pso.m in seconds.

Parameters:

• data.m
  • dimensions: no. of nodes
• pso.m
  • allowedTime: time in seconds allowed for the swarm to run
  • swarm_size: no. of particles
  • startNode: no. representing where to start the path (within the dimensions range)
  • endNode: no. representing where to end the path (within the dimensions range)
• dijkstra.m
  • accepts (costMatrix, <start_node_id>, <end_node_id>)

I'm sorry about the messy code and not using functions, but I needed to make everything inline and see all values after the code is done, or when I break.

data.m

% <<<<<<<<<<<<<<<<<< data definition >>>>>>>>>>>>>>>> %

clear;
clc;

fprintf('Generating data ...\n\n');

dimensions = 5000;

costMatrix = randi(dimensions, dimensions);

fprintf('DONE!\n\n');

pso.m

%% initialization
clc;

fprintf('Initialising swarm ...\n\n');

% parameters
% >>>>>>>>>>>>>>>>>>> edit <<<<<<<<<<<<<<<<<<< %
allowedTime = 15 * 60;
swarm_size = 50;

% SP algorithm specific.
startNode = 1;
endNode = 2;

% velocity equation params.
correction_factor_p = 2.05;
correction_factor_g = 2.05;
contrictionFactor = 0.72984;
% ^^^^^^^^^^^^^^^^^^^ end ^^^^^^^^^^^^^^^^^^^ %

gbest = 1;
oldFitness = 1000000001;
iterations = 0;

% pre-allocate arrays.
swarmPos = zeros(swarm_size, dimensions);
swarmVel = zeros(swarm_size, dimensions);
swarmBestPos = zeros(swarm_size, dimensions);
swarmBestPath = cell(swarm_size, 1);
swarmBestFit = zeros(1, swarm_size);
result = zeros(1, swarm_size);
upperBound = zeros(1, dimensions);
lowerBound = zeros(1, dimensions);

% set bounds.
for i = 1 : dimensions
    upperBound(i) = 100;
    lowerBound(i) = -100;
end

% ---- initiate swarm -----
for i = 1 : swarm_size
    for j = 2 : dimensions
        swarmPos(i,j) = lowerBound(j) + rand * (upperBound(j) - lowerBound(j));
        swarmVel(i,j) = rand * (upperBound(j) - lowerBound(j)) / 2;
        swarmBestPos(i,j) = swarmPos(i,j);
        swarmBestPath{i}(j) = -1;
        swarmBestFit(i) = 1000000000;   % best fitness so far
    end
end

% set starting node to avoid on accidental access.
for i = 1 : swarm_size
    swarmPos(i,1) = -99999999;
    swarmVel(i,1) = -99999999;
    swarmBestPos(i,1) = -99999999;
    swarmBestPath{i}(1) = startNode;
end

% ^^^^^^^^^^^^^^^^ END: initialisation ^^^^^^^^^^^^^^^^ %

% >>>>>>>>>>>>>>>>>>> START: swarming <<<<<<<<<<<<<<<<<<< %
clc;
fprintf('Swarming ...\n\n');
tic;
%% iterations
while true

    % reset results to allow summing.
    parfor i = 1 : swarm_size
        result(i) = 0;
    end

    % <<<<<<<<<<<<<<<<< START: movement and fitness >>>>>>>>>>>>>>>>> %

    for i = 1 : swarm_size
        for j = 2 : dimensions
            swarmPos(i,j) = swarmPos(i,j) + swarmVel(i,j);      % update x position

            if (swarmPos(i,j) > upperBound(j))
                swarmPos(i,j) = swarmPos(i,j) - (swarmPos(i,j) - lowerBound(j)) / 2;
            elseif (swarmPos(i,j) < lowerBound(j))
                swarmPos(i,j) = swarmPos(i,j) + (lowerBound(j) - swarmPos(i,j)) / 2;
            end
        end

        %tic;

        % <<< inline fitness function >>> %
        tempPath = [];
        tempPath(1) = startNode;
        invalidBuild = false;
        tempPos = swarmPos(i,:);

        for j = 2 : dimensions
            for k = 2 : dimensions
                [discard, maxPos] = max(tempPos);
                cost = costMatrix(tempPath(j - 1), maxPos);
                tempPos(maxPos) = -9999999 - k;

                if (cost < 100000000)
                    tempPath(j) = maxPos;
                    result(i) = result(i) + cost;
                    break;
                elseif (k == dimensions)
                    invalidBuild = true;
                end
            end

            if (invalidBuild)
                result(i) = 1000000000;
                break;
            elseif (tempPath(j) == endNode)
                break;
            end
        end
        % ^^^ END: fitness function ^^^ %


        % if new position is better
        if result(i) < swarmBestFit(i)
            for d = 1 : dimensions
                swarmBestPos(i,d) = swarmPos(i,d);  % update best x,
            end

            swarmBestPath{i} = tempPath;
            swarmBestFit(i) = result(i);        % and best value
        end
    end

    % ^^^^^^^^^ END: movement and fitness ^^^^^^^^^ %

    % <<<<<<<<<<<<<<<<< update global best >>>>>>>>>>>>>>>>> %
    for i = 1 : swarm_size
        if swarmBestFit(i) < swarmBestFit(gbest)
            gbest = i;                  % global best i.

            took = toc;     % time taken to reach this best.
        end
    end

    % <<<<<<<<<<<<<<<<< update velocity >>>>>>>>>>>>>>>>> %
    for i = 1 : swarm_size
        for j = 2 : dimensions
            swarmVel(i,j) = contrictionFactor * (swarmVel(i,j) ...
                                + correction_factor_p * rand * (swarmBestPos(i,j) - swarmPos(i,j)) ...
                                + correction_factor_g * rand * (swarmBestPos(gbest,j) - swarmPos(i,j)));

            if (swarmVel(i,j) > (upperBound(j) - lowerBound(j)) / 2)
                swarmVel(i,j) = (upperBound(j) - lowerBound(j)) / 2;
            end
        end
    end

    % <<<<<<<<<<<<<<<<< print global bests if changed >>>>>>>>>>>>>>>>> %
    if ( oldFitness ~= swarmBestFit(gbest) )
        oldFitness = swarmBestFit(gbest);

        % update display
        clc
        fprintf('Best particle position:\n');
        sizeTemp = size(swarmBestPath{gbest}, 2);
        for i = 1 : sizeTemp
            if (swarmBestPath{gbest}(i) ~= 0)
                fprintf('%d\n', swarmBestPath{gbest}(i));
            end
        end
        fprintf('\nBest fitness: %d\n\n', swarmBestFit(gbest));
    end

    iterations = iterations + 1;

    % end on timeout
    elapsedTime = toc;
    if (elapsedTime > allowedTime)
        break;
    end
end

clc;

fprintf('>>>>>>>>>>>>>>> FINISHED <<<<<<<<<<<<<<<\n\n\n');

fprintf('Best path:\n');
sizeTemp = size(swarmBestPath{gbest}, 1);
for i = 1 : sizeTemp
    if (swarmBestPath{gbest}(i) ~= 0)
        fprintf('%d\n', swarmBestPath{gbest}(i));
    end
end
fprintf('\nBest cost: %d\n\n', swarmBestFit(gbest));
fprintf('\nTook: %d iterations, and %.2f seconds.\n\n', iterations, took);

dijkstra.m

function dijkstra(matriz_costo, s, d)
% This is an implementation of the dijkstra´s algorithm, wich finds the 
% minimal cost path between two nodes. It´s supoussed to solve the problem on 
% possitive weighted instances.

% the inputs of the algorithm are:
%farthestNode: the farthest node to reach for each node after performing
% the routing;
% n: the number of nodes in the network;
% s: source node index;
% d: destination node index;

%For information about this algorithm visit:
%http://en.wikipedia.org/wiki/Dijkstra%27s_algorithm

%This implementatios is inspired by the Xiaodong Wang's implememtation of
%the dijkstra's algorithm, available at
%http://www.mathworks.com/matlabcentral/fileexchange
%file ID 5550

%Author: Jorge Ignacio Barrera Alviar. April/2007


n=size(matriz_costo,1);
S(1:n) = 0;     %s, vector, set of visited vectors
dist(1:n) = inf;   % it stores the shortest distance between the source node and any other node;
prev(1:n) = n+1;    % Previous node, informs about the best previous node known to reach each  network node 

dist(s) = 0;

iterations = 0;

tic;
while sum(S)~=n
    candidate=[];
    for i=1:n
        if S(i)==0
            candidate=[candidate dist(i)];
        else
            candidate=[candidate inf];
        end
    end
    [u_index u]=min(candidate);
    S(u)=1;
    for i=1:n
        if(dist(u)+matriz_costo(u,i))<dist(i)
            dist(i)=dist(u)+matriz_costo(u,i);
            prev(i)=u;
        end
    end

    iterations = iterations + 1;
end


sp = [d];

while sp(1) ~= s
    if prev(sp(1))<=n
        sp=[prev(sp(1)) sp];
    else
        error;
    end
end;
spcost = dist(d);
took = toc;

fprintf('Best path:\n');
fprintf('%d\n', sp);
fprintf('\nBest cost: %d\n\n', spcost);
fprintf('\nTook: %d iterations, and %.2f seconds.\n\n', iterations, took);

(1) A Nondominated Sorting Genetic Algorithm for SP Routing Problem
(2) Constriction factors and Parameters

Answer 1

I'm a newbie to the world of Artificial Intelligence. I have tried my level best to help you.

PSO is a metaheuristic approach. Problems with incomplete or imperfect information or limited computation capacity can be solved using metaheuristics. Such problems cannot be solved using Dijkstra's as it needs full topology details. Hence use of the algorithm depends on the problem domain.

Since PSO is a stochastic process, initial result will not be the optimal. As the iterations are performed, the cost function value reduces. Where as Dijkstra's finds the shortest path by one iteration.

So time consumption for PSO will be more compared to Dijkstra. Use of these algorithms are problem specific.

Hope this will help you!