I am trying to solve an exercise based on the travelling salesman problem. Basically I am provided with a list of points with their coordinates, like this:
[(523, 832), (676, 218), (731, 739), ..] (a total of 198)
and a starting one (832,500).
The goal is to find the optimal path to go through all points (starting from the (832,500) one). I have at disposal a service which will check my sequence and at the moment I have managed to get a decent solution (by applying the nearest neighbor algorithm), but it seems not to be the optimal solution. Here is my current code:
from math import sqrt
def calcDistance(start,finish):
dist = sqrt((start[0]-finish[0])*(start[0]-finish[0]) +
(start[1]-finish[1])*(start[1]-finish[1]))
return dist
def calcPath(path):
pathLen = 0
for i in range(1,len(path)):
pathLen+=calcDistance(path[i],path[i-1])
return pathLen
capital = (832, 500)
towns = [(523, 832), (676, 218), (731, 739), (803, 858),
(170, 542), (273, 743), (794, 345), (825, 569), (770, 306),
(168, 476), (198, 361), (798, 352), (604, 958), (700, 235),
(791, 661), (37, 424), (393, 815), (250, 719), (400, 183),
(468, 831), (604, 184), (168, 521), (691, 71), (304, 232),
(800, 642), (708, 241), (683, 223), (726, 257), (279, 252),
(559, 827), (832, 494), (584, 178), (254, 277), (309, 772),
(293, 240), (58, 658), (765, 300), (446, 828), (766, 699),
(407, 819), (818, 405), (626, 192), (828, 449), (758, 291),
(333, 788), (124, 219), (443, 172), (640, 801), (171, 452),
(242, 710), (496, 168), (217, 674), (785, 672), (369, 195),
(486, 168), (821, 416), (206, 654), (503, 832), (288, 756),
(789, 336), (170, 464), (636, 197), (168, 496), (832, 515),
(168, 509), (832, 523), (677, 781), (651, 796), (575, 176),
(478, 168), (831, 469), (391, 186), (735, 265), (529, 169),
(241, 292), (235, 700), (220, 321), (832, 481), (806, 629),
(176, 575), (751, 282), (511, 832), (581, 822), (708, 759),
(777, 317), (410, 180), (180, 411), (382, 189), (694, 230),
(327, 784), (177, 421), (797, 650), (742, 272), (719, 250),
(739, 731), (298, 764), (423, 177), (658, 792), (813, 611),
(667, 213), (257, 727), (178, 583), (616, 189), (342, 208),
(817, 600), (348, 205), (344, 793), (968, 541), (700, 766),
(181, 594), (633, 804), (656, 206), (831, 533), (722, 747),
(759, 708), (188, 615), (416, 822), (820, 590), (169, 529),
(172, 445), (424, 824), (687, 775), (229, 692), (597, 182),
(187, 388), (436, 826), (463, 170), (321, 220), (174, 434),
(567, 826), (224, 686), (210, 338), (608, 814), (190, 381),
(538, 170), (332, 938), (265, 735), (195, 367), (173, 562),
(270, 260), (462, 830), (192, 625), (824, 427), (781, 678),
(599, 817), (669, 786), (359, 199), (328, 216), (183, 401),
(815, 393), (827, 559), (830, 460), (215, 329), (311, 227),
(713, 755), (822, 581), (546, 829), (505, 168), (172, 554),
(748, 721), (421, 37), (184, 604), (317, 778), (286, 246),
(648, 202), (201, 645), (281, 750), (453, 171), (356, 800),
(827, 439), (491, 832), (375, 808), (807, 372), (521, 168),
(246, 286), (482, 832), (804, 365), (809, 622), (197, 637),
(232, 303), (227, 310), (362, 802), (592, 819), (533, 831),
(560, 173), (550, 171), (619, 810), (384, 811), (931, 313),
(811, 384), (168, 488), (773, 690), (781, 323), (204, 349),
(213, 667), (829, 547), (431, 175), (754, 714), (263, 267)]
#We start from the capital
start = capital
#path = ""
alreadyVisitedTowns = []
path = []
path2 = ""
while len(alreadyVisitedTowns) != 199:
#minDistanceForThisStep = {"town": towns[0], "end": calcDistance(start, towns[0]), "indices": 0}
bestTownIndex = 0
bestDistance = 999999
for i,town in enumerate(towns):
if (town not in alreadyVisitedTowns) and calcDistance(start, town)<bestDistance:
bestDistance = calcDistance(start, town)
bestTownIndex = i
path.append(bestTownIndex)
path2 += str(bestTownIndex) + " "
alreadyVisitedTowns.append(towns[bestTownIndex])
#path+= " " +str(minDistanceForThisStep["indices"])
start = towns[bestTownIndex]
print(path2)
print(" ".join(path2.split()[::-1]))
Result: 4749km
Test with networkx:
import numpy as np
import networkx as nx
from math import sqrt
#from scipy.spatial.distance import euclidean
def calcDistance(start,finish):
dist = sqrt((start[0]-finish[0])*(start[0]-finish[0]) +
(start[1]-finish[1])*(start[1]-finish[1]))
return dist
def calcPath(path):
pathLen = 0
for i in range(1,len(path)):
pathLen+=calcDistance(path[i],path[i-1])
return pathLen
# List of towns as (x, y) coordinates
towns = [(523, 832), (676, 218), (731, 739), (803, 858),
(170, 542), (273, 743), (794, 345), (825, 569), (770, 306),
(168, 476), (198, 361), (798, 352), (604, 958), (700, 235),
(791, 661), (37, 424), (393, 815), (250, 719), (400, 183),
(468, 831), (604, 184), (168, 521), (691, 71), (304, 232),
(800, 642), (708, 241), (683, 223), (726, 257), (279, 252),
(559, 827), (832, 494), (584, 178), (254, 277), (309, 772),
(293, 240), (58, 658), (765, 300), (446, 828), (766, 699),
(407, 819), (818, 405), (626, 192), (828, 449), (758, 291),
(333, 788), (124, 219), (443, 172), (640, 801), (171, 452),
(242, 710), (496, 168), (217, 674), (785, 672), (369, 195),
(486, 168), (821, 416), (206, 654), (503, 832), (288, 756),
(789, 336), (170, 464), (636, 197), (168, 496), (832, 515),
(168, 509), (832, 523), (677, 781), (651, 796), (575, 176),
(478, 168), (831, 469), (391, 186), (735, 265), (529, 169),
(241, 292), (235, 700), (220, 321), (832, 481), (806, 629),
(176, 575), (751, 282), (511, 832), (581, 822), (708, 759),
(777, 317), (410, 180), (180, 411), (382, 189), (694, 230),
(327, 784), (177, 421), (797, 650), (742, 272), (719, 250),
(739, 731), (298, 764), (423, 177), (658, 792), (813, 611),
(667, 213), (257, 727), (178, 583), (616, 189), (342, 208),
(817, 600), (348, 205), (344, 793), (968, 541), (700, 766),
(181, 594), (633, 804), (656, 206), (831, 533), (722, 747),
(759, 708), (188, 615), (416, 822), (820, 590), (169, 529),
(172, 445), (424, 824), (687, 775), (229, 692), (597, 182),
(187, 388), (436, 826), (463, 170), (321, 220), (174, 434),
(567, 826), (224, 686), (210, 338), (608, 814), (190, 381),
(538, 170), (332, 938), (265, 735), (195, 367), (173, 562),
(270, 260), (462, 830), (192, 625), (824, 427), (781, 678),
(599, 817), (669, 786), (359, 199), (328, 216), (183, 401),
(815, 393), (827, 559), (830, 460), (215, 329), (311, 227),
(713, 755), (822, 581), (546, 829), (505, 168), (172, 554),
(748, 721), (421, 37), (184, 604), (317, 778), (286, 246),
(648, 202), (201, 645), (281, 750), (453, 171), (356, 800),
(827, 439), (491, 832), (375, 808), (807, 372), (521, 168),
(246, 286), (482, 832), (804, 365), (809, 622), (197, 637),
(232, 303), (227, 310), (362, 802), (592, 819), (533, 831),
(560, 173), (550, 171), (619, 810), (384, 811), (931, 313),
(811, 384), (168, 488), (773, 690), (781, 323), (204, 349),
(213, 667), (829, 547), (431, 175), (754, 714), (263, 267)]
capital = (832, 500)
towns.append(capital)
num_towns = len(towns)
G = nx.complete_graph(num_towns)
for i in range(num_towns):
for j in range(i + 1, num_towns):
distance = calcDistance(towns[i], towns[j])
G[i][j]['weight'] = distance
G[j][i]['weight'] = distance
optimal_path_indices = list(nx.approximation.traveling_salesman_problem(G, cycle=True))
capital_index = optimal_path_indices.index(num_towns - 1)
optimal_path_indices = optimal_path_indices[capital_index:] + optimal_path_indices[:capital_index + 1]
optimal_path_indices = optimal_path_indices[:-1]
print("Optimal path indices:", optimal_path_indices)
optimal_path_indices = optimal_path_indices[::-1]
finalPath = ""
for el in optimal_path_indices:
finalPath += str(el) + " "
print(finalPath)
#print(calcPath(optimal_path_indices))
Result: 4725km (slightly better when inverting the path for some reason) Do you have any proposal on how I can improve my algorithm?
Euclidean TSP on 200 cities isn't too difficult to crack using the right tools. I used OR-Tools (as an interface to the mixed-integer program solver SCIP) and NetworkX. The code below takes a minute or two.
from collections import defaultdict
from itertools import combinations
from math import sqrt
from networkx import Graph, connected_components, dfs_preorder_nodes
from ortools.linear_solver import pywraplp
# Finds the optimal traveling salesman path from the first location, using
# integer programming.
def find_optimal_path(locations):
# Turn whatever sequence into a list.
locations = list(locations)
assert locations
# SCIP is our integer program solver of choice today.
solver = pywraplp.Solver.CreateSolver("SCIP")
# Our formulation will be to find the shortest tour minus an edge incident
# to the first location. Each edge has a Boolean variable that indicates
# whether it's in the tour.
variables = {edge: solver.BoolVar("") for edge in combinations(locations, 2)}
tour_length = sum(calcDistance(v, w) * x for ((v, w), x) in variables.items())
# One of the edges incident to the capital is free. We make variables
# indicating the identity of this edge.
exclusions = {(locations[0], w): solver.BoolVar("") for w in locations[1:]}
# One edge is excluded.
solver.Add(sum(exclusions.values()) == 1)
# It must be one of the tour edges.
for edge, x in exclusions.items():
solver.Add(x <= variables[edge])
# Minimize the length of the path.
solver.Minimize(
tour_length - sum(calcDistance(v, w) * x for ((v, w), x) in exclusions.items())
)
# Degree constraints: each node must be incident to exactly two path edges.
degrees = defaultdict(int)
for (v, w), x in variables.items():
degrees[v] += x
degrees[w] += x
for d in degrees.values():
solver.Add(d == 2)
# Let's see what's going on.
solver.EnableOutput()
while True:
status = solver.Solve()
# "Wait a second," you say. "There aren't enough constraints!" You're
# absolutely right. We're probably going to get back a collection of
# multiple disjoint paths. If we do, add *sub-path elimination
# constraints* and try again.
assert status == solver.OPTIMAL
graph = Graph()
graph.add_edges_from(
edge for (edge, x) in variables.items() if x.solution_value() > 0.5
)
graph.remove_edges_from(
edge for (edge, x) in exclusions.items() if x.solution_value() > 0.5
)
components = list(connected_components(graph))
if len(components) == 1:
return list(dfs_preorder_nodes(graph, locations[0]))
for component in components:
solver.Add(
sum(
x
for ((v, w), x) in variables.items()
if (v in component) != (w in component)
)
>= 2
)
# Your code below.
def calcDistance(start, finish):
dist = sqrt(
(start[0] - finish[0]) * (start[0] - finish[0])
+ (start[1] - finish[1]) * (start[1] - finish[1])
)
return dist
def calcPath(path):
pathLen = 0
for i in range(1, len(path)):
pathLen += calcDistance(path[i], path[i - 1])
return pathLen
capital = (832, 500)
towns = [
(523, 832),
(676, 218),
(731, 739),
(803, 858),
(170, 542),
(273, 743),
(794, 345),
(825, 569),
(770, 306),
(168, 476),
(198, 361),
(798, 352),
(604, 958),
(700, 235),
(791, 661),
(37, 424),
(393, 815),
(250, 719),
(400, 183),
(468, 831),
(604, 184),
(168, 521),
(691, 71),
(304, 232),
(800, 642),
(708, 241),
(683, 223),
(726, 257),
(279, 252),
(559, 827),
(832, 494),
(584, 178),
(254, 277),
(309, 772),
(293, 240),
(58, 658),
(765, 300),
(446, 828),
(766, 699),
(407, 819),
(818, 405),
(626, 192),
(828, 449),
(758, 291),
(333, 788),
(124, 219),
(443, 172),
(640, 801),
(171, 452),
(242, 710),
(496, 168),
(217, 674),
(785, 672),
(369, 195),
(486, 168),
(821, 416),
(206, 654),
(503, 832),
(288, 756),
(789, 336),
(170, 464),
(636, 197),
(168, 496),
(832, 515),
(168, 509),
(832, 523),
(677, 781),
(651, 796),
(575, 176),
(478, 168),
(831, 469),
(391, 186),
(735, 265),
(529, 169),
(241, 292),
(235, 700),
(220, 321),
(832, 481),
(806, 629),
(176, 575),
(751, 282),
(511, 832),
(581, 822),
(708, 759),
(777, 317),
(410, 180),
(180, 411),
(382, 189),
(694, 230),
(327, 784),
(177, 421),
(797, 650),
(742, 272),
(719, 250),
(739, 731),
(298, 764),
(423, 177),
(658, 792),
(813, 611),
(667, 213),
(257, 727),
(178, 583),
(616, 189),
(342, 208),
(817, 600),
(348, 205),
(344, 793),
(968, 541),
(700, 766),
(181, 594),
(633, 804),
(656, 206),
(831, 533),
(722, 747),
(759, 708),
(188, 615),
(416, 822),
(820, 590),
(169, 529),
(172, 445),
(424, 824),
(687, 775),
(229, 692),
(597, 182),
(187, 388),
(436, 826),
(463, 170),
(321, 220),
(174, 434),
(567, 826),
(224, 686),
(210, 338),
(608, 814),
(190, 381),
(538, 170),
(332, 938),
(265, 735),
(195, 367),
(173, 562),
(270, 260),
(462, 830),
(192, 625),
(824, 427),
(781, 678),
(599, 817),
(669, 786),
(359, 199),
(328, 216),
(183, 401),
(815, 393),
(827, 559),
(830, 460),
(215, 329),
(311, 227),
(713, 755),
(822, 581),
(546, 829),
(505, 168),
(172, 554),
(748, 721),
(421, 37),
(184, 604),
(317, 778),
(286, 246),
(648, 202),
(201, 645),
(281, 750),
(453, 171),
(356, 800),
(827, 439),
(491, 832),
(375, 808),
(807, 372),
(521, 168),
(246, 286),
(482, 832),
(804, 365),
(809, 622),
(197, 637),
(232, 303),
(227, 310),
(362, 802),
(592, 819),
(533, 831),
(560, 173),
(550, 171),
(619, 810),
(384, 811),
(931, 313),
(811, 384),
(168, 488),
(773, 690),
(781, 323),
(204, 349),
(213, 667),
(829, 547),
(431, 175),
(754, 714),
(263, 267),
]
print(calcPath(find_optimal_path([capital] + towns)))