The subtour elimination constraints are not correct, they are mostly empty. Try to debug the callback calls, as you can see also from the log, none of lazy constraints are applied. I managed to get it to work after some modifications to the callbacks to enumerate and forbid all cycles for a given forklift (maximum of `max_no_cycles` for each `k`). This is what the [last comment][1] in the thread where you got the subtour elimination code suggested. I'm not sure about the correctness or performance but here it is (disregarding the `y` constraints): ``` import math import random import gurobipy as gp from gurobipy import GRB import networkx as nx import matplotlib.pyplot as plt # Callback - use lazy constraints to eliminate all cycles def callback(model, where): if where == GRB.Callback.MIPSOL: # make a list of edges selected in the solution vals = model.cbGetSolution(model._x) for k in range(K): sub_graph = nx.DiGraph() [sub_graph.add_edge(i, j) for i, j in G.edges() if vals[i, j, k] > 0.0] # find the shortest cycle in the selected edge list count = 0 for cycle in nx.simple_cycles(sub_graph): RHS = len(cycle) - 1 cycle.append(cycle[0]) if len(cycle) < n: model.cbLazy( gp.quicksum( model._x[i, j, k] for i, j in zip(cycle, cycle[1:]) if (i, j) in G.edges() ) <= RHS ) count += 1 if count > max_no_cycles: break n = 10 # number of items to pick, equivalent to number of locations to visit K = 3 # number of fork-lifts to use # Warning changing this to a large number this will take a very long time for large graphs max_no_cycles = 100 # limit number of cycles used in callback # Create n random points points = [(0, 0)] points += [(random.randint(0, 100), random.randint(0, 100)) for i in range(n - 1)] # Dictionary of Manhattan distance between each pair of points dist = { (i, j): math.sqrt(sum((points[i][p] - points[j][p]) ** 2 for p in range(2))) for i in range(n) for j in range(n) if i != j } # Create graph G = nx.DiGraph() for k, v in dist.items(): if k[0] == 0: i = "Source" else: i = k[0] if k[1] == 0: j = "Sink" else: j = k[1] G.add_edge(i, j, dist=v) m = gp.Model() # Create variables: x_keys = {(e[0], e[1], k): e[2]["dist"] for e in G.edges(data=True) for k in range(K)} x = m.addVars( x_keys, obj=x_keys, vtype=GRB.BINARY, name="x", ) # Visit all nodes for j in G.nodes(): if j not in ["Source", "Sink"]: pred = list(G.predecessors(j)) m.addConstr(gp.quicksum(x[i, j, k] for i in pred for k in range(K)) >= 1) # Flow-balance for v in G.nodes(): if v not in ["Source", "Sink"]: m.addConstrs( gp.quicksum(x[i, v, k] for i in G.predecessors(v)) - gp.quicksum(x[v, j, k] for j in G.successors(v)) == 0 for k in range(K) ) # All k's must start at Source m.addConstrs( gp.quicksum(x["Source", j, k] for j in G.successors("Source")) == 1 for k in range(K) ) # All k's must end at Sink m.addConstrs( gp.quicksum(x[i, "Sink", k] for i in G.predecessors("Sink")) == 1 for k in range(K) ) m._x = x m.Params.LazyConstraints = 1 m.optimize(callback) for k in range(K): print(f"{k=}") for i, j in G.edges(): if x[i, j, k].X > 0.0: print(f"\t{i}->{j}") ``` PS. Sorry, I didn't read the complete thread, I thought the question was resolved. [1]: https://support.gurobi.com/hc/en-us/community/posts/4410235441681/comments/5553651253137