pythonalgorithmnearest-neighborpulptraveling-salesman

TSP with only one key to identify a road

I need to solve Traveling Salesman Problem in Python using PuLP.

import pulp
import numpy as np
import matplotlib.pyplot as plt

n = 50
np.random.seed(42)
x = 1.5*np.random.rand(n)
y = np.random.rand(n)

Roads are declared like this:

roads = pulp.LpVariable.dicts("Road", (range(n), range(n)), 0, 1, pulp.LpInteger)

Distance calculation:

xi, xj = np.meshgrid(x, x)
yi, yj = np.meshgrid(y, y)
dist = np.hypot(xi - xj, yi - yj)

Constraints:

for i in range(n):
    prob += pulp.lpSum([roads[i][j] for j in range(n)]) == 1
for j in range(n):
    prob += pulp.lpSum([roads[i][j] for i in range(n)]) == 1

for i in range(n):
    prob += roads[i][i] == 0

prob += pulp.lpSum([dist[i][j]*roads[i][j] for i in range(n) for j in range(n)]), "Objective Function"

def draw(x, y, n, path):
    plt.plot(x, y, '*', markerfacecolor = 'red', markeredgecolor = 'red')
    plt.axis('equal')
    for i in range(n):
        plt.plot((x[i], x[path[i]]), (y[i], y[path[i]]), '--b')

prob.solve()

def findPath(roads):
    path = [0]*n
    for i in range(n):
        for j in range(n):
            if pulp.value(roads[i][j]) == 1:
                path[i] = j
    return path

path = findPath(roads)
draw(x, y, n, path)

def findTours(path):
    ipath = [True]*len(path)
    tours = []
    while True in ipath:
        i0 = ipath.index(True)
        ipath[i0] = False
        tour = [i0]
        i = path[i0]
        while i != i0:
            tour.append(i)
            ipath[i] = False
            i = path[i]
        tours.append(tour)
    return tours

itCount = 0
while True:
    itCount += 1
    print('#', itCount, end = '')
    result = prob.solve()
    if result != 1:
        break
    path = findPath(roads) 
    tours = findTours(path)
    print(' => ', len(tours))
    if len(tours) == 1:
        break
    for tour in tours:
        prob += pulp.lpSum([roads[i][j] for i in tour for j in tour]) <= len(tour) - 1

draw(x, y, n, path)
print("Result = ", pulp.value(prob.objective))

This code works correctly, I need to use only one variable to identify a road, i.e.:

roads = pulp.LpVariable.dicts("Road", (range(n), range(n)), 0, 1, pulp.LpInteger)

I rewrote constraints like this:

for i in range(n):
    prob += pulp.lpSum([roads[(i,j)] for j in range(i + 1, n)]) + pulp.lpSum([roads[(j,i)] for j in range(0, i)]) == 2

prob += pulp.lpSum([dist[i][j] * roads[i,j] for i in range(n) for j in range(i + 1, n)]), "Objective Function"

And path calculation like this:

def findPath(roads):
    path = [[] for _ in range(n)]
    for i in range(n):
        for j in range(i + 1, n):
            if pulp.value(roads[i, j]) == 1:
                path[i].append(j)
    return path

Now I have a problem with function findTours which combines connected components into one graph using the Nearest neighbour algorithm - obviously it does not work with 2D arrays and I don't see an efficient way to re-implement it to support 2D arrays.

For experiments I reduced n to 10 and rewrote the function like this, but it struggles to combine everything to one graph:

def findTours(path):
    ipath = [True]*len(path)
    tours = []
    while True in ipath:
        i0 = ipath.index(True)
        ipath[i0] = False
        tour = [i0]
        u = path[i0]
        if len(u) == 0: continue
        i1, i2 = u
        while i1 != i0:
            print(i0, i1, i2)
            tour.append(i1)
            ipath[i1] = False
            u = path[i1]
            if len(u) == 0:
                i0 = i1
                i1 = i2
            elif len(u) == 1:
                i1 = u[0]
            else:
                i1, i2 = u
        tours.append(tour)
    return tours

I also changed the constraint in the loop to this, but I guess I am wrong:

prob += pulp.lpSum([roads[i,j] for i in tour for j in range(i + 1, len(tour))]) <= len(tour) - 1

Is there some correct mathematical way to use just one variable to represent a road, and where am I wrong? How should the findTours function work in this case, or do I need to come up with a way to simplify the findPath function?

Solution

  • You’re missing the subtour elimination constraints, so in general, the roads chosen will form many connected components. Probably the best way to handle these is to generate them on demand: