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

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:

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)\]

- 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).\]
- If \(\Delta E \leq 0\), jump to \(\omega'\) since it’s better (it’s a shorter route).
- 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.

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