3.4 Tre modelli classici
Classe: implementazione · Script: python/cap06_bpp.py, python/cap06_cmax.py, python/cap06_tsp.py
Bin packing, makespan e commesso viaggiatore: enunciato, modello, costruzione in
gurobipy e modello dell'istanza. Sono i tre problemi su cui il
capitolo delle euristiche costruisce next-fit, first-fit,
best-fit, LPT e vicino più vicino.
Fin qui il modello di esempio è sempre stato lo zaino. I tre problemi qui sotto ritornano nel capitolo delle euristiche, dove si costruiscono a mano le soluzioni di next-fit, first-fit, best-fit, LPT e vicino più vicino: qui si scrivono i loro modelli, così quelle euristiche hanno un ottimo con cui confrontarsi.
Bin packing: quanti contenitori bastano
Lo script è python/cap06_bpp.py.
Bin packing
Ci sono \(n\) oggetti, l'oggetto \(j\) pesa \(w_j\). I contenitori sono tutti uguali, di capacità \(c\). Si usi il minimo numero di contenitori.
Il modello
Servono due famiglie di variabili binarie: \(x_{jb} = 1\) se l'oggetto \(j\) va nel contenitore \(b\), e \(y_b = 1\) se il contenitore \(b\) viene usato.
La prima famiglia dice che ogni oggetto finisce in esattamente un contenitore. La seconda è la capacità scritta come attivazione: finché \(y_b = 0\) il contenitore \(b\) non può ricevere niente, e appena \(y_b = 1\) accoglie fino a \(c\). L'obiettivo conta i contenitori accesi.
La costruzione 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="oggetto")
m.addConstrs((gp.quicksum(w[j] * x[j, b] for j in R(n)) <= c * y[b]
for b in R(k)), name="capacita")
return m, x, y
L'istanza
Sull'istanza di quattro oggetti di peso \(w = (5, 4, 3, 3)\) e capacità \(c = 7\):
Il peso totale è \(15\), quindi nessuna soluzione può usare meno di \(\lceil 15/7 \rceil = 3\) contenitori; l'ottimo ne usa esattamente \(3\), e il conteggio è quindi stretto.
Il rilassamento del bin packing è debolissimo
Rilassando \(y_b\) a \(y_b \ge 0\) il modello compra frazioni di contenitore, e l'ottimo dell'LP scende a \(\sum_j w_j / c = 15/7 \approx 2{,}14\): il rilassamento non sa che un contenitore si apre tutto intero. È il motivo per cui su questo problema il bound duale è poco utile e le euristiche contano di più.
Makespan: il carico della macchina più carica
Lo script è python/cap06_cmax.py.
Makespan su macchine identiche
Ci sono \(n\) lavori, di durata \(d_j\), e \(k\) macchine identiche. Ogni lavoro va su una macchina sola e non si interrompe. Si minimizzi l'istante in cui l'ultima macchina finisce.
Il modello
Con \(x_{jm} = 1\) se il lavoro \(j\) va sulla macchina \(m\), e \(z \ge 0\) il carico della macchina più carica:
L'obiettivo è la sola variabile \(z\): nessun dato vi compare. Sono le \(k\) righe di carico a darle significato, dicendo che nessuna macchina lavora più a lungo di \(z\); il minimo la schiaccia allora sul carico della macchina più carica. È la tecnica min-max.
La costruzione 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="lavoro")
m.addConstrs((gp.quicksum(d[j] * x[j, mm] for j in R(n)) <= z
for mm in R(k)), name="carico")
return m, x, z
L'istanza
Sull'istanza di quattro lavori di durata \(d = (3, 4, 5, 6)\) su \(k = 2\) macchine --- la stessa su cui il capitolo delle euristiche fa correre LPT:
Il carico totale è \(18\) e le macchine sono due: nessuna soluzione può scendere sotto \(18/2 = 9\), e l'ottimo vale esattamente \(9\) — i lavori si dividono in \(6+3\) e \(5+4\). Qui il conteggio chiude il problema da solo.
Commesso viaggiatore: la formulazione MTZ
Lo script è python/cap06_tsp.py.
Commesso viaggiatore
Ci sono \(n\) città e una distanza \(d_{ij}\) fra ogni coppia. Si trovi il giro di lunghezza minima che tocca ogni città esattamente una volta e torna al punto di partenza.
Il modello
Con \(x_{ij} = 1\) se il giro va da \(i\) a \(j\), le due famiglie «si esce una volta» e «si entra una volta» non bastano: ammettono anche soluzioni fatte di sottocicli separati. La formulazione di Miller–Tucker–Zemlin aggiunge una variabile \(u_i\) per ogni città diversa dalla prima, che ne registra la posizione lungo il giro.
Il terzo gruppo è il cuore della formulazione. Se \(x_{ij} = 0\) la riga diventa \(u_i - u_j \le n - 1\), sempre vera perché le \(u\) stanno fra \(1\) e \(n-1\): non vieta niente. Se invece \(x_{ij} = 1\) diventa \(u_j \ge u_i + 1\), cioè «se vado da \(i\) a \(j\), la posizione di \(j\) è la successiva». Un sottociclo che non tocca la città \(1\) richiederebbe una catena di posizioni sempre crescenti che si richiude su se stessa, e questo è impossibile; la città \(1\) non ha la sua \(u\) proprio perché è il punto in cui il giro si chiude.
La costruzione 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="esce")
m.addConstrs((gp.quicksum(x[i, j] for i in R(n) if i != j) == 1
for j in R(n)), name="entra")
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
L'istanza
Sull'istanza di quattro città del capitolo delle euristiche il giro ottimo è \(1 \to 2 \to 4 \to 3 \to 1\) e misura \(22\). Il modello dell'istanza ha quindici colonne — dodici archi e tre posizioni — e quattordici righe:
MTZ è comoda, non è la più forte
I vincoli MTZ sono \(O(n^2)\) e si scrivono in tre righe di gurobipy, ma il
loro rilassamento lineare è debole: le \(u\) continue assorbono quasi tutto e
l'LP si avvicina poco all'ottimo intero. Le formulazioni che eliminano i
sottocicli con i tagli di connessione danno bound molto migliori, al prezzo
di un numero esponenziale di vincoli da generare a mano a mano. Per le
dimensioni di questo corso MTZ basta.
Mostra lo script completo — python/cap06_bpp.py (51 righe)
"""Bin packing: quanti contenitori bastano (capitolo 3).
Il primo dei tre problemi che il capitolo delle euristiche riprende: li' si
costruiscono a mano next-fit, first-fit e best-fit, qui si scrive il modello che
dice qual e' l'ottimo. La capacita' e' scritta come attivazione: finche' il
contenitore non si apre non puo' ricevere niente.
"""
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: il minimo numero di contenitori")
w_bpp = [5, 4, 3, 3] # peso degli oggetti
c_bpp = 7 # capacita' di un contenitore
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="oggetto")
m.addConstrs((gp.quicksum(w[j] * x[j, b] for j in R(n)) <= c * y[b] for b in R(k)),
name="capacita")
return m, x, y
# con un contenitore per oggetto si trova quanti ne servono davvero; il modello
# che si stampa usa poi solo quelli, perche' gli altri resterebbero vuoti
m_largo, _, _ = modello_bpp(w_bpp, c_bpp, n_bpp)
z_bpp = risolvi(m_largo)
minimo_teorico = -(-sum(w_bpp) // c_bpp) # arrotondamento all'insu'
print(f" Pesi {w_bpp}, capacita' {c_bpp}.")
print(f" Il peso totale e' {sum(w_bpp)}: nessuna soluzione usa meno di "
f"{sum(w_bpp)}/{c_bpp} = {minimo_teorico} contenitori, e l'ottimo ne usa {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(" Il rilassamento compra frazioni di contenitore e scende a "
f"{frazione(sum(w_bpp) / c_bpp)}: non sa che un contenitore si apre tutto intero.")
salva_dati(pd.DataFrame([{"problema": "bin packing", "z_milp": z_bpp}]), "cap06_bpp")
print("Fine.")