Ant System Implementation to solve the Traveling Salesman Problem (berlin52 dataset).

README.md

Ant System for Solving the Traveling Salesman Problem

This implementation of the Ant System (a variation of Ant Colony Optimization) [1] aims to solve the Traveling Salesman Problem. The problem is to find the shortest tour distance given a list of cities represented by its x and y coordinates where each city is visited only once. This was tested using the berlin52 dataset found [here] (http://elib.zib.de/pub/mp-testdata/tsp/tsplib/tsp/berlin52.tsp).

Files Included:

1. ant_colony.m
   This is the main file for the algorithm implementation. This calls other functions (written in separate files) in order to achieve the optimization objective.
2. ant_tour_construction.m
   For each ant in the colony, we make them traverse the whole map (tour). They first start at a randomly generated starting area (start_places), and jump to each node. The jumps are controlled by the attractiveness of the node raised to beta, and the pheromone effect raised to alpha.
3. compute_cost.m
   For each ant, we look at the path it took (also known as its tour), and then we base the distance from the distance matrix and add it to the link.
4. update_tau.m
   Here, we're using Ant-Quantity AS variation. In this case, only local information, d_{ij}, is used to update pheromone concentrations. Lower cost links are made more desirable. If d_{ij} represents the distance between links, then ant-quantity AS prefers selection of the shortest links.

Installation:

First, make sure that you have MATLAB installed. Compatibility to open-source software, such as Octave, has not yet been tested. If you're set, then just clone this repository:

$ git clone https://github.com/ljvmiranda921/tsp-ant-system.git`  

Usage:

Import the berlin52.csv file from your local path:

cd('your\local\path\to\file')

Change your Ant System parameters:

Parameter	Description
iterations	Number of iterations that the algorithm will run.
colony_size	Number of ants that will explore the search space.
evap_coeff	Rate of pheromone degradation.
alpha	Exploitation parameter, sets how the ants are attracted to pheromone concentration.
beta	Exploration parameter, sets how the ants are more attracted to try out shorter paths.
tau	Initial pheromone concentration.
el	Parameter for eliminating common costs.
Q	Positive constant for the Ant-Quantity System

Run ant_colony.m

Experiment on different values of alpha and beta in order to balance the exploitation-exploration tradeoff.

Implementation Notes:

• Parameter Tuning
  This implementation, albeit simple, was able to generate a tour distance of 7548.9927, a little bit shy of the most optimal solution for the berlin52 dataset (7542) [2]. In this case, I used the following parameters:

iterations = 100;          
colony_size = 30;          
evap_coeff = 0.15;         
alpha = 9;                            
beta = 12;                  
tau = 0.0001 * ones(52);                    
el=0.97;                    
Q = 0.2;                    

License

This project is licensed under the MIT License - see the LICENSE.txt file for details

References

[1] A.P. Engelbrecht, "Ant System" in Computational Intelligence: An Introduction, 2nd Edition, John Wiley & Sons, Ltd., pp. 368-372, 2007.
[2] "Optimal Solutions for Symmetric TSPs", Available: http://comopt.ifi.uni-heidelberg.de/software/TSPLIB95/STSP.html [Accessed: November 28, 2016].
[3] I. Sustrino, "Ant Colony Optimization to solve traveling salesman problem (TSP)," Available: https://www.mathworks.com/matlabcentral/fileexchange/51113-ant-colony-optimization--aco--to-solve-traveling-salesman-problem--tsp-, 2015. [Accessed: November 27, 2016]. Note: The partitioning of different ant-colony functions into stand-alone MATLAB m-files was adopted from Sustrino's work for this project.

ant_colony.m

%**************************************************************************
%           CI: Travelling Salesman Problem using ACO [AS Variation]
% -------------------------------------------------------------------------
% Name: MIRANDA, Lester James Validad
% ID No: 44161652-3
% -------------------------------------------------------------------------                             
clearvars;


%% =================== Part 0: Parameter setting =====================================
iterations = 500;           % Sets the number of iterations where the algorithm will run.
colony_size = 60;           % Sets the number of ants in the colony.
evap_coeff = 0.1;          % Sets the evaporation coefficient that will be used in updating the pheromone (must be between 0 and 1).

alpha = 15;                  % Sets the bias of using pheromone deposits in the decision process. 
                            % (Heuristic parameter describing how greedy the algorithm is in finding its path across the graph.)
                            % Exploitationn parameter - set this to 0 and
                            % the algo becomes greedy (random spikes). Set this high and it
                            % will just dig deep to the local optima. The
                            % problem is that it may find not good answers
                            
beta = 20;  %12                % Sets the bias for deciding if a shorter path is better.
                            % (Heuristic parameter describing how fast the ants are going to converge to a steady solution.) 
                            % Exploration parameter - set this to 0 and
                            % the attractiveness of the moves becomes
                            % neglected. It will just keep digging deep
                            % into the local optima. 

tau = 0.0001 * ones(52);    % We create a 52 x 52 tau matrix full of 1. We initialize it to a very small value (not 0 because it will 
                            % cause the system to divide by 0 later, which is
                            % bad. 
                         
el=0.97;                     % Set the parameter for eliminating common costs. 

Q = 0.2;                     % Positive constant used for updating tau. To control the level of exploration
                             % undertaken by the ants
                             % History: 1/10000 (11/25/2016) Answer: 8000 range

%% =================== Part 0.1: Load berlin52 dataset from csv file ===================
cd('C:\Users\Lj Miranda\Documents\Waseda\2016-01 Fall\110400F Computational Intelligence\reports\report-ci-final\tsp\ant-colony')
berlin52 = csvread('berlin52.csv');
x = berlin52(:,1);
y = berlin52(:,2);


%% =================== Part 0.2: Plot berlin52 dataset ==================================
figure
hold on;

subplot(2,2,1);
mycolor = [0.93,0.273,0.191];
plot(x,y,'o','MarkerFaceColor',mycolor,'MarkerEdgeColor','k','MarkerSize',5)
title('Plot of Berlin52 Dataset')
xlabel('x  (km)')
ylabel('y  (km)')
axis square

%% =================== Part 1: Initialize parameters ====================================

% Generate distance for each node

for i=1:52
    for j=1:52
        distance_matrix(i,j)=sqrt((x(i)-x(j))^2+(y(i)-y(j))^2);
    end
end

% Generate a priori attractiveness

for i=1:52
    for j=1:52
        if distance_matrix(i,j)==0
            a_priori(i,j)=0;
        else
            a_priori(i,j)=1/distance_matrix(i,j);
        end
    end
end

%% =================== Part 2: Ant System Algorithm ====================================

for i=1:iterations % We run the ant system algorithm with {iterations} number of iterations.
    % ---------- Part 2.1 Randomly Generate Places for each ant -----------
    for j = 1:colony_size
        start_places(j,1) = fix(1+rand*(52-1));     % fix() rounds each element toward 0.
    end
    % ---------- Part 2.2 Probabilistic solution construction -------------
    [tour] = ant_tour_construction(start_places, colony_size, a_priori, tau, alpha, beta);
    tour = horzcat(tour,tour(:,1));                 % horzcat() concatenates arrays horizontally.
    
    % ---------- Part 2.3 Calculate the total distance traversed-----------
    [cost,reference] = compute_cost(colony_size, distance_matrix, tour, el);
    [tau] = update_tau(colony_size, tau, tour, reference, evap_coeff, Q);
    average_cost(i) = mean(cost);
    
    % ---------- Part 2.4 Find optimized route ----------------------------
    [min_cost(i), best_index] = min(cost);
    besttour(i,:) = tour(best_index,:);
    iteration(i) = i;
    mean_cost(i) = mean(cost);
    a(i)=min(min_cost);
    
    %% =================== Part 3: Plot results ============================================
    
    % Plot Average of tour distance vs Number of Iterations
    subplot(2,2,3);
    %plot(iteration,min_cost);
    %plot(iteration,mean_cost);
    plot(iteration,a);
    title('Tour vs Iterations');
    xlabel('Iterations');
    ylabel('Tour');
    axis square
    
    % Plot Ant Behaviour
    subplot(2,2,4);
    plot(iteration,mean_cost);
    title('Ant Behaviour');
    xlabel('Iterations');
    ylabel('Tour');
    axis square
    
    
   % Plot the best route
    [k,l]=min(min_cost);
    for i=1:52+1
    X(i)=x(besttour(l,i));
    Y(i)=y(besttour(l,i));
    end
    subplot(2,2,2);
    plot(X,Y,'--ko','MarkerEdgeColor','k','MarkerFaceColor',[1 1 0],'MarkerSize',5)
    grid on;
    xlabel('x (km)');ylabel('y (km)');
    title(['Cost = ',num2str(k)]);
    %dim = [.9 .2 .6 .6];
    %str = {'Parameters',strcat('\alpha: ', num2str(alpha)), strcat('\beta: ', num2str(beta)), strcat('m: ', num2str(colony_size)), strcat('Q: ', num2str(Q)), strcat('\rho: ',num2str(evap_coeff))};
    %annotation('textbox',dim,'String',str,'FitBoxToText','on')
    axis square 
    
    pause(.1) %// pause 0.1 seconds to slow things down
    drawnow  
end

ant_tour_construction.m

function [new_places] = ant_tour_construction(start_places, colony_size, a_priori, tau, alpha, beta)
% SUMMARY: This function completes the ant tour matrix for one iteration.
%   For each ant in the colony, we make them traverse the whole map (tour). 
%   They first start at a randomly generated starting area (start_places), 
%   and jump to each node. The jumps are controlled by the attractiveness 
%   of the node raised to beta, and the pheromone effect raised to alpha. 

for i = 1:colony_size   % For each ant...
    a_priori_temp = a_priori;   % ... we track its sight// sight is what it saw, or what cities it saw. 
    for j = 1:51                % We then have it move around the 51 cities (not 52 because 1 is already starting point).
        c = start_places(i,j);      % We track its place and put it into c.
        a_priori_temp(:,c) = 0;
        temp = (tau(c,:).^alpha).*(a_priori_temp(c,:).^beta);
        s = (sum(temp));
        p =(1/s).* temp;
        r = rand; % rand just generates from 0 to 1. 
        s = 0;
        for k = 1:52
            s = s + p(k);
            if r <= s
                start_places(i,j+1)=k;
                break
            end
        end
    end
end
new_places = start_places;

compute_cost.m

function [cost, reference] = compute_cost(colony_size, distance_matrix, tour, el)
% SUMMARY: This function computes the cost of traversed by the ant, and the
% amount of tau that must be decayed.
%   For each ant, we look at the path it took <<tour>>, and then we base
%   the distance from the distance matrix and add it to the link.

for i = 1:colony_size % For each ant,...
    dist = 0;
    for j=1:52
        dist = dist + distance_matrix(tour(i,j),tour(i,j+1)); % ...we look at the tour table, then get the distance traversed with the help of distance_matrix, and store it in s.
    end
    reference(i)= dist; % ... we make a reference table, for ant 1, we have a cost of 1 (dist) and so on.
end
cost = reference; 
reference = reference - el * min(reference); % ... eliminates repeating distances

update_tau.m

function [tau]=update_tau(colony_size, tau, tour, reference, evap_coeff, Q)
%SUMMARY: The idea here is that we update the pheromones. We take into
%consideration the evaporatoin coefficient.
%   Here, we're using Ant-Quantity AS variation. In this case, only local information, d_{ij} , is used to update pheromone 
%   concentrations. Lower cost links are made more desirable. If d_{ij} represents the distance between links, then ant-quantity AS 
%   prefers selection of the shortest links.

for i= 1:colony_size % For each ant...
    for j=1:52 % ... we update certain values of tau depending on the tour table
        delta_tau = Q/reference(i); % This is an Ant-Quantity AS
        tau(tour(i,j),tour(i,j+1)) = (1-evap_coeff) * tau(tour(i,j),tour(i,j+1)) + delta_tau;
    end
end