Code covered by the BSD License

# Ant System TSP Solver

by

### Johannes (view profile)

26 Oct 2011 (Updated )

A demo of an Ant System algorithm solving classical Traveling Salesman Problems.

ant_system_tsp(map_fct, stopeval, varargin)
```function [popt, fopt, stat] = ant_system_tsp(map_fct, stopeval, varargin)
% [popt, fopt, stat] = ant_system_tsp(map_fct, stopeval, varargin)
%
% Implemenetation of an Ant System algorithm for solving TSP problems
%
% Input:
% - map            - A city map, which is a matrix of 2D city coordinates
% - stopeval       - The number of function evaluations used
%                    by the optimization algorithm
%
% Options:
% - m              - The number of ants (default: 100)
% - rho            - The evaporation rate (default: 0.1)
% - alpha          - The pheromone influence parameter (default: 1)
% - beta           - The heuristic influence parameter (default: 3)
% - visualization  - Visualize the run of the algorithm (default: true)
%
% Output:
% - popt           - A vector containing the optimal path
% - fopt           - The length of the optimal path
% - stat           - A MATLAB struct containing algorithm
%                    specific output statistics
%
% Author: Johannes W. Kruisselbrink

% Parse the inputs
p = inputParser;
p.parse(map_fct, stopeval, varargin{:});

% Set the algorithm parameters
m = p.Results.m;
rho = p.Results.rho;
alpha = p.Results.alpha;
beta = p.Results.beta;
visualization = p.Results.visualization;

% Obtain the map from the function handle
map = feval(map_fct);

% The number of cities
l = length(map(:,1));

% Compute a reference path using the nearest neighbor method
[path_nn, C_nn] = nn_shortest_path_tsp(map);

% Compute distance matrix
dmat = -1 * ones(l,l);
for i = 1:l
dmat(i,:) = sqrt(sum((map(i * ones(l,1),:) - map).^2, 2));
end

% Initialize pheromones and heuristics
pheromones = ones(l,l) * m ./ C_nn;
heuristics = ones(l,l) ./ (1e-10 * eye(l) + dmat);
heuristics(1:l+1:l*l) = 0;

evalcount = 0;
fopt = Inf;
stat.histf = zeros(1, stopeval);

% Ant System Optimization Loop
while evalcount < stopeval

% Reset the ant paths
path = NaN(m,l);

% Pre-compute the startnodes of the ants
startnode = ceil(rand(1,m) * l);

% Pre-compute the product of the pheromones and heuristics (for efficiency)
pheromones_heuristics_product = pheromones.^alpha .* heuristics.^beta;

% Construct new paths for the ants
for k = 1:m

curnode = startnode(k);
path(k,1) = curnode;

% Maintain a private pheromones and heuristics matrix (update it w.r.t. the tabu-list)
private_pheromones_heuristics_product = pheromones_heuristics_product;
private_pheromones_heuristics_product(:,curnode) = 0;

for j = 2:l

% Compute the step probabilities
current_node_attraction_scores = private_pheromones_heuristics_product(curnode,:);
step_probabilities = current_node_attraction_scores ./ sum(current_node_attraction_scores);

% Draw the next node proportional to the step probabilities
cumsum_step_probabilities = cumsum(step_probabilities);
r = rand() * cumsum_step_probabilities(l);
curnode = 1;
while (cumsum_step_probabilities(curnode) < r)
curnode = curnode + 1;
end
path(k,j) = curnode;

% Update the tabu-list: set the probabilities for this node to zero
private_pheromones_heuristics_product(:,curnode) = 0;

end

end

% Compute fitness
for k = 1:m
C(k) = tsp_evaluate_tour(map, path(k,:));

evalcount = evalcount + 1;
stat.histf(evalcount) = fopt;

end

[fopt, optindex] = min(C);
popt = path(optindex, :);

% Update pheromones
new_pheromones = zeros(l,l);
for k = 1:m
for j = 1:l-1
new_pheromones(path(k,j), path(k,j+1)) = new_pheromones(path(k,j), path(k,mod(j,l) + 1)) + (1 / C(k));
end
end
pheromones = (1 - rho) * pheromones + new_pheromones;

% Plot statistics
if (visualization)
clf

subplot(2,2,1)
plot(stat.histf(1:evalcount))
hold on
plot([1:evalcount], C_nn * ones(1,evalcount), '-.r')
grid on
hold off
title('Tour length')

subplot(2,2,2)
plot_path(map, popt, [func2str(map_fct), '.png'])
title('Currently optimal tour')

subplot(2,2,3)
image(256 * ((pheromones - min(min(pheromones))) / (max(max(pheromones)) - min(min(pheromones)))));
title('Pheromone matrix')

subplot(2,2,4)
image(256 * ((heuristics - min(min(heuristics))) / (max(max(heuristics)) - min(min(heuristics)))));
title('Heuristic desirability matrix')

drawnow()

end

end

end
```