TSP Essay

A fun study of some heuristics for the Travelling Salesman Problem.

---

I've always found the Travelling salesman problem particularly interesting. After studying a few papers on the subject (in this repository too), mainly the superb Clarke & Wright's Savings Algorithm from 1997, I felt like trying a few techniques.

To see the original ipynb format (probably better): https://github.com/rsalmei/tsp-essay/blob/master/tsp_essay.ipynb

To run this yourself, you need:

• a recent python (I'm using 3.6.5)
• a virtualenv
• pip install -r requirements.txt
• jupyter notebook tsp_essay.ipynb

Enjoy :)

Setup environment

import pandas as pd
import numpy as np
%matplotlib inline %config InlineBackend.figure_format = 'retina' import matplotlib.pyplot as plt from pylab import rcParams

/Users/rogerio/.pyenv/versions/3.6.5/lib/python3.6/importlib/_bootstrap.py:219: RuntimeWarning: numpy.dtype size changed, may indicate binary incompatibility. Expected 96, got 88 return f(args, *kwds) /Users/rogerio/.pyenv/versions/3.6.5/lib/python3.6/importlib/_bootstrap.py:219: RuntimeWarning: numpy.dtype size changed, may indicate binary incompatibility. Expected 96, got 88 return f(args, *kwds)

Customize graphics style

print(plt.style.available)

['seaborn-dark', 'seaborn-darkgrid', 'seaborn-ticks', 'fivethirtyeight', 'seaborn-whitegrid', 'classic', 'classictest', 'fast', 'seaborn-talk', 'seaborn-dark-palette', 'seaborn-bright', 'seaborn-pastel', 'grayscale', 'seaborn-notebook', 'ggplot', 'seaborn-colorblind', 'seaborn-muted', 'seaborn', 'SolarizeLight2', 'seaborn-paper', 'bmh', 'tableau-colorblind10', 'seaborn-white', 'darkbackground', 'seaborn-poster', 'seaborn-deep']

plt.style.use('classic')
rcParams['figure.figsize'] = 16, 10

Initialize data

np.random.seed(18)
size, cmin, cmax = 150, -100, 100 data = pd.DataFrame(dict(x=[0], y=[0])).append(     pd.DataFrame((np.random.randomsample(2size)(cmax-cmin)+cmin).reshape(-1,2), columns=['x', 'y']), ignoreindex=True) data

Plot problem

plt.scatter(data.x, data.y, c='r', marker='o', s=40)
plt.scatter(data.iloc[0].x, data.iloc[0].y, marker='*', c='r', s=400)
plt.axis([-100,100,-100,100])
plt.grid(True)

Compute distance matrix

distance_func = lambda a,b: np.sqrt((a.x-b.x)2 + (a.y-b.y)2)
def computedistances(dfunc=distancefunc):
    for i in range(size+1):
        current=data.iloc[i]
        data[i]=dfunc(current, data)
        #data[str(i)][i]=np.nan
%time compute_distances() data

CPU times: user 346 ms, sys: 5.93 ms, total: 352 ms Wall time: 357 ms

Helper functions

estimate_cost = lambda route: sum(data.iloc[i][j] for i,j in zip(route, route[1:]))
get_coords = lambda route: list(zip(*[(data.iloc[i].x,data.iloc[i].y) for i in route]))

Scenario: center depot, one van, nearest neighbor heuristic

Apply a greedy NN algorithm

def nearestneighbor(unserved, stopcondition=lambda route, target: False):
    current=0 #depot
    result_path=[]
    while True:
        result_path.append(current)
        unserved.remove(current)
        if not unserved:
            break
current=data.iloc[unserved,current+2].idxmin()         if stopcondition(resultpath, int(current)):             if len(result_path)>1:                 break
result_path.append(0)     return result_path
%time route = nearest_neighbor(list(range(size+1))) print('cost={}\nroute={}'.format(estimate_cost(route),route))

CPU times: user 116 ms, sys: 8.29 ms, total: 125 ms Wall time: 122 ms cost=2623.828598616268 route=[0, 53, 60, 20, 66, 148, 55, 51, 1, 120, 104, 123, 74, 93, 130, 42, 29, 19, 91, 107, 81, 11, 134, 117, 57, 30, 34, 67, 75, 46, 39, 131, 16, 114, 72, 96, 38, 22, 122, 92, 36, 126, 113, 7, 58, 111, 28, 108, 63, 2, 137, 50, 128, 69, 141, 18, 146, 37, 54, 95, 24, 121, 35, 83, 3, 49, 12, 73, 105, 129, 90, 144, 9, 86, 132, 31, 70, 127, 149, 6, 71, 41, 140, 76, 102, 32, 100, 97, 139, 88, 78, 26, 17, 135, 27, 13, 110, 89, 59, 119, 40, 15, 56, 145, 10, 106, 82, 94, 115, 87, 98, 44, 64, 52, 14, 33, 136, 124, 65, 142, 150, 147, 143, 109, 112, 23, 47, 101, 45, 138, 99, 133, 8, 48, 62, 80, 103, 43, 68, 61, 125, 21, 25, 118, 5, 79, 116, 77, 85, 84, 4, 0]

Plot result path

coords=get_coords(route)
plt.plot(coords[0], coords[1], 'b-', linewidth=2)
plt.scatter(data.x, data.y, c='r', s=30)
plt.scatter(data.iloc[0].x, data.iloc[0].y, marker='*', c='r', s=400)
plt.axis([-100,100,-100,100])
plt.grid(True)

Let's improve it with a 2-opt heuristic

def twoOptSwap(route,i,j):
    t=route[i:j+1];
    t.reverse();
    route[i:j+1]=t
def twoOpt(route):     cost=estimate_cost(route)
def inner():         nonlocal cost         for i in range(1,len(route)-2):             nodei=data.iloc[route[i]]             for j in range(i+1,len(route)-1):                 nodej=data.iloc[route[j]]                 save=nodei.iloc[route[i-1]+2]+nodej.iloc[route[j+1]+2] - (nodej.iloc[route[i-1]+2]+nodei.iloc[route[j+1]+2])                 if save>0:                     twoOptSwap(route,i,j)
cost-=save                     print('exchanging {}-{},{}-{} with {}-{},{}-{} => save={} => cost: {}'.format(                             route[i-1],route[i],route[j],route[j+1],route[i-1],route[j],route[i],route[j+1],save,cost))                     return False         return True
while True:         if inner():             break
%time twoOpt(route)

exchanging 53-64,60-52 with 53-60,64-52 => save=2.8803537373675177 => cost: 2620.948244878901 ... exchanging 1-120,81-107 with 1-81,120-107 => save=26.496627104257605 => cost: 2100.7350743873462 CPU times: user 3min 11s, sys: 1.65 s, total: 3min 13s Wall time: 3min 21s

Plot improved path

coords=get_coords(route)
plt.plot(coords[0], coords[1], 'b-', linewidth=2)
plt.scatter(data.x, data.y, c='r', s=30)
plt.scatter(data.iloc[0].x, data.iloc[0].y, marker='*', c='r', s=400)
plt.axis([-100,100,-100,100])
plt.grid(True)

Incredible! A clean path, without any crossing parts!

Scenario: center depot, one van, 2-opt heuristic

Just for fun, let's pass a random route to 2-opt heuristic

route = [x for x in range(size+1)]+[0]
print('cost={}\nroute={}'.format(estimate_cost(route),route))

cost=16332.039263088578 route=[0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 127, 128, 129, 130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142, 143, 144, 145, 146, 147, 148, 149, 150, 0]

Plot result path

coords=get_coords(route)
plt.plot(coords[0], coords[1], 'b-', linewidth=2)
plt.scatter(data.x, data.y, c='r', s=30)
plt.scatter(data.iloc[0].x, data.iloc[0].y, marker='*', c='r', s=400)
plt.axis([-100,100,-100,100])
plt.grid(True)

Let's see what 2-opt heuristic can do

%time twoOpt(route)

exchanging 0-4,1-5 with 0-1,4-5 => save=11.966404491447804 => cost: 16320.07285859713 ... exchanging 47-23,112-8 with 47-112,23-8 => save=11.556726776627453 => cost: 2181.5435075388827 CPU times: user 19min 11s, sys: 8.18 s, total: 19min 19s Wall time: 19min 44s

Plot improved path

coords=get_coords(route)
plt.plot(coords[0], coords[1], 'b-', linewidth=2)
plt.scatter(data.x, data.y, c='r', s=30)
plt.scatter(data.iloc[0].x, data.iloc[0].y, marker='*', c='r', s=400)
plt.axis([-100,100,-100,100])
plt.grid(True)

WOW, very cool, a clean path too, almost as good as the previous one.

But extremely slower... It took ~19 minutes, while NN + 2-opt took only ~3 minutes.

Scenario: center depot, several vans, no clustering

Let's apply a sequential NN algorithm, but including a max cost constraint

def sequentialnn(maxcost):
    unserved=list(range(size+1))
    result=[]
    while True:
        resultpath = nearestneighbor(unserved, lambda route, target: estimatecost(route+[target,0])>=maxcost)
        result.append(result_path)
        print('van {}: {}'.format(len(result), estimatecost(resultpath)))
if not unserved:             break         unserved.append(0)     return result
%time resultpaths = sequentialnn(400)

van 1: 389.0742611200774 van 2: 398.3533734081084 van 3: 397.4195886051539 van 4: 390.9917686566429 van 5: 371.6762623788105 van 6: 320.29683749922486 van 7: 364.8394134778564 van 8: 320.99556479051967 van 9: 210.1189842618746 van 10: 276.59421133431306 CPU times: user 1.48 s, sys: 24.4 ms, total: 1.51 s Wall time: 1.54 s

Plot result paths

for path in result_paths:
    coords=get_coords(path)
    plt.plot(coords[0], coords[1], linewidth=2)
plt.scatter(data.x, data.y, c='r', s=30)
plt.scatter(data.iloc[0].x, data.iloc[0].y, marker='*', c='r', s=400)
plt.axis([-100,100,-100,100])
plt.grid(True)

Let's improve them with a 2-opt heuristic

def improveseqnn():
    for i,path in enumerate(result_paths):
        print('Starting path {}/{}'.format(i+1, len(result_paths)))
        twoOpt(path)
%time improveseqnn()

Starting path 1/10 exchanging 55-19,51-91 with 55-51,19-91 => save=2.6873147586143844 => cost: 386.386946361463 ... exchanging 67-39,75-0 with 67-75,39-0 => save=4.275672113835263 => cost: 368.50842741951476 Starting path 2/10 exchanging 17-97,78-0 with 17-78,97-0 => save=0.09782392307953103 => cost: 398.25554948502884 ... exchanging 88-78,26-0 with 88-26,78-0 => save=1.6805548383893267 => cost: 394.249851279636 Starting path 3/10 Starting path 4/10 exchanging 124-127,65-0 with 124-65,127-0 => save=0.5258127779938491 => cost: 390.465955878649 exchanging 14-65,33-0 with 14-33,65-0 => save=3.343148155406432 => cost: 387.1228077232426 Starting path 5/10 exchanging 0-18,141-69 with 0-141,18-69 => save=6.7520080457503155 => cost: 364.9242543330602 exchanging 50-121,146-0 with 50-146,121-0 => save=41.908681513869396 => cost: 323.0155728191908 Starting path 6/10 exchanging 47-45,101-138 with 47-101,45-138 => save=1.4275147235708587 => cost: 318.869322775654 Starting path 7/10 exchanging 3-35,83-77 with 3-83,35-77 => save=3.633104522812438 => cost: 361.20630895504394 ... exchanging 90-129,144-0 with 90-144,129-0 => save=5.489882024385494 => cost: 347.690672350374 Starting path 8/10 Starting path 9/10 Starting path 10/10 CPU times: user 3.25 s, sys: 37.7 ms, total: 3.29 s Wall time: 3.31 s

Plot improved paths

for path in result_paths:
    coords=get_coords(path)
    plt.plot(coords[0], coords[1], linewidth=2)
plt.scatter(data.x, data.y, c='r', s=30)
plt.scatter(data.iloc[0].x, data.iloc[0].y, marker='*', c='r', s=400)
plt.axis([-100,100,-100,100])
plt.grid(True)

Great! Several clean paths, but with many superpositions, let's improve them once more.

Scenario: center depot, several vans, KMeans clustering

Compute delivery clusters

from sklearn.cluster import KMeans
workpervan=20 #how many deliveries each van will start, ie granularity; a better approach would be capacity. model = KMeans(nclusters=1+size//workpervan, init='k-means++', randomstate=0) %time model.fit(data[['x', 'y']])

CPU times: user 75.7 ms, sys: 4.05 ms, total: 79.7 ms Wall time: 92.8 ms

KMeans(algorithm='auto', copyx=True, init='k-means++', maxiter=300, nclusters=8, ninit=10, njobs=1, precomputedistances='auto', random_state=0, tol=0.0001, verbose=0)

Plot clusters

plt.scatter(data.x, data.y, marker='o', c=model.labels_, s=200)
plt.scatter(model.clustercenters[:,0], model.clustercenters[:,1], marker='+', c='r', s=500)
plt.scatter(data.iloc[0].x, data.iloc[0].y, marker='*', c='r', s=400)
plt.axis([-100,100,-100,100])
plt.grid(True)

Let's apply a greedy NN algorithm sequentially, on each of the clusters

def greedynnclusters():
    result=[]
    for i in range(model.n_clusters):
        unserved=data[model.labels_ == i].index.tolist()
        if 0 not in unserved:
            unserved.append(0)
        resultpath = nearestneighbor(unserved)
        result.append(result_path)
        print('van {}: {}'.format(i+1, estimatecost(resultpath)))
    return result
%time resultpaths = greedynn_clusters()

van 1: 409.3080162039939 van 2: 266.4745473994821 van 3: 413.3606064604762 van 4: 446.80922078430706 van 5: 457.66477785115126 van 6: 455.80304195030294 van 7: 459.76894476196577 van 8: 376.5008949566861 CPU times: user 164 ms, sys: 10.1 ms, total: 174 ms Wall time: 174 ms

Plot result paths

plt.scatter(data.x, data.y, marker='o', c=model.labels_, s=200)
plt.scatter(model.clustercenters[:,0], model.clustercenters[:,1], marker='+', c='r', s=500)
for path in result_paths:
    coords=get_coords(path)
    plt.plot(coords[0], coords[1], 'k-', linewidth=2)
plt.scatter(data.iloc[0].x, data.iloc[0].y, marker='*', c='r', s=400)
plt.axis([-100,100,-100,100])
plt.grid(True)

Let's improve them with a 2-opt heuristic

def improveresultpaths():
    for i, path in enumerate(result_paths):
        print('Starting path {}/{}'.format(i+1, len(result_paths)))
        twoOpt(path)
%time improveresultpaths()

Starting path 1/8 exchanging 48-112,99-138 with 48-99,112-138 => save=11.5162692633929 => cost: 397.791746940601 exchanging 45-118,25-5 with 45-25,118-5 => save=13.984640981159586 => cost: 383.8071059594414 Starting path 2/8 exchanging 148-52,136-65 with 148-136,52-65 => save=1.6057804518926915 => cost: 264.8687669475894 exchanging 148-124,52-142 with 148-52,124-142 => save=6.904413311483367 => cost: 257.964353636106 Starting path 3/8 exchanging 63-24,108-111 with 63-108,24-111 => save=11.657922849361384 => cost: 401.7026836111148 ... exchanging 113-36,126-0 with 113-126,36-0 => save=7.772149456496223 => cost: 378.1902685773841 Starting path 4/8 exchanging 147-109,150-0 with 147-150,109-0 => save=5.601665728075147 => cost: 441.2075550562319 ... exchanging 62-109,143-150 with 62-143,109-150 => save=0.5147852919035643 => cost: 394.93412566866925 Starting path 5/8 exchanging 9-149,6-144 with 9-6,149-144 => save=0.0011015802705607314 => cost: 457.66367627088067 ... exchanging 105-129,144-90 with 105-144,129-90 => save=3.918490824545831 => cost: 407.2051777421843 Starting path 6/8 exchanging 30-69,34-116 with 30-34,69-116 => save=1.1688372549307786 => cost: 454.63420469537215 ... exchanging 22-116,79-38 with 22-79,116-38 => save=5.039989423215701 => cost: 380.39750649818404 Starting path 7/8 exchanging 0-55,1-107 with 0-1,55-107 => save=14.962441964752642 => cost: 444.8065027972131 ... exchanging 0-107,55-42 with 0-55,107-42 => save=6.044195398723886 => cost: 403.197665386089 Starting path 8/8 exchanging 0-71,26-59 with 0-26,71-59 => save=17.738760993960298 => cost: 358.76213396272584 ... exchanging 17-78,26-88 with 17-26,78-88 => save=0.17155006280330554 => cost: 319.53217682612376 CPU times: user 2.81 s, sys: 64.8 ms, total: 2.87 s Wall time: 3.03 s

Plot improved paths

plt.scatter(data.x, data.y, marker='o', c=model.labels_, s=200)
plt.scatter(model.clustercenters[:,0], model.clustercenters[:,1], marker='+', c='r', s=500)
for path in result_paths:
    coords=get_coords(path)
    plt.plot(coords[0], coords[1], 'k-', linewidth=2)
plt.scatter(data.iloc[0].x, data.iloc[0].y, marker='*', c='r', s=400)
plt.axis([-100,100,-100,100])
plt.grid(True)

Cool clean paths, with clearly defined cluster subsections!

Let's change the depot to the farthest point

Find farthest point and set it as new depot

maxid=data[0].idxmax()
print('farthest point={}; distance={}; coords={}'.format(maxid, data.iloc[maxid,2], tuple(data.iloc[maxid][['x','y']].tolist())))
data.drop(data.columns.difference(['x', 'y']), axis=1, inplace=True)
data.iloc[0]=data.iloc[maxid]
data.drop(maxid, inplace=True)
data.reset_index(drop=True, inplace=True)
size-=1
data

farthest point=84; distance=138.29710566715653; coords=(99.07656719022526, 96.48794364952255)