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 in the thread where you got the subtour elimination code suggested.You can:
- use 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 in the thread where you got the subtour elimination code suggested. This is very expensive for large graphs.
- Use MTZ constraints.
I'm not sure about the correctness or performance but hereHere it is using MTZ constraints (disregarding the y constraints): code modified from Kuifje02/vrpy/: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
# 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]
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 = 10100 # 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)
)
# MTZ
M = len(G.nodes())
# Add rank varibles
z = m._xaddVars(
[(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 != x"Sink"),
name="Sink_after_%s" % v,
)
m.Paramsparams.LazyConstraintstimelimit = 160
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.
