Traveling Salesperson Problem

Setup

First we generate a bunch of points in the unit square [0,1]×[0,1].

numPoints=7;
points=RandomReal[1,{numPoints,2}]

It’s helpful to plot the points.

ListPlot[points->Range[numPoints]]

Let a tour be represented by a list of indexes that specify the order in which we visit the points. The simplest tour is to visit the points in the order that they appear in our list. We’ll call this the trivial tour, and it is probably not close to optimal.

trivialTour=Range[numPoints]

Previously, we wrote a module that draws the tour:

plotTour[points_,tour_]:=Module[
{n=Length[points],tourPts},
tourPts=Table[0,n+1]; (*initialize list of size n+1*)

(* fill in the cities in order in tourPts *)
Do[
tourPts[[i]]=points[[tour[[i]]]],
{i,1,n}];
tourPts[[n+1]]=points[[tour[[1]]]] ;(* return to starting point *)

(* return the plot *)
ListPlot[tourPts,Joined->True]
]

plotTour[points,trivialTour]

We also wrote a module to compute the total distance along the tour.

distTour[points_,tour_]:=Module[
{dist=0},

(* loop over all edges along the tour *)
Do[
dist+=EuclideanDistance[points[[tour[[i]]]],points[[tour[[i+1]]]]],
{i,1,Length[tour]-1}];
dist+=EuclideanDistance[points[[tour[[-1]]]],points[[tour[[1]]]]]; (* return to starting point*)

(* return the distance *)
dist
]

distTour[points,trivialTour]

Exhaustively consider all tours

Suppose we start at point 1. Here are all the ways of visiting the other nine points:

Permutations[Range[2,numPoints]]

Exhaustively compute distances for all tours

Last time, we wrote a module that finds the minium-length tour through a given set of points. Since this module uses an exhaustive search, it’s only feasible for small point sets (say ten or fewer points).

exhaustiveMinDist[points_]:=Module[
{n=Length[points],
perms, (* this variable will hold all permutations of {2, 3, ..., n} *)
currTour, (* the current tour *)
currDist,(* length of current tour *)
minTour, (* min tour found so far *)
minDist=∞ (* min distance found so far; initialize to a large value *)
},

(* compute all permutations *)
perms=Permutations[Range[2,n]];

(* loop over all tours *)
Do[
currTour=Prepend[p,1];
currDist=distTour[points,currTour];
(*Print["tour ",currTour,"; distance ",currDist];*)
If[currDist<minDist,minDist=currDist;minTour=currTour]
,{p,perms}];

(* return the minimum distance*)
{minDist,minTour}
]

Test the previous module:

numPoints=5;
points=RandomReal[1,{numPoints,2}];
{minDist,minTour}=exhaustiveMinDist[points]
plotTour[points,minTour]

Average length of a tour

What is the average length of a tour of n random points? How does this depend on n? For how large a value of n can we estimate this average?

The following module estimates the average length of a tour through numPoints points selected uniformly at random from the square [0,1]×[0,1]. It does this by generating numTrials sets of points and exhaustively finding the minimum-length tour for each.

avgDist[numPoints_,numTrials_]:=Module[
{points,distList=Table[0,numTrials]},
Do[
points=RandomReal[1,{numPoints,2}];
distList[[i]]=exhaustiveMinDist[points][[1]]
,{i,numTrials}];
(*Print[distList];*)
Mean[distList]
]

Try it out:

avgDist[8,20]

Estimate the minimum tour length for up to ten points:

avgVals=Table[avgDist[n,50],{n,2,10}]

Plot the tour lengths. Does the graph resemble a square root function?

Show[
ListPlot[avgVals,DataRange->{2,10},PlotRange->{0,3.3},AxesLabel->{"number of points","min tour distance"}],
Plot[Sqrt[n],{n,2,10},PlotStyle->Orange]
]

Simulated  Annealing

Let’s use simulated annealing to find (approximately) minimum tours for larger values of n.

Objective Function

We want to minimize the length of the tour. We already implemented our distTour function that computes this length.

Transition  function

One option is to simply swap two indexes in the tour:

proposeMoveSwap[currTour_]:=Module[
{n=Length[currTour],
propTour=currTour},




(* return the proposed tour *)
propTour
]

Unfortunately, we can get into situations where any individual swap will lengthen the tour, but some combination of simultaneous swaps will shorten the tour. In this case, it can be helpful to reverse a randomly-selected subsequence of the tour. (See page 312 of our Computational Mathematics text for a diagram.)

proposeMoveReverse[currTour_]:=Module[
{n=Length[currTour],
propTour=currTour},



(* return the proposed tour *)
propTour
]

Testing:

The simulated annealing algorithm

Now we can implement the full simulated annealing algorithm. The following module requires the following parameters:

a list of points

a function that returns a proposal transition

a starting value of the β parameter

the multiplicative increment for the β parameter

the number of iterations to compute

findMinTour[points_,proposeMove_,βstart_,βinc_,numIter_]:=Module[
{n=Length[points], (* number of points *)
currTour, currDist,
β=βstart, (* value of β we can update *)
propTour, (* variables for generating a new state *)
currF,propF,Δf ,fVals=Table[0,numIter],(* values of the objective function *),
ρ (* acceptance ratio *)
},

 (* start with a trivial tour *)
currTour=Range[n];

(* compute the starting tour distance *)
currDist=distTour[points,currTour];
fVals[[1]]=currDist;

(* main loop *)
Do[
(* generate a proposal tour *)
propTour=proposeMove[currTour];

(* compute change in the tour distance *) 
Δf=distTour[points,currTour]-distTour[points,propTour];

(* select the next β in the annealing schedule *)
β*=βinc;

(* compute the acceptance ratio *) 
ρ=Exp[β*Δf];

(* decide whether to accept the proposed transition *) 
If[RandomReal[]<ρ,
currTour=propTour;
];

(* store the tour distance *)
currDist=distTour[points,currTour];
fVals[[i]]=currDist;
,{i,2,numIter}];

(* return values *) 
{fVals,β,currTour,currDist}
]

Try the simulated annealing algorithm, first on small sets of points, then on larger sets of points.