Skip to content

Callback: User cuts

Consider again the Vehicle Routing Problem with Time Windows and let's add some valid inequalities to the path MIP.

Recall that did_i is the demand of customer iCi \in C and the QQ is the vehicle capacity. The rounded capacity inequalities are

(i,j)δ(S)xijiSdiQ \sum_{(i,j) \in \delta(S)} x_{ij} \geq \bigg \lceil \frac {\sum_{i \in S} d_i} Q \bigg \rceil

for some set SCS \subseteq C with S2|S| \geq 2. The set δ(S)\delta(S) are the connecting edges between SS and it's complement set S\overline{S}

Use the callback routine to add cuts when violated. In this example we do a super simple enumeration of sets SS where S=2|S| = 2.

# Vehicle Routing Problem with Time Windows

import math
from flowty import Model, xsum, CallbackModel, Where
from flowty.datasets import vrp_rep

bunch = vrp_rep.fetch_vrp_rep("solomon-1987-r1", instance="R102_025")
name, n, es, c, d, Q, t, a, b, x, y = bunch["instance"]

m = Model()

# one graph, it is identical for all vehicles
g = m.addGraph(obj=c, edges=es, source=0, sink=n - 1, L=1, U=n - 2, type="B")

# adds resources variables to the graph.
# demand and capacity
m.addResourceDisposable(
    graph=g, consumptionType="V", weight=d, boundsType="V", lb=0, ub=Q, name="d"
)

# travel time and customer time windows
m.addResourceDisposable(
    graph=g, consumptionType="E", weight=t, boundsType="V", lb=a, ub=b, name="t"
)


# The callback for adding rounded capacity cuts
# x(delta(S)) >= floor( sum_(i \in S) d_i / Q ) ,  S \subset C, |S| >= 2
#
# We do a pretty stupid enumeration of all sets S with |S| = 2
# Do something smarter like CVRPSEP
def callback(cb: CallbackModel, where: Where):
    if where == Where.PathMIPCuts:
        epsilon = 1e-4
        relax = cb.x

        # Enumerate cuts for |S| = 2
        for i in range(n)[1:-1]:
            for j in range(n)[1:-1]:

                # not depot vertices in S
                if i == 0 or j == n - 1 or i == j:
                    continue

                rhs = math.ceil((d[i] + d[j]) / Q)

                lhs = 0
                lhsIdxs = []
                for e, edge in enumerate(es):
                    # inside S
                    if (i, j) == edge or (j, i) == edge:
                        continue

                    # if either end of edge is in S
                    if i == edge[0] or i == edge[1] or j == edge[0] or j == edge[1]:
                        lhs += relax[e]
                        lhsIdxs += [e]

                if lhs < rhs - epsilon:
                    cb.addCut(xsum(g.vars[e] for e in lhsIdxs) >= rhs)


m.setCallback(callback)

# set partition constriants
for i in range(n)[1:-1]:
    m += xsum(x * 1 for x in g.vars if i == x.source) == 1

# packing set
for i in range(n)[1:-1]:
    m.addPackingSet([x for x in g.vars if i == x.source])

status = m.optimize()
print(f"ObjectiveValue {round(m.objectiveValue, 1)}")

# get the variable values
for var in m.vars:
    if var.x > 0:
        print(f"{var.name} = {round(var.x, 1)}")