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.
In the previous post, we used heuristics to solve the TSP problem. In this post we will improve our solutions using metaheuristics.
1. Metaheuristics
As we saw in the previous post, computing optimal solutions is intractable for many optimization problems of industrial and scientific importance. In practice, we must find alternative methods such as heuristic or metaheuristic algorithms and be satisfied with “good” resulting solutions.
Metaheuristic search methods can be defined as “upper level” general methodologies used to solve specific optimization problems. There are two common design questions related to all iterative metaheuristics: the representation of solutions handled by algorithms and the definition of the objective function that will guide the search.
Representation
When designing a metaheuristics, we need to think of a way to encode a solution. Although many representations may exists for a specific problem, a representation must follow common characteristics :
- Completeness: All solutions associated with the problem must be represented.
- Connexity: A search path must exist between any two solutions of the search space, such that any solution can be attained.
- Efficiency: The representation must be easy to manipulate by the search operators. The time and space complexities of the operators dealing with the representation must be reduced.
For the TSP problem, we will use a simple list encoding, with the i-th element of the list being the i-th city to be visited, in order.
Cost function
The objective function or cost function formulates the goal to achieve. It associates with each solution of the search space a real value that describes the quality or the fitness of the solution. It is an important element in designing a metaheuristic. It will guide the search toward “good” solutions of the search space. If the objective function is improperly defined, it can lead to nonacceptable solutions whatever metaheuristic is used.
For instance, we will use the following objective function :
where π represents a permutation encoding a tour and n the number of cities.
There exist two families of metaheuritics algorithms : population-based, where a family of solution is generated and modified during the optimization process, and trajectory-based, were the initial solution is modified. Here we will use the Simulated Annealing algorithm, which is a simple trajectory-based algorithm.
Metaheuristics work in a similar way :
- At first, an initial solution (or family of solutions) is generated
- During the optimization, the algorithm “navigates” in the search space by modifying the solution with operators, while trying to avoid local optima.
At the end of the optimization process, the best solution is returned.
Helper functions
import csv
import numpy as np
import matplotlib.pyplot as plt
plt.style.use('seaborn')
import time
import random
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 combinatorialConverter(x):
return [elt[0] for elt in sorted(zip(range(len(x)),x), key=lambda t: t[1], reverse=False)]
def f(x):
L = combinatorialConverter(x)
routeDistance = 0
routeDistance += computeDistanceBetweenCities(instance[L[-1]],instance[L[0]])
for i in range(len(x)-1):
routeDistance += computeDistanceBetweenCities(instance[L[i]],instance[L[i+1]])
return routeDistance
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()
#------Optimization functions---------#
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 sample_random(end):
"""sample 2 different numbers"""
i = random.randrange(0,end)
j = random.randrange(0,end)
while i==j:
i = random.randrange(0,end)
j = random.randrange(0,end)
return i,j
We use the same code from the last post :
- a function to build a distance matrix between the cities
- an evaluation function to check solution quality
- a function to plot the results
Search operators
def twoOpt(x,i,j):
"""perform permutation
x : current solution
i<j is necessary"""
if i==0:
return x[:i]+x[j::-1]+x[j+1:]
return x[:i]+x[j:i-1:-1]+x[j+1:]
def randomPermutation(x,n,matDist):
"""Find neighbour using 2-opt exchange operator and return its score
x : current solution
n : len of x"""
i,j = sample_random(n)
if i>j:
i,j=j,i
oldScore = dScore(x,n,i,j,matDist)
x = twoOpt(x,i,j)
deltaScore = modifiedObjFunc(x,n,i,j,oldScore,matDist)
return x, deltaScore, i, j
How to navigate trough the search space ? We use operators that modify the solution by returning a “neighbour” solution. We must then define what is a neighbour solution for the TSP solution.
A straightforward approach consists in choosing two cities and changing their position in the list. Is is called the swap operator. The problem however with this operator is that we may observe a large variation generated depending on the cities chosen (called weak locality).
A better approach is the λ-opt operator, and especially the 2-opt operator which consists in swapping edges instead of cities. This operator ensures a strong locality.
Objective function
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 dScore(x,n,i,j,matDist):
s=0
if i>0:
s+=matDist[x[i-1],x[i]]
else:
#i=0, need to update distance between x[0] and x[-1]
s+=matDist[x[-1]][x[0]]
if j<n-1:
s+=matDist[x[j],x[j+1]]
else:
#j=n-1, need to update distance between x[0] and x[-1]
s+=matDist[x[-1]][x[0]]
return s
def modifiedObjFunc(x,n,i,j,oldScore,matDist):
"""compute the delta of score in an optimized way
OldScore : sum of the distance between x[i-1]&x[i] and between x[j]&x[j+1]"""
return dScore(x,n,i,j,matDist)-oldScore
We can use the simple function formulation as sum of the length of the cycle. However, if we think about it, at each iteration, only a part of the solution is modified, thus we could only compute the difference in the cost induced by the search operator. That is exactly what the function dScore is doing.
Simulated Annealing
def SA(x0, iTemp, fTemp, nPerm, cRate, distMat):
"""solve tsp instance using optimized simulated annealing"""
n = len(x0)
T = iTemp
niter = 0
SCORES = []
x = x0.copy()
score = fcomb(x, distMat)
best_score=score
bestx = x0.copy()
ptime = time.time()
titer = np.log(fTemp/iTemp)/np.log(cRate)
while T>fTemp:
# if niter%20000==0:
# print(f'Time elapsed : {round(time.time()-ptime,1)}s, {round(100*(niter/nPerm)/titer,1)} %')
for _ in range(nPerm):
x, deltaScore, i, j = randomPermutation(x,n,distMat)
if deltaScore<0 or random.random()<np.exp(-deltaScore/T):
score+=deltaScore #we need the true score for best solution calculation
if score<best_score:
bestx = x.copy()
best_score = score
else:
#reverse operator
x = twoOpt(x,i,j)
niter+=1
SCORES.append(score)
T*=cRate
return bestx, best_score, SCORES
Finally, we can implement the optimization procedure. In the early 1980s three IBM researchers, Kirkpatrick et al., introduced the concepts of annealing in combinatorial optimization. These concepts are based on a strong analogy with the physical annealing of materials. This process involves bringing a solid to a low energy state after raising its temperature. In the context of the optimization algorithm, the function to be minimized represents the energy of the solid.
The procedure goes as follows :
Initialization (i := istart, k := 0, ck = c0, Lk := L0)
while not(ck≃0):
for l = 0 to Lk do:
Generate a solution j from the neighborhood Si of the current solution i
If f(j) < f(i) then i := j (j becomes the current solution)
Else, j becomes the current solution with probability exp[(f(i)−f(j))/Tk]
k := k+1
Compute(Lk, ck)
When the temperature is high, the algorithm performs an exploration phase to find promising parts of the solution space (the probability of accepting a non-improving solution is around 1).
When the temperature is low, the algorithm intensifies its search, and non-improving solutions are less likely to be accepted.
At each temperature stage, we perform a fixed number of iteration to sufficiently navigate the search space, and we decrease the temperature with a geometric cooling scheme.
How to choose the initial temperature ?
def ITAcceptance(solution,acceptance_rate,n_samples,matDist):
"""compute initial temperature such that a certain acceptance rate is achieved
T->0 if R->0.5 and T->oo if R->1"""
m1 = 0
m2 = 0
DF = 0
n = len(solution)
for _ in range(n_samples):
solution, deltaScore, i, j = randomPermutation(solution,n,matDist)
solution = twoOpt(solution,i,j)
if deltaScore<0:
DF-=deltaScore
m1+=1
else:
m2+=1
mDF = DF/m1
return mDF/np.log(m2/(abs(m2*acceptance_rate-m1*(1-acceptance_rate))))
It can be difficult to choose an initial temperature, as it depends on the type of problem and the objective function. Thus, we choose the initial temperature such that X% of the initial neighbours are accepted. In fact, we simulate the algorithm a certain number of time, and compute the initial temperature needed to meet the “acceptance requirements”. For the experiments, the ending temperature is set as T0/1000.
instance = loadFixedInstance("TSPInstance150.csv")
matDist = np.array(distMatrix(instance))
initx = list(range(len(instance)))
T0 = ITAcceptance(initx,0.8,10000,matDist)
print(f'T0 : {round(T0,1)}')
ptime = time.time()
bestx, best_score, SCORES = SA(initx, T0, T0/1000, 500, 0.995, matDist)
print(f'Optimization done in {round(time.time()-ptime,1)} s, the best score was : {round(best_score,1)}')
plt.plot(list(range(len(SCORES))),SCORES)
plt.show()
T0 : 120.7
Optimization done in 25.4 s, the best score was : 1908.2
plotRoute(bestx)
In only 25 seconds we were able to find a pretty good solution with no overlapping trajectories. We can tweak the parameters to gain even more speed (as we can see the optimization was almost finished at 800 temperature iterations).
Let’s try the instance with 400 cities !
instance = loadFixedInstance("TSPInstance400.csv")
matDist = np.array(distMatrix(instance))
initx = list(range(len(instance)))
T0 = ITAcceptance(initx,0.8,10000,matDist)
print(f'T0 : {round(T0,1)}')
ptime = time.time()
bestx, best_score, SCORES = SA(initx, T0, T0/1000, 500, 0.996, matDist)
print(f'Optimization done in {round(time.time()-ptime,1)} s, the best score was : {round(best_score,1)}')
plt.plot(list(range(len(SCORES))),SCORES)
plt.show()
T0 : 126.6
Optimization done in 28.7 s, the best score was : 3184.2
plotRoute(bestx)
We did it ! The instance was solved in less than 30 seconds. This experiment shows that, using the right search operators and an optimized objective function, the simulated annealing algorithm is a great optimization algorithm.
And then ?
In this tutorial we saw that using metaheuristics enable us to solve NP-hard cominbatorial problems in a reasonable time. Here we used the Simulated Annealing algorithm but many, many more algorithms do exist. If you are interested in these other algorithms, I recommand the Handbook of Metaheuristics by Springer (third version as of may 2022, can be found online) which contains a complete overview of classical metaheuristics algorithms.
Source :
- Talbi (2009). Metaheuristics - From Design To Implementation. Wiley.
- Gendreau, Potvin (2019). Handbook of Metaheuristics (3rd edition). Springer.
Credits : Photo by CHUTTERSNAP on Unsplash