Travelling Salesman Problem - Comparison of solutions¶

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

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)

# 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

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]

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

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

distance_GA = best_fit_list[-1]
route_GA = []
final_sol = best_sol + best_sol[0:1]
route_GA = final_sol

def subplot_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)