Pulp constraints conflict

Question

I'm trying to use PuLP's GLPK solver to solve VRP with additional time cost at each node, which is added directly to the travel time matrix. Here is the code with the constraints to ensure each customer is visited once by 1 staff and all route for each staff start and end at node 0.

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)

the "test.txt" file contains:

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

When I tried to execute the code, it return this message:

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

which indicates the problem is infeasible. I have tried to isolate each constraint to evaluate them as in the image:

.

It is very likely that constraint 1 has a problem, but I can't find what it is. Can someone help me identify the problem here ? I appreciate all your help and opinion.

Answer 1

Your first constraint is incorrect. You should be summing over i or j along with k inside of the summation statement. Right now it says that one staff member should go from each i to each j, which is obviously not what you want. You want one staff member to arrive at each j from some i.

Also, your MTZ constraint needs to be indexed by k as well as each member of the staff needs their own sequence numbering.

Try those things out. If it is still infeasible, try commenting out individual constraints and running it, that should give you some clues.

If you are still stuck after trying that, edit your post to include the updated effort, and comment back.