The problem

The traveling salesman problem is an example of a combinatorial optimization problem where the goal is to minimize the distance of a tour. A tour is a route through all the cities (represented as points in the plane) in a given set, visiting each city once and returning to the start. The distance between cities is Euclidean distance in the plane. This is a natural problem that, for example, shipping companies face when deciding how to route their drivers as they make deliveries. The issue is that it’s in a class of problems for which algorithms are provably hard to design.

If we have \(n\) cities or vertices to visit, then a tour can be thought of mathematically as a permutation of \(\{1,\ldots,n\}\). The following code defines a random set of points to be the cities and plots a tour chosen uniformly at random.

n <- 100
x <- rnorm(n = n, mean = 0, sd = 10)
y <- rnorm(n = n, mean = 0, sd = 10)
tour <- sample(x = 1:n, size = n, replace = F)
route <- rbind(c(x[tour], x[tour[1]]), c(y[tour], y[tour[1]]))
qplot(x = route[1,], y = route[2,], geom = "path", xlab = "", ylab = "")

An algorithm

Our goal is to search randomly for the tour which minimizes Euclidean distance. Picking uniformly at random among the set of tours is hopeless since it’s a massive set. Instead, we’ll construct a random process which jumps through the set of possible tours one at a time in a hopefully more efficient way. Our algorithm is defined as follows:

  1. Suppose we are currently at tour \(\omega\). Choose two cities from \(\omega\) at random and interchange them. That is, if our current tour is \[\omega = (a_1,\ldots,a_{i-1},a_i,a_{i+1},\ldots,a_{j-1},a_j,a_{j+1},\ldots,a_n)\] and we choose \(i<j\) to interchange, then we propose jumping to \[\omega' = (a_1,\ldots,a_{i-1},a_j,a_{i+1},\ldots,a_{j-1},a_i,a_{j+1},\ldots,a_n)\]

  2. Let \(E(\omega)\) and \(E(\omega')\) denote the total distance traveled on tour \(\omega\) and \(\omega'\) respectively and let \[\Delta E = E(\omega')-E(\omega).\]
    1. If \(\Delta E \leq 0\), jump to \(\omega'\) since it’s better (it’s a shorter route).
    2. If \(\Delta E > 0\), either jump to \(\omega'\) or stay at \(\omega\). The probability of jumping is \[e^{-\Delta E/(kT)}\] where \(k, T > 0\) are given parameters (called the Boltzmann constant and the temperature respectively; we typically suppose \(k = 1\)). The probability of staying is the complementary probability.
  3. Repeat.


The following code outlines this process.

N <- 10000 # number of steps of chain
Temperature <- 1
Ehistory <- array(0, N)
tours <- array(0,c(N+1, n))
tours[1,] <- tour

for (k in 1:N) {
  swap <- sample(1:n, size = 2) # choose two cities to swap (this is i and j)
  i <- min(swap)
  j <- max(swap)
  proposedTour <- tours[k,] # make a copy of current tour
  proposedTour[c(i,j)] <- proposedTour[c(j,i)] # swap the two cities in our current tour
  oldRoute <- rbind(c(x[tours[k,]], x[tours[k,1]]), c(y[tours[k,]], y[tours[k,1]])) # xy-coordinates of current tour
  proposedRoute <- rbind(c(x[proposedTour], x[proposedTour[1]]), c(y[proposedTour], y[proposedTour[1]])) # xy-coordinates of proposed tour
  oldE <- 0
  proposedE <- 0
  for (l in 1:n) { # find total distance of the two tours
    oldE <- oldE + dist(oldRoute[1:2,l:(l+1)], method = "euclidean")
    proposedE <- proposedE + dist(proposedRoute[1:2,l:(l+1)], method = "euclidean")
  DeltaE <- proposedE - oldE
  if (DeltaE <= 0) { # if proposed tour is shorter accept
    tours[k + 1,] <- proposedTour
    Ehistory[k] <- proposedE
  else {
    if (runif(1) < exp(-DeltaE/Temperature)) {
      tours[k + 1,] <- proposedTour
      Ehistory[k] <- proposedE
    else {
      tours[k + 1,] <- tours[k,]
     Ehistory[k] <- oldE

qplot(x = 1:N, y = Ehistory, xlab = "step", ylab = "tour length")