[Python] Wrong structure of tours Location Routing Problem gurobipy

I need to program a code for the LRP in Gurobipy. Unfortunately it is not working right but I can't find the mistake. It assigns my customers correctly to the depots but it does not route the vehicles right. In my output graph the vehicles approach the customers individually but there should be several customers served in one tour as long as the capacity of the vehicle is not exceeded (Which wouldn't be the case for some tours with my data). Can somebody maybe help me fix my code?

The dataset I am using has 3 depots (capacity 10000, opening costs 7500) and 10 customers with demands and transportation costs from the depots to the customers given. It looks like this:

3 10
10000 7500.
10000 7500.
10000 7500.
146
6739.72500 10355.05000 7650.40000
87
3204.86250 5457.07500 3845.40000
672
4914.00000 26409.60000 19622.40000
1337
32372.11250 29982.22500 21024.32500
31
1715.46250 2152.17500 1577.90000
559
6421.51250 23701.60000 16197.02500
2370
81972.37500 28499.25000 43134.00000
1089
33391.46250 26544.37500 6370.65000
33
2020.83750 2480.77500 1869.45000
32
1459.60000 1995.20000 1402.40000


And here is the code:

from gurobipy import Model, GRB, quicksum
import math
import networkx as nx
import matplotlib.pyplot as plt

# ------------------------- Data input ----------------------------
def parse_allocation_format(filepath):
with open(filepath, 'r') as f:
lines = [line.strip() for line in f if line.strip()]

m, n = map(int, lines[0].split())
I = list(range(m)) # Depots
J = list(range(n)) # Kunden

W, O = {}, {}
for i, line in enumerate(lines[1:1 + m]):
capacity, fixed_cost = map(float, line.split())
W = capacity
O = fixed_cost

d = {}
allocation_costs = {}
j_index = 0
i = 1 + m
while i < len(lines):
if len(lines.split()) == 1:
d[j_index] = float(lines)
allocation_costs[j_index] = list(map(float, lines[i+1].split()))
j_index += 1
i += 2
else:
i += 1

return I, J, d, W, O, allocation_costs

# ------------------------- Modell ----------------------------
filepath = "C:/Users" # Adjust path
I, J, d, W, O, allocation_costs = parse_allocation_format(filepath)

m = Model("CLRP")

Q = 2500 # Vehicle capacity
F = 1500 # Fixed costs per vehicle

V = list(range(len(I))) + [j + len(I) for j in J] # Depots + Customers

# Kost matrix
c = {(i, j + len(I)): allocation_costs[j] for j in J for i in I}

# Decision variables
y = m.addVars(I, vtype=GRB.BINARY, name="y") # Indicates that depot i is opened
w = m.addVars(J, I, vtype=GRB.BINARY, name="w") # Indicates that customer j is served by depot i
x = m.addVars(((i, j) for i in V for j in V if i != j), vtype=GRB.BINARY, name="x") # Indicates that vehicle travels from i to j
u = m.addVars(J, vtype=GRB.CONTINUOUS, name="u") # Additional variable for MTZ

# Objective function
m.setObjective(
quicksum(x[i, j] * c[i, j] if (i in I and j in [k + len(I) for k in J]) else 0 for i, j in x.keys()) +
quicksum(x[i, j] * F for i, j in x.keys()) +
quicksum(y * O for i in I),
GRB.MINIMIZE
)

# Constraints
for j in J:
cust = j + len(I)
m.addConstr(quicksum(x[i, cust] for i in V if i != cust) == 1) # Constraint: 2 -> Every customer served by exactly on vehicle
m.addConstr(quicksum(x[cust, i] for i in V if i != cust) == 1) # Constraint: 3 -> Vehicle leaves every customer again
m.addConstr(quicksum(w[j, i] for i in I) == 1) # Constraint: 4 -> Every customer assigned to exactly one depot

# Constraint: 5 -> Customers only assigned to open depots
for j in J:
for i in I:
m.addConstr(w[j, i] <= y)

# Constraint: 6 -> Capacty constraint for depot
for i in I:
m.addConstr(quicksum(w[j, i] * d[j] for j in J) <= W)

# Constraints: 7-10:
for i in I:
for j in J:
cust = j + len(I)
m.addConstr(x[i, cust] <= w[j, i]) # Constraint: 9 -> arc from i to j only exists if customer assigned to depot (route & starts ends at depot)
m.addConstr(x[cust, i] <= w[j, i]) # Constraint: 10 -> arc from j to i only exists if customer assigned to depot

for i in I:
for j in J:
for k in J:
if j != k:
m.addConstr(x[j + len(I), k + len(I)] <= 1 - w[j, i] + w[k, i]) # Constraint: 7 -> only one arc between customers if both served by same depot
m.addConstr(x[j + len(I), k + len(I)] <= 1 + w[j, i] - w[k, i]) # Constraint: 8 -> vehicles only serve customers assigned to same depot


# MTZ subtour elimination + vehicle capacity
for j in J:
m.addConstr(u[j] >= d[j])
m.addConstr(u[j] <= Q)

for j in J:
for k in J:
if j != k:
m.addConstr(u[j] + d[k] - Q * (1 - x[j + len(I), k + len(I)]) <= u[k])

# Initialisation MTZ from Depot
for i in I:
for j in J:
m.addConstr(u[j] >= d[j] * x[i, j + len(I)])

# Rout starts and ends at depot
for j in J:
m.addConstr(
quicksum(x[i, j + len(I)] for i in I) ==
quicksum(x[j + len(I), i] for i in I)
)

# ------------------------- Optimisation ----------------------------
m.optimize()

# ------------------------- Visualisation ----------------------------
G = nx.DiGraph()
for i, j in x.keys():
if x[i, j].X > 0.5:
G.add_edge(i, j)

pos = nx.spring_layout(G, seed=42)

active_nodes = set(G.nodes)
depots = [i for i in I if i in active_nodes]
customers = [j + len(I) for j in J if (j + len(I)) in active_nodes]

nx.draw_networkx_nodes(G, pos, nodelist=depots, node_color='red', node_shape='s', node_size=700)
nx.draw_networkx_nodes(G, pos, nodelist=customers, node_color='skyblue', node_shape='o', node_size=700)
nx.draw_networkx_edges(G, pos, arrows=True)
nx.draw_networkx_labels(G, pos)

plt.title("CLRP-Solution: Depots (red), Customers (blue)")
plt.axis('off')
plt.show()

Continue reading...