在Python中定义旅行商线性规划模型

6
使用Python和PuLP库,我们如何创建线性规划模型来解决旅行商问题(TSP)? 从维基百科得知,目标函数和约束条件为:

enter image description here

问题:这是我的部分尝试,我卡住了。

  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 代码

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

数据文件可以在这里找到

结果

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



...

更新的代码

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

你的代码出了什么问题? - kindall
@kindall 我没有包含最终的限制条件,因为我不知道如何定义它。此外,解决当前模型会得到决策变量,例如 x1_1 等于 1.0,这显然是错误的...我无法弄清楚为什么会这样。 - Nyxynyx
@kindall 更新了问题以澄清问题,并在当前代码中包含了错误的结果。 - Nyxynyx
您无法使用此方法解决更大的问题,因为由于子旅行限制条件,公式将变得非常庞大。您可以考虑另一种动态分离无效解决方案的方法。这里有一个Python示例:https://github.com/SCIP-Interfaces/PySCIPOpt/blob/master/tests/test_tsp.py - mattmilten
这可能对您有所帮助:https://nbviewer.jupyter.org/github/cochoa0x1/intro-mathematical-programming/blob/master/05-routes-and-schedules/traveling_salesman.ipynb - Akavall
显示剩余3条评论
1个回答

2
最后一个约束条件并不是单一的约束条件。对于每对满足该条件的索引,您应该添加一个约束条件:
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

然而,你需要首先将u_i变量声明为整数,通过将cat='Integer'参数传递给LpVariable

u = []
for i in range(n):
    u.append(pulp.LpVariable('u_' + str(i), cat='Integer'))

谢谢,我正在研究你的建议。目前的代码看起来还好吗?最后一个约束条件的缺失是导致 x0_0x1_1 等被求解成 1.0 而不是 0.0 的原因吗? - Nyxynyx
@Nyxynyx 最后一个限制条件对于这个问题来说是绝对基本的。如果没有它,这个问题甚至不是一个整数线性问题,而只是一个线性问题,因此无法以任何有意义的方式表示TSP(请记住,TSP无法在多项式时间内解决,甚至无法在多项式时间内以“良好的方式”近似,因此每个非整数线性规划问题(总是可以在多项式时间内解决)都不能代表它或任何良好的近似)。 - Bakuriu
由于我正在尝试吸收最后一个约束条件的代码,所以我将它们添加到当前代码中并尝试解决它,但结果仍然看起来很奇怪,其中变量x [i,j]的值为1.0,其中 i == j - Nyxynyx
@Nyxynyx 你应该尝试使用for i in range(n+1): for j in range(n+1)来定义你的x变量,并修改所有约束条件以从1开始,所以即使是最后一个也应该是for i in range(1, n): for j in range(i+1, n+1): ...。此外,没有u_0变量,因此在定义u向量时需要进行调整...只需在开头添加一个虚拟条目:u = [None]; for i in range(1, n+1): u.append(...)这样u[i]就指代索引为i的城市,而你不必在最终约束条件中添加-1 - Bakuriu
@Nyxynyx 如果您仔细查看目标函数,您会发现 j 索引的范围从 1n,并且还有一个约束条件 j != i,因此您应该根本不会计算出首先的 x[0,0] - Bakuriu
显示剩余5条评论

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接