PuLP求解带节点时间成本的VRP时约束冲突致不可行问题排查

问题背景 #

使用PuLP调用GLPK求解含节点额外时间成本的VRP问题，已将节点时间成本直接并入旅行时间矩阵。模型约束包含：每个客户仅被一名员工访问一次、每位员工路线首尾为节点0，但求解器返回问题不可行。经逐一隔离约束排查，确认约束1是导致不可行的核心原因。

原代码 #

import pulp
from pulp import LpMaximize, LpProblem, LpStatus, lpSum, LpVariable, LpMinimize, LpBinary, LpInteger, GLPK
input = []
with open("test.txt", "r") as f:
    for line in f:
        input.append(line)

n = [int(x) for x in input[0].split()]
N = n[0] #number of customer
K = n[1] #number of staff
c = [int(x) for x in input[1].split()] #time cost at each customer vector 
c.insert(0, 0)
poi_ma = [] # distance matrix
for i in range(N+1):
    poi_ma.append([int(x) for x in input[i+2].split()])
    for y in range(0, N + 1):
        if y == i:
            continue
        else:
            poi_ma[i][y] = poi_ma[i][y] + c[y]

model = LpProblem("myprob",sense=LpMinimize)
X =[]
X.append(0)
for k in range(1, K+1):
    X.append([])
    for i in range(N+1):
        X[k].append([])
        for j in range(N+1):
            X[k][i].append(LpVariable(f"X[{k}][{i}][{j}]", cat=LpBinary))

y = LpVariable("y", cat= LpInteger, lowBound = 0)
U = [0]
for i in range(1, N+1):
    U.append(LpVariable(f"U[{i}]",lowBound = 1, upBound = 10000))


# 1 - Each customer is served by exactly one staff member
for i in range(N+1):
    for j in range(N+1):
        if i != j:
            model += (lpSum([X[k][i][j] for k in range(1, K+1)]) == 1)
        else:
            model += (lpSum([X[k][i][j] for k in range(1, K+1)]) == 0)


# 2 - Each staff member starts and ends at the headquarters
for k in range(1, K+1):
    model += (lpSum([X[k][0][j] for j in range(N+1)]) == 1)
    model += (lpSum([X[k][x][0] for x in range(N+1)]) == 1)

# 3 - Enforce flow conservation at each customer location
for k in range(1, K+1):
    for i in range(N + 1):
        tmp1 = lpSum(X[k][i][j] for j in range (N + 1))
        tmp2 = lpSum(X[k][j][i] for j in range (N + 1))
        model += (tmp1 <= 1)
        model += (tmp2 <= 1)
        model += (tmp1 - tmp2 == 0)


# 4 - Miller-Tucker-Zemlin formulation
for i in range(1, N+1):
    #model += (1 <= U[i])
    for j in range(1, N+1):
        if j != i:
            for k in range(1, K+1):
                model += (1 - 10000 + 10000 * X[k][i][j]) <= (U[j] - U[i])

# 5 - Calculate work time for each staff member and ensure it's less than or equal to y
for k in range(1, K+1):
    work_time = lpSum([poi_ma[i][j] * X[k][i][j] for i in range(N + 1) for j in range(N + 1) if i != j])
    model += (work_time <= y)

model += y
status = model.solve(pulp.PULP_CBC_CMD(msg=True, cuts=True)) 

print(model.status)
for k in range(1, K+1):
    print(f"{k}", end = ':')
    for i in range(N+1):
        for j in range(N+1):
            if X[k][i][j].varValue == 1:
                print(f"{i},{j}", end = ' ')
    print()
print(y.varValue)

test.txt输入内容 #

5 2
6 8 7 1 9
0 5 10 6 4 8
5 0 5 4 2 6
10 5 0 5 7 4
6 4 5 0 6 2
4 2 7 6 0 8
8 6 4 2 8 0

运行报错信息 #

Local\Temp\086601480d8341f89343bb8327bacee3-pulp.mps gomory on knapsack on probing on timeMode elapsed branch printingOptions all solution C:\Users\hungn\AppData\Local\Temp\086601480d8341f89343bb8327bacee3-pulp.sol (default strategy 1)
At line 2 NAME          MODEL
At line 3 ROWS
At line 123 COLUMNS
At line 825 RHS
At line 944 BOUNDS
At line 1028 ENDATA
Problem MODEL has 118 rows, 78 columns and 542 elements
Coin0008I MODEL read with 0 errors
Option for gomoryCuts changed from ifmove to on
Option for knapsackCuts changed from ifmove to on
Option for timeMode changed from cpu to elapsed
Problem is infeasible - 0.00 seconds
Option for printingOptions changed from normal to all
Total time (CPU seconds):       0.01   (Wallclock seconds):       0.00

-1
1:0,3 1,0 2,0 2,3 3,0 4,3 5,0 5,2 5,3
2:0,1 0,2 0,4 0,5 1,2 1,3 1,4 1,5 2,1 2,4 3,1 3,2 3,4 3,5 4,0 4,1
346.0

排查结果 #

逐一隔离约束测试后，仅保留约束1时问题仍不可行，确认约束1是导致模型不可行的关键因素。

问题定位与修复 #

约束1的错误分析 #

原约束1的逻辑完全错误：

for i in range(N+1):
    for j in range(N+1):
        if i != j:
            model += (lpSum([X[k][i][j] for k in range(1, K+1)]) == 1)

该约束要求每一对不同的节点(i,j)，都必须有且仅有一名员工走这条弧。对于N=5的场景，共有6*5=30对i≠j的节点组合，这意味着模型需要生成30条运输弧，远超过VRP实际需要的5次客户访问量，直接导致矛盾，模型不可行。

修正后的约束1 #

正确的约束逻辑是：每个客户节点（i∈1~N）被恰好访问一次，即所有员工到达该客户的弧之和为1（或离开该客户的弧之和为1，结合流量守恒可自动保证入出平衡）。同时保留禁止自环的约束：

# 修正后的约束1 - 每个客户仅被访问一次，禁止自环
# 禁止所有员工走自环弧
for i in range(N+1):
    model += (lpSum([X[k][i][i] for k in range(1, K+1)]) == 0)

# 每个客户节点（1~N）被恰好访问一次（通过入弧之和约束）
for i in range(1, N+1):
    model += (lpSum([X[k][j][i] for k in range(1, K+1) for j in range(N+1) if j != i]) == 1)

额外优化建议 #

1. MTZ约束的大M值优化：原代码中M取10000过大，建议取N+1（如6），避免数值计算问题，修正后的MTZ约束：

M = N + 1  # 取节点总数+1即可
for i in range(1, N+1):
    for j in range(1, N+1):
        if j != i:
            for k in range(1, K+1):
                model += U[j] >= U[i] + 1 - M * (1 - X[k][i][j])

2. 流量守恒约束简化：原代码中tmp1 <=1和tmp2 <=1可省略，因为客户节点的入弧之和已被约束为1，结合流量守恒tmp1=tmp2，自然保证入出弧数量不超过1；总部节点0的入出弧数量由约束2控制为1，无需额外限制。

内容的提问来源于stack exchange，提问作者alksdhalksjdb