3.4 Three classic models
Class: implementation · Script: python/cap06_bpp.py, python/cap06_cmax.py, python/cap06_tsp.py
Bin packing, makespan and travelling salesman: statement, model, construction in
gurobipy and model of the instance. They are the three problems on which the
heuristics chapter builds next-fit, first-fit, best-fit, LPT
and nearest neighbour.
So far the example model has always been the knapsack. The three problems below come back in the heuristics chapter, where the solutions of next-fit, first-fit, best-fit, LPT and nearest neighbour are built by hand: here their models are written, so that those heuristics have an optimum to compare themselves with.
Bin packing: how many containers are enough
The script is python/cap06_bpp.py.
Bin packing
There are \(n\) items, item \(j\) weighs \(w_j\). The containers are all alike, of capacity \(C\). Use the smallest number of containers.
The model
Two families of binary variables are needed: \(x_{jb} = 1\) if item \(j\) goes into container \(b\), and \(y_b = 1\) if container \(b\) is used.
The first family says that every item ends up in exactly one container. The second is the capacity written as an activation: while \(y_b = 0\) container \(b\) can receive nothing, and as soon as \(y_b = 1\) it takes up to \(c\). The objective counts the containers switched on.
Building it in gurobipy
def modello_bpp(w, c, k):
n = len(w)
m = nuovo_modello("bin_packing")
x = m.addVars(n, k, vtype=GRB.BINARY, name="x")
y = m.addVars(k, vtype=GRB.BINARY, name="y")
m.setObjective(y.sum(), GRB.MINIMIZE)
m.addConstrs((x.sum(j, "*") == 1 for j in R(n)), name="item")
m.addConstrs((gp.quicksum(w[j] * x[j, b] for j in R(n)) <= c * y[b]
for b in R(k)), name="capacity")
return m, x, y
The instance
On the instance of four items of weight \(w = (5, 4, 3, 3)\) and capacity \(c = 7\):
The total weight is \(15\), so no solution can use fewer than \(\lceil 15/7 \rceil = 3\) containers; the optimum uses exactly \(3\), and the count is therefore tight.
The relaxation of bin packing is very weak
Relaxing \(y_b\) to \(y_b \ge 0\) the model buys fractions of a container, and the LP optimum drops to \(\sum_j w_j / c = 15/7 \approx 2.14\): the relaxation does not know that a container opens whole. This is why on this problem the dual bound is of little use and the heuristics matter more.
Makespan: the load of the busiest machine
The script is python/cap06_cmax.py.
Makespan on identical machines
There are \(n\) jobs, of duration \(d_j\), and \(k\) identical machines. Every job goes on a single machine and is not interrupted. Minimise the load of the busiest machine.
The model
With \(x_{jm} = 1\) if job \(j\) goes on machine \(m\), and \(z \ge 0\) the load of the busiest machine:
The objective is the single variable \(z\): no datum appears in it. It is the \(k\) load rows that give it meaning, saying that no machine works longer than \(z\); the minimisation then presses it down onto the load of the busiest machine. It is the min-max technique.
Building it in gurobipy
def modello_cmax(d, k):
n = len(d)
m = nuovo_modello("makespan")
x = m.addVars(n, k, vtype=GRB.BINARY, name="x")
z = m.addVar(name="z")
m.setObjective(z, GRB.MINIMIZE)
m.addConstrs((x.sum(j, "*") == 1 for j in R(n)), name="job")
m.addConstrs((gp.quicksum(d[j] * x[j, mm] for j in R(n)) <= z
for mm in R(k)), name="load")
return m, x, z
The instance
On the instance of four jobs of duration \(d = (3, 4, 5, 6)\) on \(k = 2\) machines — the very one the heuristics chapter runs LPT on:
The total load is \(18\) and the machines are two: no solution can go below \(18/2 = 9\), and the optimum is exactly \(9\) — the jobs split into \(6+3\) and \(5+4\). Here the count settles the problem on its own.
Travelling salesman: the MTZ formulation
The script is python/cap06_tsp.py.
Travelling salesman
There are \(n\) cities and a distance \(d_{ij}\) between every pair. Find the shortest tour that visits every city exactly once and returns to the starting point.
The model
With \(x_{ij} = 1\) if the tour goes from \(i\) to \(j\), the two families "one leaves once" and "one enters once" are not enough: they also admit solutions made of separate subtours. The Miller–Tucker–Zemlin formulation adds a variable \(u_i\) for every city other than the first, recording its position along the tour.
The third group is the heart of the formulation. If \(x_{ij} = 0\) the row becomes \(u_i - u_j \le n - 1\), always true because the \(u\) lie between \(1\) and \(n-1\): it forbids nothing. If instead \(x_{ij} = 1\) it becomes \(u_j \ge u_i + 1\), that is "if I go from \(i\) to \(j\), the position of \(j\) is the next one". A subtour missing city \(1\) would need a chain of ever-growing positions closing on itself, which is impossible; city \(1\) has no \(u\) of its own precisely because it is where the tour closes.
Building it in gurobipy
def modello_tsp(D):
n = len(D)
m = nuovo_modello("tsp")
x = m.addVars(((i, j) for i in R(n) for j in R(n) if i != j),
vtype=GRB.BINARY, name="x")
u = m.addVars(R(1, n), lb=1, ub=n - 1, name="u")
m.setObjective(gp.quicksum(D[i][j] * x[i, j] for i, j in x), GRB.MINIMIZE)
m.addConstrs((gp.quicksum(x[i, j] for j in R(n) if j != i) == 1
for i in R(n)), name="out")
m.addConstrs((gp.quicksum(x[i, j] for i in R(n) if i != j) == 1
for j in R(n)), name="in")
m.addConstrs((u[i] - u[j] + n * x[i, j] <= n - 1
for i in R(1, n) for j in R(1, n) if i != j), name="mtz")
return m, x, u
The instance
On the four-city instance of the heuristics chapter the optimal tour is \(1 \to 2 \to 4 \to 3 \to 1\) and measures \(22\). The model of the instance has fifteen columns — twelve arcs and three positions — and fourteen rows:
MTZ is convenient, not the strongest
The MTZ constraints are \(O(n^2)\) and take three lines of gurobipy, but their
linear relaxation is weak: the continuous \(u\) absorb almost everything and the
LP stays far from the integer optimum. Formulations that remove subtours with
connectivity cuts give much better bounds, at the price of an exponential
number of constraints to be generated as they are needed. For the sizes of
this course MTZ is enough.
Show the complete script — python/cap06_bpp.py (51 lines)
"""Bin packing: how many containers are enough (chapter 3).
The first of the three problems the heuristics chapter takes up again: there
next-fit, first-fit and best-fit are built by hand, here the model that says what
the optimum is. The capacity is written as an activation: until the container
opens it can receive nothing.
"""
import gurobipy as gp
import pandas as pd
from gurobipy import GRB
from esteso import salva_modello
from mip import frazione, nuovo_modello, risolvi
from stile import intestazione, salva_dati
R = range
intestazione("Bin packing: the smallest number of containers")
w_bpp = [5, 4, 3, 3] # weight of the items
c_bpp = 7 # capacity of one container
n_bpp = len(w_bpp)
def modello_bpp(w, c, k):
n = len(w)
m = nuovo_modello("bin_packing")
x = m.addVars(n, k, vtype=GRB.BINARY, name="x")
y = m.addVars(k, vtype=GRB.BINARY, name="y")
m.setObjective(y.sum(), GRB.MINIMIZE)
m.addConstrs((x.sum(j, "*") == 1 for j in R(n)), name="item")
m.addConstrs((gp.quicksum(w[j] * x[j, b] for j in R(n)) <= c * y[b] for b in R(k)),
name="capacity")
return m, x, y
# with one container per item one finds how many are really needed; the model that
# gets printed then uses only those, because the others would stay empty
m_largo, _, _ = modello_bpp(w_bpp, c_bpp, n_bpp)
z_bpp = risolvi(m_largo)
minimo_teorico = -(-sum(w_bpp) // c_bpp) # rounding up
print(f" Weights {w_bpp}, capacity {c_bpp}.")
print(f" The total weight is {sum(w_bpp)}: no solution uses fewer than "
f"{sum(w_bpp)}/{c_bpp} = {minimo_teorico} containers, and the optimum uses {int(z_bpp)}.")
assert z_bpp == minimo_teorico
m_bpp, x_bpp, y_bpp = modello_bpp(w_bpp, c_bpp, int(z_bpp))
risolvi(m_bpp)
salva_modello(m_bpp, "cap06_bpp")
print(" The relaxation buys fractions of a container and drops to "
f"{frazione(sum(w_bpp) / c_bpp)}: it does not know a container opens whole.")
salva_dati(pd.DataFrame([{"problem": "bin packing", "z_milp": z_bpp}]), "cap06_bpp")
print("Done.")