Transport Analytics Training Series - Last Revision: October 2022
Transport Analytics Training Series - Last Revision: October 2022
It's Christmas Eve! And Father Christmas is running out of time! And he has so many ore homes to visit!
What's the best way to find the shortest tour visiting each home only once?
This notebook compares the solutions for travelling salesperson problem derived by using the following algorithms:
The comparison is shown here.
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
from matplotlib import collections
import matplotlib.cm as cm
import seaborn as sns
from sklearn.datasets import make_blobs
from scipy.spatial import distance
import pandas as pd
from pulp import *
from deap import base, creator, tools
import copy
from time import time
%matplotlib inline
We will first create a map of the 15 homes that all the five alogirthms can use as an input.
(The colors don't mean anything -- the multi-color plot is to make it more festive! )
num_homes = 15
sns.set_style("darkgrid", {"axes.facecolor": "honeydew"})
plot_size = 16
plot_width = 20
plot_height = 10
params = {'legend.fontsize': 'large',
'figure.figsize': (plot_width,plot_height),
'axes.labelsize': plot_size,
'axes.titlesize': plot_size,
'xtick.labelsize': plot_size*0.75,
'ytick.labelsize': plot_size*0.75,
'axes.titlepad': 25}
plt.rcParams.update(params)
def create_map(homes_coord, homes_coord_dict):
colors = cm.rainbow(np.linspace(0, 1, len(homes_coord[:, 0])))
plt.scatter(homes_coord[:, 0],
homes_coord[:, 1],
s=plot_size*3,
color = colors)
for i in range(num_homes):
plt.annotate(i, (homes_coord[i, 0]+1, homes_coord[i, 1]+0.5), size=14)
homes_ids = [f'Node {i}' for i in range(num_homes)]
center_box = (100, 200)
homes_coord, _ = make_blobs(n_samples=num_homes,
centers=4, # this is to make the home more spatially distributed
cluster_std=20,
center_box=center_box,
random_state = 10)
homes_coord_dict = {name: coord for name,coord in zip(homes_ids, homes_coord)}
create_map(homes_coord, homes_coord_dict)
plt.show()
# create the distance matrix of the homes
dist_matrix = distance.cdist(homes_coord, homes_coord, 'euclidean')
dist_matrix = np.where(dist_matrix==0, np.Inf, dist_matrix)
# Generic plot solution function
def plot_solution(_N, _path, _homes_coord, _homes_coord_dict):
create_map(_homes_coord, _homes_coord_dict)
lines = []
for p in range(len(_path) - 1):
i = _path[p]
j = _path[p+1]
colour = 'black'
plt.arrow(_homes_coord_dict[i][0], _homes_coord_dict[i][1],_homes_coord_dict[j][0] - _homes_coord_dict[i][0], _homes_coord_dict[j][1] - _homes_coord_dict[i][1], color=colour)
plt.show()
start_node = homes_ids[0] # start from a random home - in our case the first home
t0 = time()
name_to_index = {i:idx for idx,i in enumerate(homes_ids)}
visited = [False if start_node is not _ else True for _ in homes_ids]
route_NN = [start_node]
curr = start_node
distance_NN = 0
dist_matrix_copy = copy.deepcopy(dist_matrix)
while False in visited:
row = dist_matrix_copy[name_to_index[curr]]
for idx in range(len(row)):
if visited[idx]:
row[idx] = np.Inf
closest = np.argmin(row)
visited[closest] = True
curr = homes_ids[closest]
route_NN.append(curr)
distance_NN = distance_NN + dist_matrix[[int(s) for s in curr.split() if s.isdigit()][0]][[int(s) for s in route_NN[-2].split() if s.isdigit()][0]]
route_NN.append(start_node)
distance_NN = distance_NN + dist_matrix[[int(s) for s in start_node.split() if s.isdigit()][0]][[int(s) for s in route_NN[-2].split() if s.isdigit()][0]]
time_NN = time() - t0
# for full details see notebook n10_3 - Travelling Salesman Problem - Nearest Neighbour Heuristic
This is the same with nearest neighbour algorithm with N steps but repeating the method for all nodes as starting node.
t0 = time()
distance_N2 = {}
route_N2 = {}
for i in range(len(homes_ids)):
start_node = homes_ids[i]
name_to_index = {i:idx for idx,i in enumerate(homes_ids)}
visited = [False if start_node is not _ else True for _ in homes_ids]
route_N2[i] = [start_node]
curr = start_node
distance_N2[i] = 0
dist_matrix_copy = copy.deepcopy(dist_matrix)
while False in visited:
row = dist_matrix_copy[name_to_index[curr]]
for idx in range(len(row)):
if visited[idx]:
row[idx] = np.Inf
closest = np.argmin(row)
visited[closest] = True
curr = homes_ids[closest]
distance_N2[i] = distance_N2[i] + dist_matrix[[int(s) for s in curr.split() if s.isdigit()][0]][[int(s) for s in route_N2[i][-1].split() if s.isdigit()][0]]
route_N2[i].append(curr)
distance_N2[i] = distance_N2[i] + dist_matrix[[int(s) for s in route_N2[i][0].split() if s.isdigit()][0]][[int(s) for s in route_N2[i][-1].split() if s.isdigit()][0]]
route_N2[i].append(route_N2[i][0])
time_N2 = time() - t0
shortest_dist_N2 = 9999999999999999999999
for k,v in distance_N2.items():
if v < shortest_dist_N2:
shortest_dist_N2 = v
shortest_route_N2 = route_N2[k]
Take the shortestroute calculated from nearest neighbour with $|N|^{2}$ steps
dist_matrix_copy = copy.deepcopy(dist_matrix)
shortest_route_N2_copy = copy.deepcopy(shortest_route_N2)
shortest_route_N2_copy.pop()
'Node 13'
def route_eval(_dist_matrix_copy, _route):
dist = 0
for p in range(len(_route) - 1):
_i = _route[p]
_j = _route[p+1]
dist += _dist_matrix_copy[[int(s) for s in _i.split() if s.isdigit()][0]][[int(s) for s in _j.split() if s.isdigit()][0]]
dist += dist_matrix_copy[[int(s) for s in _route[-1].split() if s.isdigit()][0]][[int(s) for s in _route[0].split() if s.isdigit()][0]]
return dist
def two_opt(_dist_matrix_copy, _route):
stay_in_loop = True
best_route = _route
route_length = len(best_route)
best_fit = route_eval(_dist_matrix_copy, best_route)
while stay_in_loop:
stay_in_loop = False
for i in range(route_length):
for k in range(i, route_length):
new_route = best_route[0:i] + list(reversed(best_route[i:k])) + best_route[k:]
new_route_cost = route_eval(_dist_matrix_copy, new_route)
if new_route_cost < best_fit:
best_fit = new_route_cost
best_route = new_route
stay_in_loop = True
break
if stay_in_loop:
break
best_route.append(best_route[0])
return best_route, best_fit
t0 = time()
two_opt_route, two_opt_dist = two_opt(dist_matrix_copy, shortest_route_N2_copy)
two_opt_time = time() - t0
dist_matrix_copy = copy.deepcopy(dist_matrix)
prob = LpProblem('prob', LpMinimize)
N = [('Node '+str(i)) for i in range(num_homes)] # the route will start from Node 1
c = {n_f:{n_t: round(dist_matrix_copy[i][j],2) for j,n_t in enumerate(N)} for i,n_f in enumerate(N)}
x = LpVariable.dicts('x',(N,N), 0,1,LpBinary)
prob += lpSum([lpSum([x[i][j]*c[i][j] for j in N if i!=j]) for i in N])
for i in N:
prob += lpSum([x[i][j] for j in N if i!=j]) == 1
for j in N:
prob += lpSum([x[i][j] for i in N if i!=j]) == 1
# we need to keep track of the order in the tour to eliminate the possibility of subtours
u = LpVariable.dicts('u', N, 0, num_homes-1, LpInteger)
U=num_homes
for i in N:
for j in N:
if i != j and (i != 'Node 0' and j!= 'Node 0'):
prob += u[i] - u[j] <= U*(1-x[i][j]) - 1
t0 = time()
%time status = prob.solve()
print(f'\nSTATUS\n{LpStatus[status]}\n')
time_MTZ = time() - t0
CPU times: user 6.28 ms, sys: 9.54 ms, total: 15.8 ms Wall time: 1.29 s STATUS Optimal
route_MTZ = []
distance_MTZ = 0
N_unvisited = N.copy()
current = 'Node 0'
route_MTZ.append(N_unvisited.pop(N_unvisited.index(current)))
while len(N_unvisited) > 0:
for k in N_unvisited:
if x[current][k].varValue ==1:
route_MTZ.append(N_unvisited.pop(N_unvisited.index(k)))
distance_MTZ = distance_MTZ + dist_matrix_copy[[int(s) for s in route_MTZ[-1].split() if s.isdigit()][0]][[int(s) for s in route_MTZ[-2].split() if s.isdigit()][0]]
current=k
break
distance_MTZ = distance_MTZ + dist_matrix_copy[[int(s) for s in route_MTZ[-1].split() if s.isdigit()][0]][[int(s) for s in route_MTZ[0].split() if s.isdigit()][0]]
route_MTZ.append('Node 0')
dist_matrix_copy = copy.deepcopy(dist_matrix)
np.random.seed(3) # we are setting the same random seed to be able to reproduce the answers.
def chromo_create(_homes_ids):
chromo = copy.deepcopy(_homes_ids)
np.random.shuffle(chromo)
return chromo
def chromo_eval(_dist_matrix_copy, _chromo):
dist = 0
for p in range(len(_chromo) - 1):
_i = _chromo[p]
_j = _chromo[p+1]
dist += _dist_matrix_copy[[int(s) for s in _i.split() if s.isdigit()][0]][[int(s) for s in _j.split() if s.isdigit()][0]]
dist += dist_matrix_copy[[int(s) for s in _chromo[-1].split() if s.isdigit()][0]][[int(s) for s in _chromo[0].split() if s.isdigit()][0]]
return dist,
tb = base.Toolbox()
creator.create('Fitness_Func', base.Fitness, weights=(-1.0,))
creator.create('Individual', list, fitness=creator.Fitness_Func)
num_population = 200
num_generations = 2000
prob_crossover = .4
prob_mutation = .6
tb.register('indexes', chromo_create, homes_ids)
tb.register('individual', tools.initIterate, creator.Individual, tb.indexes)
tb.register('population', tools.initRepeat, list, tb.individual)
tb.register('evaluate', chromo_eval, dist_matrix_copy)
tb.register('select', tools.selTournament)
tb.register('mate', tools.cxPartialyMatched)
tb.register('mutate', tools.mutShuffleIndexes)
population = tb.population(n=num_population)
fitness_set = list(tb.map(tb.evaluate, population))
for ind, fit in zip(population, fitness_set):
ind.fitness.values = fit
t0 = time()
best_fit_list = []
best_sol_list = []
best_fit = np.Inf
for gen in range(0, num_generations):
offspring = tb.select(population, len(population), tournsize=3)
offspring = list(map(tb.clone, offspring))
for child1, child2 in zip(offspring[0::2], offspring[1::2]):
child1_no = [[int(s) for s in i.split() if s.isdigit()][0] for i in child1]
child2_no = [[int(s) for s in i.split() if s.isdigit()][0] for i in child2]
if np.random.random() < prob_crossover:
tb.mate(child1_no, child2_no)
del child1.fitness.values
del child2.fitness.values
for chromo in offspring:
if np.random.random() < prob_mutation:
tb.mutate(chromo, indpb=0.01)
del chromo.fitness.values
invalid_ind = [ind for ind in offspring if not ind.fitness.valid]
fitness_set = map(tb.evaluate, invalid_ind)
for ind, fit in zip(invalid_ind, fitness_set):
ind.fitness.values = fit
population[:] = offspring
curr_best_sol = tools.selBest(population, 1)[0]
curr_best_fit = curr_best_sol.fitness.values[0]
if curr_best_fit < best_fit:
best_sol = curr_best_sol
best_fit = curr_best_fit
best_fit_list.append(best_fit)
best_sol_list.append(best_sol)
time_GA = time() - t0
plt.plot(best_fit_list)
plt.show()
distance_GA = best_fit_list[-1]
route_GA = []
final_sol = best_sol + best_sol[0:1]
route_GA = final_sol
# Nearest neighbor
plot_solution(homes_ids, route_NN, homes_coord, homes_coord_dict)
print(distance_NN)
print(route_NN)
print(time_NN)
335.1060231126455 ['Node 0', 'Node 13', 'Node 10', 'Node 6', 'Node 5', 'Node 14', 'Node 9', 'Node 12', 'Node 11', 'Node 8', 'Node 1', 'Node 7', 'Node 4', 'Node 2', 'Node 3', 'Node 0'] 0.0008549690246582031
# Nearest neighbor with N2
plot_solution(homes_ids, shortest_route_N2, homes_coord, homes_coord_dict)
print(shortest_route_N2)
print(shortest_dist_N2)
print(time_N2)
['Node 13', 'Node 0', 'Node 12', 'Node 9', 'Node 14', 'Node 11', 'Node 8', 'Node 1', 'Node 7', 'Node 4', 'Node 2', 'Node 3', 'Node 6', 'Node 5', 'Node 10', 'Node 13'] 297.3262161693636 0.0034859180450439453
# Nearest neighbor with 2-opt
plot_solution(homes_ids, two_opt_route, homes_coord, homes_coord_dict)
print(two_opt_route)
print(two_opt_dist)
print(two_opt_time)