5

Using Python and PuLP library, how can we create the linear programming model to solve the Traveling Salesman Problem (TSP)?

From Wikipedia, the objective function and constraints are

enter image description here

Problem: Here is my partial attempt where I am stuck.

  1. I did not include the final constraint in the code because I dont know how to define it. I believe this constraint with u variables are for preventing sub-cycles in the solution

  2. Also, solving for the current model gives decision variables such as x0_0 and x1_1 being equal to 1.0 which is definitely wrong... I can't figure out why this is so even though I had

        if i == j:
        upperBound = 0

Python Code

import pulp

def get_dist(tsp):
    with open(tsp, 'rb') as tspfile:
        r = csv.reader(tspfile, delimiter='\t')
        d = [row for row in r]

    d = d[1:] # skip header row
    locs = set([r[0] for r in d]) | set([r[1] for r in d])
    loc_map = {l:i for i, l in enumerate(locs)}
    idx_map = {i:l for i, l in enumerate(locs)}
    dist = [(loc_map[r[0]], loc_map[r[1]], r[2]) for r in d]
    return dist, idx_map


def dist_from_coords(dist, n):
    points = []
    for i in range(n):
        points.append([0] * n)
    for i, j, v in dist:
        points[i][j] = points[j][i] = float(v)
    return points


def find_tour():
    tsp_file = '/Users/test/' + 'my-waypoints-dist-dur.tsv'
    coords, idx_map = get_dist(tsp_file)
    n = len(idx_map)
    dist = dist_from_coords(coords, n)


    # Define the problem
    m = pulp.LpProblem('TSP', pulp.LpMinimize)


    # Create variables
    # x[i,j] is 1 if edge i->j is on the optimal tour, and 0 otherwise
    # Also forbid loops
    x = {}
    for i in range(n):
        for j in range(n):
            lowerBound = 0
            upperBound = 1

            # Forbid loops
            if i == j:
                upperBound = 0
                # print i,i

            x[i,j] = pulp.LpVariable('x' + str(i) + '_' + str(j), lowerBound, upperBound, pulp.LpBinary)
            # x[j,i] = x[i,j]


    # Define the objective function to minimize
    m += pulp.lpSum([dist[i][j] * x[i,j] for i in range(n) for j in range(n)])


    # Add degree-2 constraint
    for i in range(n):
        m += pulp.lpSum([x[i,j] for j in range(n)]) == 2


    # Solve and display results
    status = m.solve()
    print pulp.LpStatus[status]
    for i in range(n):
        for j in range(n):
            if pulp.value(x[i,j]) >0:
                print str(i) + '_' + str(j) + ': ' + str( pulp.value(x[i,j]) )

find_tour()

my-waypoints-dist-dur.tsv

The data file can be found here.

Result

0_0: 1.0
0_5: 1.0
1_1: 1.0
1_15: 1.0
2_2: 1.0
2_39: 1.0
3_3: 1.0
3_26: 1.0
4_4: 1.0
4_42: 1.0
5_5: 1.0
5_33: 1.0
6_6: 1.0
6_31: 1.0
7_7: 1.0
7_38: 1.0
8_8: 1.0
8_24: 1.0
9_9: 1.0
9_26: 1.0
10_4: 1.0
10_10: 1.0
11_11: 1.0
11_12: 1.0
12_11: 1.0
12_12: 1.0
13_13: 1.0
13_17: 1.0
14_14: 1.0
14_18: 1.0
15_1: 1.0
15_15: 1.0
16_3: 1.0
16_16: 1.0
17_13: 1.0
17_17: 1.0
18_14: 1.0
18_18: 1.0
19_19: 1.0
19_20: 1.0
20_4: 1.0
20_20: 1.0
21_21: 1.0
21_25: 1.0
22_22: 1.0
22_27: 1.0
23_21: 1.0
23_23: 1.0
24_8: 1.0
24_24: 1.0
25_21: 1.0
25_25: 1.0
26_26: 1.0
26_43: 1.0
27_27: 1.0
27_38: 1.0
28_28: 1.0
28_47: 1.0
29_29: 1.0
29_31: 1.0
30_30: 1.0
30_34: 1.0
31_29: 1.0
31_31: 1.0
32_25: 1.0
32_32: 1.0
33_28: 1.0
33_33: 1.0
34_30: 1.0
34_34: 1.0
35_35: 1.0
35_42: 1.0
36_36: 1.0
36_47: 1.0
37_36: 1.0
37_37: 1.0
38_27: 1.0
38_38: 1.0
39_39: 1.0
39_44: 1.0
40_40: 1.0
40_43: 1.0
41_41: 1.0
41_45: 1.0
42_4: 1.0
42_42: 1.0
43_26: 1.0
43_43: 1.0
44_39: 1.0
44_44: 1.0
45_15: 1.0
45_45: 1.0
46_40: 1.0
46_46: 1.0
47_28: 1.0
47_47: 1.0



...

Updated Code

import csv
import pulp


def get_dist(tsp):
    with open(tsp, 'rb') as tspfile:
        r = csv.reader(tspfile, delimiter='\t')
        d = [row for row in r]

    d = d[1:] # skip header row
    locs = set([r[0] for r in d]) | set([r[1] for r in d])
    loc_map = {l:i for i, l in enumerate(locs)}
    idx_map = {i:l for i, l in enumerate(locs)}
    dist = [(loc_map[r[0]], loc_map[r[1]], r[2]) for r in d]
    return dist, idx_map


def dist_from_coords(dist, n):
    points = []
    for i in range(n):
        points.append([0] * n)
    for i, j, v in dist:
        points[i][j] = points[j][i] = float(v)
    return points


def find_tour():
    tsp_file = '/Users/test/' + 'my-waypoints-dist-dur.tsv'
    coords, idx_map = get_dist(tsp_file)
    n = len(idx_map)
    dist = dist_from_coords(coords, n)


    # Define the problem
    m = pulp.LpProblem('TSP', pulp.LpMinimize)


    # Create variables
    # x[i,j] is 1 if edge i->j is on the optimal tour, and 0 otherwise
    # Also forbid loops
    x = {}
    for i in range(n+1):
        for j in range(n+1):
            lowerBound = 0
            upperBound = 1

            # Forbid loops
            if i == j:
                upperBound = 0
                # print i,i

            # Create the decision variable and First constraint
            x[i,j] = pulp.LpVariable('x' + str(i) + '_' + str(j), lowerBound, upperBound, pulp.LpBinary)
            # x[j,i] = x[i,j]


    # Define the objective function to minimize
    m += pulp.lpSum([dist[i][j] * x[i,j] for i in range(1,n+1) for j in range(1,n+1)])


    # Add degree-2 constraint (3rd and 4th)
    for i in range(1,n+1):
        m += pulp.lpSum([x[i,j] for j in range(1,n+1)]) == 2


    # Add the last (5th) constraint (prevents subtours)
    u = []
    for i in range(1, n+1):
        u.append(pulp.LpVariable('u_' + str(i), cat='Integer'))
    for i in range(1, n-1):
        for j in range(i+1, n+1):
            m += pulp.lpSum([ u[i] - u[j] + n*x[i,j]]) <= n-1


    # status = m.solve()
    # print pulp.LpStatus[status]
    # for i in range(n):
    #   for j in range(n):
    #       if pulp.value(x[i,j]) >0:
    #           print str(i) + '_' + str(j) + ': ' + str( pulp.value(x[i,j]) )

find_tour()
Ami Tavory
  • 66,807
  • 9
  • 114
  • 153
Nyxynyx
  • 52,075
  • 130
  • 401
  • 707
  • What problem are you having with your code? – kindall Oct 06 '16 at 06:00
  • @kindall I did not include the final constraint because I dont know how to define it. Also, solving for the current model gives decision variables such as `x1_1` being equal to `1.0` which is definitely wrong... I can't figure out why this is so. – Nyxynyx Oct 06 '16 at 06:10
  • @kindall Updated the question to clarify the problem, and included the incorrect results with the current code. – Nyxynyx Oct 06 '16 at 06:15
  • You won't be able to solve larger problems with this approach because the formulation will grow very big due to the sub tour elimination constraints. You might consider another approach that dynamically separates invalid solutions. There is an example in Python here: https://github.com/SCIP-Interfaces/PySCIPOpt/blob/master/tests/test_tsp.py – mattmilten Oct 07 '16 at 05:51
  • @mattmilten Thanks for the suggestion, very useful because I am new to TSP and LP and do not know what are the more efficient approaches, especially the ones used in real life. – Nyxynyx Oct 07 '16 at 19:04
  • @mattmilten Is there a name for the approach that you linked to? – Nyxynyx Oct 07 '16 at 19:13
  • @Nyxynyx Look for this paper: Solution of a Large-Scale Traveling-Salesman Problem by Dantzig, Fulkerson and Johnson: www.cs.uleth.ca/~benkoczi/OR/read/tsp-dantzig-fulkerson-johnson-54.pdf – mattmilten Oct 07 '16 at 19:52
  • This will probably be helpful to you: https://nbviewer.jupyter.org/github/cochoa0x1/intro-mathematical-programming/blob/master/05-routes-and-schedules/traveling_salesman.ipynb – Akavall Apr 09 '19 at 22:16

1 Answers1

1

The last constraint is not a single constraint. You should add one constraint for each pair of indices i, j that satisfy that condition:

for i in range(n-1):
    for j in range(i+1, n):
        m += pulp.lpSum([ u[i] - u[j] + n*x[i,j]]) <= n-1

However you first have to declare the u_i variables as integers, passing the cat='Integer' argument to LpVariable:

u = []
for i in range(n):
    u.append(pulp.LpVariable('u_' + str(i), cat='Integer'))
Bakuriu
  • 85,459
  • 18
  • 168
  • 202
  • Thanks, I'm looking into your suggestion. Is the current code looking ok so far? Is the lack of the last constraint the reason `x0_0`, `x1_1` etc are being solved to have values of `1.0` instead of `0.0`? – Nyxynyx Oct 06 '16 at 06:19
  • @Nyxynyx The last constraint is absolutely fundamental for the problem. Without it the problem is not even an *integer* linear problem, but just a linear problem and as such it cannot represent TSP in any meaningful way (remember that TSP cannot be solved in polynomial time and cannot even be approximated in "a good way" in polynomial time, hence every non-integer linear programming problem (which can always be solved in polynomial time) cannot represent it or any kind of good approximation) – Bakuriu Oct 06 '16 at 06:22
  • As I am trying to assimilate the code for the last constraint, I added them to the current code and tried to solve for it, but the results still look strange, getting values of `1.0` for variables `x[i,j]` where `i == j`. – Nyxynyx Oct 06 '16 at 06:28
  • @Nyxynyx I believe there could be an error in your code. In the original model, if I'm not mistaken, the node with index `0` should be node that "closes the cycle". You can see that in the objective function you have the outer sum from `i=0` to `n`. However the other constraints start at `i=1`. In your code it seems like you thought that the model was simply `1`-indexed and thus changed the constraints to `0`-indexed... – Bakuriu Oct 06 '16 at 06:32
  • @Nyxynyx You should try to define your `x` variables using `for i in range(n+1): for j in range(n+1)`, and modify all the constraints to start at `1`, so even this last one should be `for i in range(1, n): for j in range(i+1, n+1): ...`. Also there is no `u_0` variable, so you'll have to adjust for this when defining the `u` vector... just add a fake entry at the beginning: `u = [None]; for i in range(1, n+1): u.append(...)` in this way `u[i]` refers to city with index `i` and you don't have to add `-1` in the final constraint. – Bakuriu Oct 06 '16 at 06:36
  • Thanks for pointing out the indexing problem. If I were to "Label the cities with the numbers 1, …, n ", but create the `x` variables with `for i in range(n+1): for j in range(n+1)`?, which city does `x[0,0]` corresponds to? I found this out because `dist[i][j]` gives an index out of range error when we are generating the objective function – Nyxynyx Oct 06 '16 at 06:53
  • @Nyxynyx If you look carefully at the objective function you'll see that the `j` index ranges from `1` to `n` and has also the constraint `j != i`, so you should *never* end up computing `x[0,0]` in the first place. – Bakuriu Oct 06 '16 at 08:26
  • @Nyxynyx I believe you should read [this article](https://www.unc.edu/~pataki/papers/teachtsp.pdf) about the formulations of the TSP problem. The first model is the same of wikipedia (although it's written in slightly a different way), but it has a lot more explanation about the model etc. Right now I don't really have time to look further into this. – Bakuriu Oct 06 '16 at 08:33
  • I think I'm almost there! Just cant figure out how to define city with index `0`, which is the city node that closes the loop of the tour... any suggestions? – Nyxynyx Oct 06 '16 at 22:29
  • On the topic of the city indexed `0`, if it is the city that closes the loop of `n` cities in the tour, must this city be manually chosen (possibly randomly too) and added to the list of `n` cities indexed `1` to `n`, such that our list of cities to be toured now contains `n+1` cities, where only `n` are unique ? – Nyxynyx Oct 06 '16 at 22:33