Return to Revisions

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 subtours for a given forklift. I'm not sure about the correctness or perfomance but here 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 sub-tours
def subtourelim(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):
            selected = gp.tuplelist(
                (i, j, k) for i, j in G.edges() if vals[i, j, k] > 0.0
            )
            # find the shortest cycle in the selected edge list
            tours = allsubtours(selected)
            for tour in tours:
                tourlen = len(tour)
                tour.append(tour[0])
                if tourlen < n:
                    model.cbLazy(
                        gp.quicksum(
                            model._x[i, j, k]
                            for i, j in zip(tour, tour[1:])
                            if (i, j) in G.edges()
                        )
                        <= tourlen - 1
                    )
                    print("added lazy cut", expr)


# Given a tuplelist of edges, find all subtours
def allsubtours(edges):
    unvisited = list(range(n))
    cycles = []  # list(range(n + 1))  # initial length has 1 more location
    while unvisited:
        thiscycle = []
        neighbors = unvisited
        while neighbors:
            current = neighbors[0]
            thiscycle.append(current)
            unvisited.remove(current)
            neighbors = [
                j for i, j, _ in edges.select(current, "*", "*") if j in unvisited
            ]
        if len(thiscycle) > 1:
            cycles.append(thiscycle)
    print(f"{cycles=}")
    return cycles


n = 10  # 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 ["Sink"]:
        pred = list(G.predecessors(j))
        if len(pred) > 0:
            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(subtourelim)

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.