Una macchina, classi di lavori con setup
Classe: BIP · Legami: attivazione disaggregata, CNF · Script: python/fam07_5_classisetup.py
Difficoltà: ★★★☆☆ · Tempo: 30–45 min
Problema 7.5
Un'azienda ha \(n\) lavori eseguibili su una macchina con disponibilità \(a \in \mathbb{Q}_{>0}\) minuti. Per ogni lavoro \(j\), \(t_j\) è il tempo e \(r_j\) il ricavo se eseguito. I lavori sono partizionati in \(q \ge 2\) classi \(\mathscr{J}_1, \dots, \mathscr{J}_q\). Se la macchina esegue lavori di una classe \(c\), si paga un costo di setup \(f_c \ge 0\) e si consuma un tempo di setup \(s_c \ge 0\). La macchina non esegue lavori in parallelo. Massimizzare il profitto.
Il problema a parole. Decidiamo quali lavori eseguire e quali classi attivare. L'obiettivo: ricavi meno costi di setup. Il vincolo: tempi dei lavori più tempi di setup entro la disponibilità. Uno zaino con costi fissi per gruppo: il legame di attivazione, stavolta disaggregato fin dall'inizio.
Modello
| Simbolo | Tipo | Significato |
|---|---|---|
| \(n\), \(a\) | numero di lavori, disponibilità | |
| \(t_j\), \(r_j\) | \(\in \mathbb{Q}_{>0}\) | tempo e ricavo del lavoro \(j\) |
| \(q\), \(\mathscr{J}_c\) | numero di classi e lavori della classe \(c\) (partizione) | |
| \(f_c\), \(s_c\) | \(\in \mathbb{Q}_{\ge 0}\) | costo e tempo di setup della classe \(c\) |
Variabili. \(n + q\) binarie: \(x_j = 1\) se il lavoro \(j\) è eseguito; \(y_c = 1\) se almeno un lavoro della classe \(c\) è eseguito.
- l'obiettivo massimizza ricavi meno setup;
- il vincolo di disponibilità (\(1\) vincolo lineare);
- i vincoli di link: se un lavoro di una classe è eseguito, la classe è attivata (\(n\) vincoli lineari, uno per lavoro);
- i vincoli di dominio definiscono le variabili.
Legame fra le variabili: la CNF diventa vincolo
Dal vincolo. «Se almeno un lavoro della classe \(c\) è eseguito, la classe è attivata»: \((x_j \,\mathtt{OR}\, x_s \,\mathtt{OR}\, \dots) \Rightarrow y_c\), contronominale \(\mathtt{NOT}\,y_c \Rightarrow (\mathtt{NOT}\,x_j \,\mathtt{AND}\, \dots)\). L'espressione \(\mathtt{NOT}(x_j \,\mathtt{OR}\, \dots) \,\mathtt{OR}\, y_c\) diventa, con De Morgan e distributività, la CNF \((\mathtt{NOT}\,x_j \,\mathtt{OR}\, y_c) \,\mathtt{AND}\, (\mathtt{NOT}\,x_s \,\mathtt{OR}\, y_c) \,\mathtt{AND}\, \dots\), cioè \(1 - x_j + y_c \ge 1\): esattamente i vincoli di link \(x_j \le y_c\). Verifica nei due versi: \(x_j = 1\) forza \(y_c = 1\); \(y_c = 0\) forza tutti gli \(x_j\) della classe a \(0\).
Dall'ottimo. «Se nessun lavoro della classe è eseguito, la classe non è attivata»: non imposta, segue senza perdita di ottimalità: porre \(y_c = 0\) resta ammissibile, libera \(s_c\) minuti e non diminuisce l'obiettivo perché \(f_c \ge 0\). Poiché \(f_c\) può essere nullo, la conclusione corretta è «esiste un ottimo in cui…», non «in ogni ottimo».
Il modello in gurobipy
m = gp.Model("classi_setup"); m.Params.OutputFlag = 0
x = m.addVars(n, vtype=GRB.BINARY, name="x")
y = m.addVars(q, vtype=GRB.BINARY, name="y")
m.setObjective(gp.quicksum(r[j] * x[j] for j in range(n))
- gp.quicksum(f[c] * y[c] for c in range(q)), GRB.MAXIMIZE)
m.addConstr(gp.quicksum(t[j] * x[j] for j in range(n))
+ gp.quicksum(s[c] * y[c] for c in range(q)) <= a, name="disponibilita")
m.addConstrs((x[j] - y[c] <= 0 for c in range(q) for j in J[c]), name="link")
m.optimize()
L'istanza
\(n = 7\), \(q = 3\): \(\mathscr{J}_1 = \{1, 2\}\), \(\mathscr{J}_2 = \{3, 4\}\), \(\mathscr{J}_3 = \{5, 6, 7\}\), \(a = 50\).
| \(j=1\) | \(j=2\) | \(j=3\) | \(j=4\) | \(j=5\) | \(j=6\) | \(j=7\) | |
|---|---|---|---|---|---|---|---|
| \(r_j\) | 10 | 6 | 8 | 6 | 7 | 9 | 5 |
| \(t_j\) | 5 | 10 | 8 | 6 | 9 | 5 | 6 |
| \(c=1\) | \(c=2\) | \(c=3\) | |
|---|---|---|---|
| \(f_c\) | 10 | 5 | 4 |
| \(s_c\) | 10 | 12 | 6 |
Il modello scritto sui dati dell'istanza:
Euristica costruttiva: il bound primale
Classe per classe: il primo lavoro paga anche il setup, se ci sta.
- Passo 1. Classe 1: \(s_1 + t_1 = 15 \le 50\); \(y[1] = x[1] = 1\), \(ra = 35\).
- Passo 2. \(t_2 = 10 \le 35\); \(x[2] = 1\), \(ra = 25\).
- Passo 3. Classe 2: \(s_2 + t_3 = 20 \le 25\); \(y[2] = x[3] = 1\), \(ra = 5\).
- Passo 4. \(t_4 = 6 > 5\): saltato. Passi 5–7. Classe 3: \(s_3 + t_j > 5\): saltati.
Profitto \(10 + 6 + 8 - 10 - 5 = 9\): \(z(\mathit{MILP}) \ge 9\).
Rilassamento LP e duale: il bound duale
Con \(\pi \ge 0\) (disponibilità) e \(\lambda_j \ge 0\) (link):
Lo stesso duale, scritto sui dati dell'istanza:
Una soluzione duale a mano. \(\bar\lambda = 0\) e \(\bar\pi = \max_j r_j/t_j = \tfrac{10}{5} = 2\): valore \(100\). Quindi \(9 \le z(\mathit{MILP}) \le 100\): un bound grossolano, come spesso i bound «di zaino», che ignora setup e costi.
Quello che dice il solver. \(z(\mathit{LP}) = 425/13 = 32{,}7\) (con \(\pi = \tfrac{17}{26}\) e alcuni \(\lambda_j > 0\)); \(z(\mathit{LP}^+) = 329/13\). Ottimo intero \(21\): classi 2 e 3, lavori \(3, 4, 5, 6\), profitto \(30 - 9\). L'euristica resta a \(9\) (gap \(57\%\)): l'ordine di scansione conta.
| \(LB\) | \(UB\) (duale a mano) | \(z(\mathit{LP})\) | \(z(\mathit{LP}^+)\) | \(z(\mathit{MILP})\) | gap euristica |
|---|---|---|---|---|---|
| 9 | 100 | \(425/13\) | \(329/13\) | 21 | \(57{,}1\%\) |
Considerazioni aggiuntive
- \(y_c \le 1\) rafforza il rilassamento; \(x_j \le 1\) è implicato.
- \(\sum_{j \in \mathscr{J}_c} x_j \ge y_c\) (\(q\) vincoli) non è valido ma preserva l'ottimo.
- La forma aggregata \(\sum_{j \in \mathscr{J}_c} x_j \le |\mathscr{J}_c|\, y_c\) ha lo stesso insieme intero e un rilassamento più debole.
Domande di modellazione aggiuntive
7.5.1 — Una classe subordinata a un'altra
La classe 3 si può attivare solo se si attiva anche la classe 1.
Una variante svolta: una sola classe
Cambiare classe richiede una pulizia della macchina che l'azienda vuole evitare del tutto: si può attivare al più una classe.
Un vincolo di set packing sulle attivazioni, \(\sum_{c=1}^{q} y_c \le 1\) (un vincolo lineare): al più una \(y_c\) vale \(1\), e con essa, tramite, solo i lavori di quella classe possono essere eseguiti. Sull'istanza l'ottimo scende a \(17\): la classe migliore da sola è la 3, con tutti i suoi lavori (\(9 + 5 + 6 + 6 = 26 \le 50\)), profitto \(21 - 4 = 17\).
«Al più una classe attivata» aggiunge al duale una variabile \(\theta \ge 0\) con termine noto \(1\): entra nell'obiettivo e allenta le colonne delle \(y_c\). Fissato il prezzo del tempo \(\pi\), le \(\lambda_j\) e \(\theta\) sono forzate dai vincoli duali, quindi la ricetta si riduce a cercare \(\pi\) fra i rapporti \(r_j/t_j\) e tenere il valore più basso. Il bound migliora di molto rispetto alla ricetta del problema base.
| valore | che cos'è | |
|---|---|---|
| \(\mathit{UB}\) | \(\frac{157}{5}\) | certificato duale costruito a mano |
| \(\mathit{LB}\) | \(17\) | soluzione euristica |
| \(z(\mathit{LP})\) | \(17\) | rilassamento senza i bound |
| \(z(\mathit{LP}^+)\) | \(17\) | rilassamento con i bound |
| \(z(\mathit{MILP})\) | \(17\) | ottimo del MILP |
Codice
Script completo: python/fam07_5_classisetup.py;
notebook: notebooks/fam07_5_classisetup.ipynb.
Mostra lo script completo — python/fam07_5_classisetup.py (214 righe)
"""Problema 7.5 -- Una macchina, classi di lavori con setup.
Il legame di attivazione disaggregato dedotto passo passo dalla CNF di
un'implicazione booleana: (OR di lavori) => classe attivata.
"""
import gurobipy as gp
import numpy as np
import pandas as pd
from gurobipy import GRB
from euristiche import best_fit, first_fit, matrice, next_fit
from mip import (ammissibile, dualita_forte, due_rilassamenti, frazione,
nuovo_modello, registra_bound, rilassamenti, risolvi,
stampa_soluzione, valuta)
from stile import CICLO, ROSSO, intestazione, plt, salva_dati, salva_figura
from esteso import salva_modello
R = range
# ---------- 1. MODELLO E ISTANZA ----------
intestazione("5. Classi di lavori con costo e tempo di setup: y_c attiva la classe")
r5 = [10, 6, 8, 6, 7, 9, 5]
t5 = [5, 10, 8, 6, 9, 5, 6]
J5 = [[0, 1], [2, 3], [4, 5, 6]] # classi (0-based)
f5 = [10, 5, 4]
s5 = [10, 12, 6]
a5 = 50
salva_dati(pd.DataFrame({"lavoro": R(1, 8), "r": r5, "t": t5,
"classe": [c + 1 for j in R(7) for c in R(3) if j in J5[c]]}), "fam07_5_lavori")
salva_dati(pd.DataFrame({"classe": R(1, 4), "f": f5, "s": s5}), "fam07_5_classi")
def modello_5(r, t, J, f, s, a):
n, q = len(r), len(J)
m = nuovo_modello("classi_setup")
x = m.addVars(n, vtype=GRB.BINARY, name="x")
y = m.addVars(q, vtype=GRB.BINARY, name="y")
m.setObjective(gp.quicksum(r[j] * x[j] for j in R(n)) - gp.quicksum(f[c] * y[c] for c in R(q)),
GRB.MAXIMIZE)
m.addConstr(gp.quicksum(t[j] * x[j] for j in R(n)) + gp.quicksum(s[c] * y[c] for c in R(q)) <= a,
name="disponibilita")
m.addConstrs((x[j] - y[c] <= 0 for c in R(q) for j in J[c]), name="link")
return m, x, y
def duale_5(r, t, J, f, s, a):
"""min a pi; t_j pi + lam_j >= r_j; s_c pi - sum_{j in J_c} lam_j >= -f_c; pi, lam >= 0."""
n, q = len(r), len(J)
d = nuovo_modello("duale_classi_setup")
pi = d.addVar(name="pi")
lam = d.addVars(n, name="lam")
d.setObjective(a * pi, GRB.MINIMIZE)
d.addConstrs((t[j] * pi + lam[j] >= r[j] for j in R(n)), name="rc_x")
d.addConstrs((s[c] * pi - gp.quicksum(lam[j] for j in J[c]) >= -f[c] for c in R(q)), name="rc_y")
return d
def euristica_5(r, t, J, f, s, a):
"""Classe per classe: il primo lavoro paga anche il setup, se ci sta."""
n, q = len(r), len(J)
x, y, ra, passi = [0] * n, [0] * q, a, []
for c in R(q):
for j in J[c]:
if y[c] == 0:
if s[c] + t[j] <= ra:
y[c], x[j] = 1, 1
passi.append(f"Classe {c + 1} non attiva: s[{c + 1}] + t[{j + 1}] = {s[c]} + {t[j]} = "
f"{s[c] + t[j]} <= ra = {ra}; y[{c + 1}] = 1, x[{j + 1}] = 1, ra = {ra - s[c] - t[j]}.")
ra -= s[c] + t[j]
else:
passi.append(f"Classe {c + 1} non attiva: s[{c + 1}] + t[{j + 1}] = {s[c] + t[j]} > ra = {ra}; "
f"il lavoro {j + 1} viene saltato.")
else:
if t[j] <= ra:
x[j] = 1
passi.append(f"Classe {c + 1} attiva: t[{j + 1}] = {t[j]} <= ra = {ra}; x[{j + 1}] = 1, ra = {ra - t[j]}.")
ra -= t[j]
else:
passi.append(f"Classe {c + 1} attiva: t[{j + 1}] = {t[j]} > ra = {ra}; il lavoro {j + 1} viene saltato.")
return x, y, passi
m5, x5, y5 = modello_5(r5, t5, J5, f5, s5, a5)
salva_modello(m5, "fam07_5_primale")
# ---------- 2. IL RILASSAMENTO LP ----------
zlp5, zlp5r, _ = rilassamenti(m5)
# ---------- 3. IL DUALE DEL RILASSAMENTO (UPPER BOUND: E' UN MASSIMO) ----------
d5 = duale_5(r5, t5, J5, f5, s5, a5)
salva_modello(d5, "fam07_5_duale")
pi_mano = max(r5[j] / t5[j] for j in R(7))
ub5, viol = valuta(d5, {"pi": pi_mano})
assert viol <= 1e-9
print(f"Soluzione duale a mano: lam = 0, pi = max_j r_j/t_j = {frazione(pi_mano)} -> ub = {frazione(ub5)}")
dualita_forte(d5, zlp5)
# ---------- 4. EURISTICA COSTRUTTIVA (LOWER BOUND: E' UN MASSIMO) ----------
xe, ye, passi = euristica_5(r5, t5, J5, f5, s5, a5)
print("Euristica classe per classe:")
for i, s in enumerate(passi, 1):
print(f" Passo {i}. {s}")
lb5 = sum(r5[j] * xe[j] for j in R(7)) - sum(f5[c] * ye[c] for c in R(3))
print(f" lb = {lb5} (x = {xe}, y = {ye})")
# ---------- 5. SOLUZIONE OTTIMA DEL MILP ----------
z5 = risolvi(m5)
print("Soluzione ottima del MILP:")
stampa_soluzione(m5, solo_non_nulle=True)
riga = registra_bound("5 classi setup", ub5, lb5, zlp5, zlp5r, z5, senso="max")
salva_dati(pd.DataFrame([riga]), "fam07_5_bound")
# ---------- 6. DOMANDE DI MODELLAZIONE AGGIUNTIVE ----------
varianti = {}
def variante(nome, m):
z = risolvi(m)
print(f" {nome:70s} z = {frazione(z)}")
return z
# 5a: una sola classe attiva
m, x, y = modello_5(r5, t5, J5, f5, s5, a5)
m.addConstr(y.sum() <= 1, name="una_classe")
varianti["5a"] = variante("5a. Al più una classe attivata (sum y_c <= 1)", m)
# 5b: la classe 3 solo se la classe 1
m, x, y = modello_5(r5, t5, J5, f5, s5, a5)
m.addConstr(y[2] <= y[0], name="3_solo_se_1")
varianti["5b"] = variante("5b. La classe 3 si attiva solo se si attiva la classe 1 (y_3 <= y_1)", m)
salva_dati(pd.DataFrame({"variante": list(varianti), "z": list(varianti.values())}), "fam07_5_varianti")
# ---------- 7. IL SANDWICH SULLA VARIANTE 5a ----------
intestazione("5a. Il sandwich sulla variante: al piu' una classe attivata")
def modello_5a(r, t, J, f, s, a):
mm_, xx, yy = modello_5(r, t, J, f, s, a)
mm_.addConstr(yy.sum() <= 1, name="una_classe")
return mm_, xx, yy
def duale_5a(r, t, J, f, s, a):
"""Al duale di 7.5 si aggiunge theta >= 0 per il vincolo sum_c y_c <= 1.
Il termine noto e' 1, quindi theta entra nell'obiettivo; e compare con
segno piu' nelle colonne delle y_c, che cosi' si allentano."""
nn, q = len(r), len(J)
d = nuovo_modello("duale_classi_setup_5a")
pi = d.addVar(name="pi")
lam = d.addVars(nn, name="lam")
th = d.addVar(name="theta")
d.setObjective(a * pi + th, GRB.MINIMIZE)
d.addConstrs((t[j] * pi + lam[j] >= r[j] for j in R(nn)), name="rc_x")
d.addConstrs((s[c] * pi - gp.quicksum(lam[j] for j in J[c]) + th >= -f[c] for c in R(q)),
name="rc_y")
return d
m5a, x5a, y5a = modello_5a(r5, t5, J5, f5, s5, a5)
salva_modello(m5a, "fam07_5a_primale")
# -- euristica ammissibile: la base, ristretta alla classe migliore --
print("Euristica costruttiva: una classe per volta, si tiene la migliore. Dentro la classe")
print("i lavori entrano per rapporto r_j/t_j decrescente finche' il tempo lo permette.")
migliore_5a = (0.0, None, [])
for c in R(len(J5)):
residuo = a5 - s5[c]
presi = []
for j in sorted(J5[c], key=lambda j: -r5[j] / t5[j]):
if t5[j] <= residuo:
presi.append(j); residuo -= t5[j]
valore = sum(r5[j] for j in presi) - f5[c]
print(f" classe {c + 1}: setup {s5[c]} minuti e costo {f5[c]}; lavori "
f"{[j + 1 for j in sorted(presi)]} -> {frazione(valore)}")
if valore > migliore_5a[0]:
migliore_5a = (valore, c, presi)
lb5a, classe_5a, presi_5a = migliore_5a
sol_5a = {f"x[{j}]": 1 for j in presi_5a} | {f"y[{classe_5a}]": 1}
assert ammissibile(m5a, sol_5a), "la soluzione euristica della variante deve essere ammissibile"
print(f" la migliore e' la classe {classe_5a + 1} -> lb = {frazione(lb5a)}")
# -- certificato duale: il prezzo del tempo si cerca fra i rapporti r_j/t_j --
d5a = duale_5a(r5, t5, J5, f5, s5, a5)
salva_modello(d5a, "fam07_5a_duale")
def valore_duale_5a(pi_val):
"""Dato il prezzo del tempo, le altre variabili duali sono forzate."""
lam_v = {j: max(0.0, r5[j] - t5[j] * pi_val) for j in R(len(r5))}
th_v = max([0.0] + [sum(lam_v[j] for j in J5[c]) - s5[c] * pi_val - f5[c]
for c in R(len(J5))])
return a5 * pi_val + th_v, lam_v, th_v
candidati = sorted({r5[j] / t5[j] for j in R(len(r5))})
scelto = min(candidati, key=lambda p: valore_duale_5a(p)[0])
ub5a, lam_5a, th_5a = valore_duale_5a(scelto)
mano_5a = {"pi": scelto, "theta": th_5a} | {f"lam[{j}]": lam_5a[j] for j in R(len(r5))}
ub5a_val, viol_5a = valuta(d5a, mano_5a)
assert viol_5a <= 1e-9, viol_5a
print("Soluzione duale a mano: il prezzo del tempo pi si cerca fra i rapporti r_j/t_j;")
print(" fissato pi, le lam_j e theta sono forzate dai vincoli duali. Si tiene il pi che")
print(f" da' il valore piu' basso: pi = {frazione(scelto)}, theta = {frazione(th_5a)}")
print(f" -> ub = {frazione(ub5a_val)} (la ricetta del problema base, pi = max_j r_j/t_j,")
print(f" darebbe {frazione(a5 * max(r5[j] / t5[j] for j in R(len(r5))))})")
zlp5a, zlp5ar, _ = due_rilassamenti(m5a, d5a)
z5a = risolvi(m5a)
riga_5a = registra_bound("5a al piu' una classe", ub5a_val, lb5a, zlp5a, zlp5ar, z5a, senso="max")
salva_dati(pd.DataFrame([riga_5a]), "fam07_5a_bound")
assert lb5a <= z5a <= zlp5a + 1e-9 <= ub5a_val + 1e-9
print("Fine.")