How to extend TSP to MTSP using Pulp

Question

We've studied TSP and now we're tasked to extend it for multiple salespersons. Below code using PULP with my added logic which unfortunately does not work. It does not correctly identify the correct tours.

Eg. using below cost matrix:

cost_matrix = [[ 0, 1, 3, 4],
     [ 1, 0, 2, 3 ], 
     [ 3, 2, 0, 4 ], 
     [ 4, 3, 4, 0]]

and n = len(cost_matrix)

for salespersons (k = 3) the results should be:

The path of SP_1 is: 0 => 1 => 0

The path of SP_2 is: 0 => 2 => 3 => 0

Can someone help me solve this problem?

    # create encoding variables
    bin_vars = [ # add a binary variable x_{ij} if i not = j else simply add None
        [ LpVariable(f'x_{i}_{j}', cat='Binary') if i != j else None for j in range(n)] 
        for i in range(n) ]
    time_stamps = [LpVariable(f't_{j}', lowBound=0, upBound=n, cat='Continuous') for j in range(1, n)]
    # create add the objective function
    objective_function = lpSum( [ lpSum([xij*cj if xij != None else 0 for (xij, cj) in zip(brow, crow) ])
                           for (brow, crow) in zip(bin_vars, cost_matrix)] )
    
    prob += objective_function 

    # add constraints
    for i in range(n):
        # Exactly one leaving variable
        prob += lpSum([xj for xj in bin_vars[i] if xj != None]) == 1
        # Exactly one entering
        prob += lpSum([bin_vars[j][i] for j in range(n) if j != i]) == 1
    
    # add timestamp constraints
    for i in range(1,n):
        for j in range(1, n):
            if i == j: 
                continue
            xij = bin_vars[i][j]
            ti = time_stamps[i-1]
            tj = time_stamps[j -1]
            prob += tj >= ti + xij - (1-xij)*(n+1)

    
    # Binary variables to ensure each node is visited by a salesperson
    visit_vars = [LpVariable(f'u_{i}', cat='Binary') for i in range(1, n)]
    
    # Salespersons constraints
    prob += lpSum([bin_vars[0][j] for j in range(1, n)]) == k
    prob += lpSum([bin_vars[i][0] for i in range(1, n)]) == k

    for i in range(1, n):
        prob += lpSum([bin_vars[i][j] for j in range(n) if j != i]) == visit_vars[i - 1]
        prob += lpSum([bin_vars[j][i] for j in range(n) if j != i]) == visit_vars[i - 1]
    

    # Done: solve the problem
    status = prob.solve(PULP_CBC_CMD(msg=False))

Answer 1

The function needs fixing and does not [require] total refactoring.

Incorrect. Among other issues,

• you're failing to use LpVariable.matrix when appropriate;
• you're failing to name your constraints (critical during debugging); and
• rather than None for diagonal elements, you should just set 0 and then treat the elements unconditionally.

As I said in the comments, the number one issue with the code is that the sales constraints don't make sense and render the problem infeasible. I demonstrate without them. The following does run and solve, though it's dubious whether it produces the results you want - there's not particularly a representation of separate salespeople, and if there needs to be, you need to re-design the problem construction.

import pulp

cost_matrix = (
    (0, 1, 3, 4),
    (1, 0, 2, 3),
    (3, 2, 0, 4),
    (4, 3, 4, 0),
)
n = len(cost_matrix)
k = 3

encoding = [
    [
        0 if i == j else pulp.LpVariable(name=f'x_{i}_{j}', cat=pulp.LpBinary)
        for j in range(n)
    ]
    for i in range(n)
]
time_stamps = pulp.LpVariable.matrix(
    name='t', lowBound=0, upBound=n, cat=pulp.LpContinuous,
    indices=range(1, n),
)

# Binary variables to ensure each node is visited by a salesperson
visited = pulp.LpVariable.matrix(
    name='u', cat=pulp.LpBinary, indices=range(1, n),
)

objective_function = pulp.lpSum(
    pulp.lpDot(brow, crow)
    for (brow, crow) in zip(encoding, cost_matrix)
)

prob = pulp.LpProblem()
prob.setObjective(objective_function)

for i, brow in enumerate(encoding):
    # Exactly one leaving variable
    prob.addConstraint(
        name=f'leaving_{i}_j', constraint=pulp.lpSum(brow) == 1)
    # Exactly one entering
    prob.addConstraint(
        name=f'entering_i_{i}',
        constraint=pulp.lpSum(encoding[j][i] for j in range(n)) == 1,
    )

for i, brow in enumerate(encoding[1:], start=1):
    for j, xij in enumerate(brow[1:], start=1):
        ti = time_stamps[i - 1]
        tj = time_stamps[j - 1]
        prob.addConstraint(
            name=f'time_{i}_{j}',
            constraint=tj >= ti + xij - (1 - xij)*(n + 1),
        )

# These constraints render the problem infeasible.
# prob.addConstraint(
#     name='sales_0_j',
#     constraint=pulp.lpSum(encoding[0][j] for j in range(1, n)) == k,
# )
# prob.addConstraint(
#     name='sales_i_0',
#     constraint=pulp.lpSum(encoding[i][0] for i in range(1, n)) == k,
# )

for i in range(1, n):
    prob.addConstraint(
        name=f'visit_{i}_j',
        constraint=pulp.lpSum(encoding[i][j] for j in range(n)) == visited[i - 1],
    )
    prob.addConstraint(
        name=f'visit_i_{i}',
        constraint=pulp.lpSum(encoding[j][i] for j in range(n)) == visited[i - 1],
    )

print(prob)
prob.solve()
if prob.status != pulp.LpStatusOptimal:
    raise ValueError(prob.status)

print('Encoding:')
for row in encoding:
    print([int(pulp.value(e)) for e in row])

print('\nTimestamp:')
print([pulp.value(t) for t in time_stamps])

print('\nVisited:')
print([int(pulp.value(v)) for v in visited])

MINIMIZE
1*x_0_1 + 3*x_0_2 + 4*x_0_3 + 1*x_1_0 + 2*x_1_2 + 3*x_1_3 + 3*x_2_0 + 2*x_2_1 + 4*x_2_3 + 4*x_3_0 + 3*x_3_1 + 4*x_3_2 + 0
SUBJECT TO
leaving_0_j: x_0_1 + x_0_2 + x_0_3 = 1

entering_i_0: x_1_0 + x_2_0 + x_3_0 = 1

leaving_1_j: x_1_0 + x_1_2 + x_1_3 = 1

entering_i_1: x_0_1 + x_2_1 + x_3_1 = 1

leaving_2_j: x_2_0 + x_2_1 + x_2_3 = 1

entering_i_2: x_0_2 + x_1_2 + x_3_2 = 1

leaving_3_j: x_3_0 + x_3_1 + x_3_2 = 1

entering_i_3: x_0_3 + x_1_3 + x_2_3 = 1

time_1_1: 0 t_1 >= -5

time_1_2: - t_1 + t_2 - 6 x_1_2 >= -5

time_1_3: - t_1 + t_3 - 6 x_1_3 >= -5

time_2_1: t_1 - t_2 - 6 x_2_1 >= -5

time_2_2: 0 t_2 >= -5

time_2_3: - t_2 + t_3 - 6 x_2_3 >= -5

time_3_1: t_1 - t_3 - 6 x_3_1 >= -5

time_3_2: t_2 - t_3 - 6 x_3_2 >= -5

time_3_3: 0 t_3 >= -5

visit_1_j: - u_1 + x_1_0 + x_1_2 + x_1_3 = 0

visit_i_1: - u_1 + x_0_1 + x_2_1 + x_3_1 = 0

visit_2_j: - u_2 + x_2_0 + x_2_1 + x_2_3 = 0

visit_i_2: - u_2 + x_0_2 + x_1_2 + x_3_2 = 0

visit_3_j: - u_3 + x_3_0 + x_3_1 + x_3_2 = 0

visit_i_3: - u_3 + x_0_3 + x_1_3 + x_2_3 = 0

VARIABLES
t_1 <= 4 Continuous
t_2 <= 4 Continuous
t_3 <= 4 Continuous
0 <= u_1 <= 1 Integer
0 <= u_2 <= 1 Integer
0 <= u_3 <= 1 Integer
0 <= x_0_1 <= 1 Integer
0 <= x_0_2 <= 1 Integer
0 <= x_0_3 <= 1 Integer
0 <= x_1_0 <= 1 Integer
0 <= x_1_2 <= 1 Integer
0 <= x_1_3 <= 1 Integer
0 <= x_2_0 <= 1 Integer
0 <= x_2_1 <= 1 Integer
0 <= x_2_3 <= 1 Integer
0 <= x_3_0 <= 1 Integer
0 <= x_3_1 <= 1 Integer
0 <= x_3_2 <= 1 Integer

Welcome to the CBC MILP Solver 
Version: 2.10.3 
Build Date: Dec 15 2019 

Result - Optimal solution found

Objective value:                11.00000000
Enumerated nodes:               0
Total iterations:               0
Time (CPU seconds):             0.01
Time (Wallclock seconds):       0.00

Option for printingOptions changed from normal to all
Total time (CPU seconds):       0.01   (Wallclock seconds):       0.01

Encoding:
[0, 0, 1, 0]
[1, 0, 0, 0]
[0, 0, 0, 1]
[0, 1, 0, 0]

Timestamp:
[4.0, 0.0, 1.0]

Visited:
[1, 1, 1]