### 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.

set.seed(339)
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.

### Code

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")