The traveling salesman problem (TSP) asks the following question: “Given a list of n cities, among which an origin city, and the distances between each pair of cities, what is the shortest possible route that visits each city exactly once and returns to the origin city?” This is an important NP-hard combinatorial optimization problem, particularly in the fields of operations research and theoretical computer science. The problem was first formulated in 1930 and is one of the most intensively studied problems in discrete optimization.
We are going to solve this problem using different approaches : in this post, we will use heuristics, and in an upcoming post we will improve our solutions using metaheuristics.
1. Exact methods
Many optimization problems can be solved by a linear programming approach, which is amazing because linear problems (LP) can be solved in a short amount of time, even for huge instances. And this is not limited to continuous variable problems : under mathematical conditions, we can perform a constraint relaxation and solve integer linear program (ILP) with for instance the simplex algorithm, and still obtain the optimal solution.
Unfortunatly for us, altough the TSP problem does have a linear formulation, we can’t solve it optimally with a continuous constraint relaxation method.
And it is now that the real difficulties begin to appear, because the naive solution’s complexity is 
(that’s 77 years to solve an instance with 20 cities with a frequency of 1GHz).
The DP (dynamic programming) version algorithm (Bellman-Held-Karp algorithm) will have the complexity of 
. By reducing the complexity from factorial to exponential, if the size of n is relatively small, the problem can be solvable (this time we can solve a 50 cities-sized instance in only 90 years… maybe we can do better 😅).
So instead of searching an exact solution, we will try to find an approximate one in a much shorter timespan. I am talking about solving a instance with 400 cities in less than 30 seconds. Yes, that’s quite the improvement !
2. Heuristics
Heuristics are simple methods used to solve the problem. They are simple to understand, and produce an approximate solution in a really short time. In fact, even humans use heuristics for decision making in real life. We will study four of them : first fit, best fit, nearest point and a special ‘circle’ heuristic.
We will use a simple structure to represent our solutions : a list of length $n$ which stores the order in which we visit the different cities.
Helper functions
import csv
import numpy as np
import matplotlib.pyplot as plt
import time
import random
from math import atan2
def loadFixedInstance(fileName):
folderName = "Instances/"
if not (".csv" in fileName):
fileName += ".csv"
if not (folderName in fileName):
fileName = folderName + fileName
with open(fileName) as file:
lines = csv.reader(file)
instance = []
for line in lines:
instance.append((int(line[0]), int(line[1])))
return instance
def computeDistanceBetweenCities(city1, city2):
return np.linalg.norm(np.array(city2) - np.array(city1), ord = 2)
def possiblePermutations(x):
n = len(x)
perm = []
for i in range(n):
for j in range(i+1,n):
perm.append((i,j))
return perm
def combinatorialConverter(x):
return [elt[0] for elt in sorted(zip(range(len(x)),x), key=lambda t: t[1], reverse=False)]
def distMatrix(cities):
M=[]
n = len(cities)
for i in range(n):
L = []
for j in range(n):
L.append(computeDistanceBetweenCities(cities[i],cities[j]))
M.append(L)
return M
def fcomb(x, distMat):
routeDistance = 0
routeDistance += distMat[x[-1]][x[0]]
for i in range(len(x)-1):
routeDistance += distMat[x[i]][x[i+1]]
return routeDistance
def permuteElements(x,i,j,copy=False):
if copy == True:
x_ = x.copy()
x_[i],x_[j] = x_[j],x_[i]
return x_
x[i],x[j] = x[j],x[i]
return x
def plotRoute(route):
xs = [instance[i][0] for i in route] + [instance[route[0]][0]]
ys = [instance[i][1] for i in route] + [instance[route[0]][1]]
plt.plot(xs, ys, 'o-')
plt.ylabel('x')
plt.xlabel('y')
plt.show()
We write some code that will be used by our algorithms :
- a function to build a distance matrix between the cities
- an evaluation function to check solution quality
- a function to plot the results
Hill climbing - First Fit
def stoch_descent(x,permutations,distMat):
n = len(permutations)
si = list(range(n))
score = fcomb(x,distMat)
random.shuffle(si)
i = 0
while i<n:
x_ = x.copy()
x_ = permuteElements(x_,*permutations[si[i]])
i+=1
score_ = fcomb(x_,distMat)
if score_<score:
improvement = True
x = x_
score = score_
i = 0
return x,score
instance = loadFixedInstance("TSPInstance50.csv")
permutations = possiblePermutations([0]*len(instance))
distMat = distMatrix(instance)
initx = list(range(len(instance)))
ptime = time.time()
bestx, bestScore = stoch_descent(initx,permutations,distMat)
print(f'Minimal cost found : {round(bestScore,1)} in {round(time.time()-ptime,1)} s')
plotRoute(bestx)
Minimal cost found : 1491.8 in 1.5 s
Here, we loop over every possible permutation that can be made after having shuffled the permutation list (result quality may vary). Every time we improve the solution, we start again from the beggining of the permutation list. If we loop through the entire loop without seing any improvement, the algorithm is terminated.
With this algorithm, we can find a solution in one second, with a mean score of 1530 over 20 solutions. The best score was 1333 and the worst was 1815. The score are not optimal because the algorithm was stuck in a local minimum. Running the algorithm multiple times is a way to find a better quality solution.
Best Fit
def complete_descent(x,permutations,distMat):
continueLoop = True
n = len(permutations)
score = fcomb(x,distMat)
while continueLoop:
scores = [fcomb(permuteElements(x,*permutations[i],copy=True),distMat) for i in range(n)]
bestidx = np.argmin(scores)
if scores[bestidx]<score:
score = scores[bestidx]
x = permuteElements(x,*permutations[bestidx])
else:
continueLoop = False
return x, score
initx = list(range(len(instance)))
ptime = time.time()
bestx, bestScore = complete_descent(initx,permutations,distMat)
print(f'Minimal cost found : {round(bestScore,1)} in {round(time.time()-ptime,1)} s')
plotRoute(bestx)
Minimal cost found : 1589.5 in 2.8 s
This time, we solve the problem using a best fit approach. At each iteration, we choose the permutation which improves the solution the most, until we observe no further improvement.
As expected, this algorithm is a little bit slower, and the end result is not particularly great… maybe some randomness is neeeded 😀
Nearest city
def nearestHeuristic(dist_mat):
x = [0]
n = len(distMat)
for _ in range(n-1):
nextCities = [elt[0] for elt in sorted(zip(list(range(n)),dist_mat[x[-1]]), key = lambda x:x[1]) if elt[0] not in x]
x.append(nextCities[0])
return x, fcomb(x,dist_mat)
ptime = time.time()
bestx, bestScore = nearestHeuristic(distMat)
print(f'Minimal cost found : {round(bestScore,1)} in {round(time.time()-ptime,3)} s')
plotRoute(bestx)
Minimal cost found : 1351.3 in 0.006 s
Here we use a greedy algorithm by always choosing the nearest available city. This is much better than previous solutions : in only 6ms, the algorithm was able to find a solution that beats the previous ones ! But still we can observe some crossing, which means that the solution is still sub-optimal
Circle heuristic
def circleHeuristic(cities):
minx = min(cities,key=lambda x:x[0])[0]
maxx = max(cities,key=lambda x:x[0])[0]
miny = min(cities,key=lambda x:x[1])[1]
maxy = max(cities,key=lambda x:x[1])[1]
center = ((minx+maxx)/2,(miny+maxy)/2)
angles = []
for city in cities:
angles.append(atan2(city[1]-center[1],city[0]-center[0]))
return combinatorialConverter(angles)
ptime = time.time()
bestx = circleHeuristic(instance)
bestScore = fcomb(bestx, distMat)
print(f'Minimal cost found : {round(bestScore,1)} in {round(time.time()-ptime,3)} s')
plotRoute(bestx)
Minimal cost found : 1667.6 in 0.001 s
This is once again a very interesting approach, here we project every city on a circle and then join them in a circular way
Can we do better ?
We saw that using simple heuristics, we can find approximate solution to the TSP in a very short amount of time. However, the quality of the solutions is oftentimes disappointing. To find better solutions, we can do multiple things. Firstly, we can combine heuristics : For instance, we can build an initial solution with our circle heuristic, an then improve the solution with a first fit hill climbing algorithm.
initx = circleHeuristic(instance)
ptime = time.time()
bestx, bestScore = complete_descent(initx,permutations,distMat)
print(f'Minimal cost found : {round(bestScore,1)} in {round(time.time()-ptime,1)} s')
plotRoute(bestx)
Minimal cost found : 1329.7 in 1.1 s
Otherwise, we can use other approximate algorithms, such as the Lin-Kernighan algorithm (LK), which uses an operator called λ-opt to modify the solution into a better one. Finally, in the next post we will implement an algorithm called Simulated Annealing, that is a type of metaheuristic algorithm.
Credits : Photo by CHUTTERSNAP on Unsplash