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.

You can: 
- use callbacks to enumerate and forbid all cycles for a given forklift. This is what the [last comment][1] in the thread where you got the subtour elimination code suggested. This is very expensive for large graphs. 
- Use MTZ constraints.

For more in-depth study of the SDVRP, please see [Dror et al. (1994)][2].

Here it is using MTZ constraints (disregarding the `y` constraints): code modified from [Kuifje02/vrpy/:vrpy/subproblem_lp.py#L227](https://github.com/Kuifje02/vrpy/blob/master/vrpy/subproblem_lp.py#L227): 

```
import math
import random

import gurobipy as gp
from gurobipy import GRB
import networkx as nx
import matplotlib.pyplot as plt

n = 100  # number of items to pick, equivalent to number of locations to visit
K = 3  # number of fork-lifts to use

# 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)
)

# MTZ
M = len(G.nodes())

# Add rank varibles
z = m.addVars(
    [(v, k) for v in G.nodes() for k in range(K)],
    name="z",
    lb=0,
    ub=len(G.nodes()),
    vtype=GRB.INTEGER,
)
# Add big-M constraints
m.addConstrs(
    (
        z[i, k] + 1 <= z[j, k] + M * (1 - x[i, j, k])
        for k in range(K)
        for (i, j) in G.edges()
    ),
    name="elementary",
)

# Source is first, Sink is last (optional)
m.addConstrs((z["Source", k] == 0 for k in range(K)), name="Source_is_first")
m.addConstrs(
    (z[v, k] <= z["Sink", k] for k in range(K) for v in G.nodes() if v != "Sink"),
    name="Sink_after_%s" % v,
)

m.params.timelimit = 60
m.optimize()

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}")
```


  [1]: https://support.gurobi.com/hc/en-us/community/posts/4410235441681/comments/5553651253137
  [2]: https://www.sciencedirect.com/science/article/pii/0166218X9200172I?ref=cra_js_challenge&fr=RR-1