Travelling Salesman Problem - Best path to go through all points

Question

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?

Answer 1

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

Answer 2

Great job! Your code looks really good!

I am conducting research on a modified TSP problem so I decided to plug it into our algorithms. This is the route (I gave each point an id, 1 - 200):

[1, 32, 79, 72, 153, 44, 171, 144, 57, 42, 151, 191, 174, 178, 190, 13, 8, 61, 194, 86, 10, 38, 45, 82, 94, 74, 29, 95, 27, 15, 90, 28, 3, 101, 113, 166, 24, 63, 43, 104, 22, 125, 33, 70, 186, 187, 136, 75, 175, 159, 52, 56, 71, 128, 169, 48, 162, 198, 98, 87, 20, 73, 89, 55, 148, 107, 105, 149, 129, 155, 25, 36, 165, 30, 141, 200, 34, 176, 76, 47, 181, 182, 78, 154, 133, 195, 12, 139, 135, 126, 150, 88, 92, 130, 17, 121, 50, 62, 11, 192, 64, 66, 23, 120, 6, 160, 140, 81, 103, 111, 163, 37, 117, 143, 180, 167, 58, 196, 53, 132, 124, 77, 51, 19, 102, 138, 7, 168, 60, 97, 35, 164, 91, 46, 108, 170, 183, 137, 173, 189, 18, 41, 118, 122, 127, 39, 142, 21, 177, 172, 59, 83, 2, 185, 158, 31, 131, 14, 84, 184, 146, 134, 188, 112, 49, 69, 99, 147, 68, 123, 110, 85, 156, 5, 115, 4, 96, 161, 199, 116, 40, 193, 145, 54, 16, 93, 26, 80, 179, 100, 106, 119, 157, 9, 152, 197, 109, 114, 67, 65, 1]

4708.252935788349 km was the result in 2.599 seconds!

I solved this using an old C solution, Concorde, that just has a Python wrapper around it. Honestly, using your data just helped me find a bug in the algorithm I wrote myself! Thank you!

I can't exactly post Concorde's code sooo... and my wrapper code is currently private, sorry :)

Matplotlib graph of TSP Solution