Localizzazione capacitata
Classe: MILP · Legami: attivazione aggregata (anche vincolo di capacità) · Script: python/fam08_1_capacitata.py
Difficoltà: ★★★☆☆ · Tempo: 30–45 min
Problema 8.1
Un'azienda deve servire \(n \in \mathbb{Z}_{\ge 1}\) clienti e ha individuato \(m \in \mathbb{Z}_{\ge 1}\) sedi candidate. Per ogni cliente \(c\), \(d_c \in \mathbb{Q}_{>0}\) è la domanda in litri. Per ogni sede \(l\) e cliente \(c\), \(t_{lc} \in \mathbb{Q}_{>0}\) è il costo di trasporto per litro. Per ogni sede \(l\), \(u_l \in \mathbb{Q}_{>0}\) è la capacità e \(i_l \in \mathbb{Q}_{>0}\) il costo di installazione. Si vuole decidere dove installare e come servire i clienti, a costo minimo.
Il problema a parole. Decidiamo dove installare le strutture e quanto spedire da ciascuna sede a ciascun cliente. L'obiettivo: costo totale (installazione più trasporto) minimo. I vincoli: da una sede non installata non parte nulla, e una installata non supera la capacità; la domanda va soddisfatta esattamente. È la localizzazione capacitata.
Modello
Dati.
| Simbolo | Tipo | Significato |
|---|---|---|
| \(m\) | \(\in \mathbb{Z}_{\ge 1}\) | numero di sedi, \(l \in \{1, 2, \dots, m\}\) |
| \(n\) | \(\in \mathbb{Z}_{\ge 1}\) | numero di clienti, \(c \in \{1, 2, \dots, n\}\) |
| \(t_{lc}\) | \(\in \mathbb{Q}_{>0}\) | costo di trasporto dalla sede \(l\) al cliente \(c\) |
| \(u_l\) | \(\in \mathbb{Q}_{>0}\) | capacità della sede \(l\) |
| \(i_l\) | \(\in \mathbb{Q}_{>0}\) | costo di installazione della sede \(l\) |
| \(d_c\) | \(\in \mathbb{Q}_{>0}\) | domanda del cliente \(c\) |
Variabili decisionali. \(m\) binarie \(x_l\) (sede \(l\) installata) e \(m\,n\) continue non negative \(y_{lc}\) (litri spediti da \(l\) a \(c\)):
Modello MILP:
- l'obiettivo minimizza il costo totale (installazione più trasporto);
- il primo vincolo lega trasporto e installazione e impone la capacità (\(m\) vincoli lineari);
- il secondo soddisfa la domanda di ogni cliente (\(n\) vincoli lineari);
- i vincoli restanti definiscono le variabili.
Il legame. Se una quantità positiva parte dalla sede \(l\), la sede deve essere installata; dalla contronominale, una sede chiusa non spedisce nulla. Entrambi i versi sono imposti direttamente dal primo vincolo. Il verso opposto — una sede installata spedisce qualcosa — non è imposto ma segue dall'obiettivo: poiché \(i_l > 0\), un ottimo non lascia mai una sede aperta inutilizzata. Una sola famiglia di vincoli fa dunque sia da legame di attivazione sia da vincolo di capacità.
Il modello in gurobipy
mod = gp.Model("localizzazione_capacitata")
x = mod.addVars(m, vtype=GRB.BINARY, name="x")
y = mod.addVars(m, n, name="y")
mod.setObjective(gp.quicksum(i[l] * x[l] for l in range(m))
+ gp.quicksum(t[l][c] * y[l, c] for l in range(m) for c in range(n)), GRB.MINIMIZE)
mod.addConstrs((u[l] * x[l] - gp.quicksum(y[l, c] for c in range(n)) >= 0
for l in range(m)), name="capacita")
mod.addConstrs((gp.quicksum(y[l, c] for l in range(m)) == d[c] for c in range(n)), name="domanda")
L'istanza
\(m = 2\) sedi, \(n = 3\) clienti:
| \(t_{lc}\) | \(c=1\) | \(c=2\) | \(c=3\) |
|---|---|---|---|
| \(l=1\) | 4 | 5 | 6 |
| \(l=2\) | 6 | 4 | 3 |
| \(l=1\) | \(l=2\) | |
|---|---|---|
| \(u_l\) | 50 | 50 |
| \(i_l\) | 60 | 90 |
| \(c=1\) | \(c=2\) | \(c=3\) | |
|---|---|---|---|
| \(d_c\) | 8 | 25 | 27 |
Il modello scritto sui dati dell'istanza:
Euristica costruttiva: il bound primale
Si scandiscono le sedi in ordine; per ciascuna, i clienti, spedendo il minimo fra capacità residua e domanda residua.
Esecuzione: la sede 1 spedisce \(8\) al cliente 1, \(25\) al cliente 2, \(17\) al cliente 3 (capacità esaurita); la sede 2 spedisce i restanti \(10\) al cliente 3. Valore: \(60+90 + (4{\cdot}8+5{\cdot}25+6{\cdot}17+3{\cdot}10) = 150+289 = 439\). Quindi \(z(\mathit{MILP}) \le \mathit{UB} = 439\).
Rilassamento LP e duale: il bound duale
Il duale del rilassamento lineare, una variabile per vincolo del primale:
Lo stesso duale, scritto sui dati dell'istanza:
Con \(\bar\mu_l = i_l/u_l\) (spalma il costo fisso sulla capacità) e \(\bar\pi_c = \min_l(t_{lc}+\bar\mu_l)\):
di valore \(8{\cdot}26/5 + 25{\cdot}29/5 + 27{\cdot}24/5 = 1581/5\). Per la dualità debole, \(\mathit{LB} = 1581/5 \le z(\mathit{LP}) \le z(\mathit{MILP}) \le \mathit{UB} = 439\).
Quello che dice il solver. \(z(\mathit{LP}) = 1581/5\) esattamente: la soluzione duale a mano è già ottima. Rafforzando con \(x_l \le 1\), \(z(\mathit{LP}^+) = 317\). \(z(\mathit{MILP}) = 365\), con entrambe le sedi aperte: la sede 1 serve il cliente 1 e parte del cliente 2, la sede 2 il resto del cliente 2 e tutto il cliente 3. Gap euristica \(20{,}3\%\).
| \(UB\) | \(LB\) (duale) | \(z(\mathit{LP})\) | \(z(\mathit{LP}^+)\) | \(z(\mathit{MILP})\) | gap dell'euristica |
|---|---|---|---|---|---|
| 439 | \(1581/5\) | \(1581/5\) | 317 | 365 | \(20{,}3\%\) |
Considerazioni aggiuntive
- Se \(u_l < d_c\) nessuna sede da sola può soddisfare il cliente \(c\): non il caso qui, ma va verificato.
- \(y_{lc} \le d_c\, x_l\) è valida ma implicata dai due vincoli insieme.
Domande di modellazione aggiuntive
8.1.1 — Lotto minimo per ogni sede aperta
Ogni sede aperta deve spedire almeno \(5\) litri. Come cambia il modello? Qual è il nuovo ottimo?
Una variante svolta: apertura condizionata
La sede 2 può essere installata solo se lo è anche la sede 1 (ad esempio, un vincolo logistico di supervisione).
È un legame fra due variabili della stessa famiglia, imposto dal singolo vincolo lineare
(un vincolo lineare): se \(x_2 = 1\) allora necessariamente \(x_1 = 1\). Sull'istanza l'ottimo del problema 8.1 apre già entrambe le sedi, quindi il vincolo aggiuntivo non è vincolante e l'ottimo resta \(365\).
«La sede 2 si apre solo se si apre la sede 1» aggiunge al duale \(\rho \le 0\), che sposta costo fisso fra le due sedi: \(\mu_1 = (i_1 + \rho)/u_1\) e \(\mu_2 = (i_2 - \rho)/u_2\). Si prova \(\rho\) fra i valori che fanno cambiare il minimo che definisce \(\pi_c\); su questa istanza il migliore è \(\rho = 0\), perché la sede 1 ha già il costo per litro più basso.
| valore | che cos'è | |
|---|---|---|
| \(\mathit{UB}\) | \(439\) | soluzione euristica |
| \(\mathit{LB}\) | \(\frac{1581}{5}\) | certificato duale costruito a mano |
| \(z(\mathit{LP})\) | \(325\) | rilassamento senza i bound |
| \(z(\mathit{LP}^+)\) | \(325\) | rilassamento con i bound |
| \(z(\mathit{MILP})\) | \(365\) | ottimo del MILP |
Codice
Script completo —
python/fam08_1_capacitata.py
(riproducibile con python3 python/fam08_1_capacitata.py dalla cartella
python/). Notebook —
notebooks/fam08_1_capacitata.ipynb
— che si apre in Colab dal badge in cima alla pagina.
Mostra lo script completo — python/fam08_1_capacitata.py (247 righe)
"""Problema 8.1 -- Localizzazione capacitata (costo minimo).
Attivazione aggregata fra la variabile binaria x_l (apri la sede l) e le
variabili continue di flusso y_lc: il legame si dimostra nei due versi
esattamente come nel problema 7.2, ma qui il vincolo di link è anche un
vincolo di capacità (una sola famiglia di vincoli fa entrambe le cose).
"""
import gurobipy as gp
import pandas as pd
from gurobipy import GRB
from mip import (ammissibile, dualita_forte, due_rilassamenti, frazione,
nuovo_modello, registra_bound, rilassamenti, risolvi,
stampa_soluzione, valuta)
from stile import CICLO, intestazione, plt, salva_dati, salva_figura
from esteso import salva_modello
R = range
# ---------- 1. MODELLO E ISTANZA ----------
intestazione("1. Localizzazione capacitata: dove aprire, quanto spedire")
t1 = [[4, 5, 6], [6, 4, 3]] # costo di trasporto sede l -> cliente c
u1 = [50, 50] # capacita' delle sedi
i1 = [60, 90] # costo di apertura
d1 = [8, 25, 27] # domanda dei clienti
m, n = 2, 3
salva_dati(pd.DataFrame([{"sede": l + 1, "cliente": c + 1, "t": t1[l][c]}
for l in R(m) for c in R(n)]), "fam08_1_costi")
salva_dati(pd.DataFrame({"sede": R(1, m + 1), "u": u1, "i": i1}), "fam08_1_sedi")
salva_dati(pd.DataFrame({"cliente": R(1, n + 1), "d": d1}), "fam08_1_clienti")
def modello_1(t, u, i, d):
m, n = len(u), len(d)
mod = nuovo_modello("localizzazione_capacitata")
x = mod.addVars(m, vtype=GRB.BINARY, name="x")
y = mod.addVars(m, n, name="y")
mod.setObjective(gp.quicksum(i[l] * x[l] for l in R(m))
+ gp.quicksum(t[l][c] * y[l, c] for l in R(m) for c in R(n)), GRB.MINIMIZE)
mod.addConstrs((u[l] * x[l] - gp.quicksum(y[l, c] for c in R(n)) >= 0 for l in R(m)),
name="capacita")
mod.addConstrs((gp.quicksum(y[l, c] for l in R(m)) == d[c] for c in R(n)), name="domanda")
return mod, x, y
def duale_1(t, u, i, d):
"""min sum d_c pi_c; u_l mu_l <= i_l; -mu_l + pi_c <= t_lc; mu >= 0, pi libere."""
m, n = len(u), len(d)
dl = nuovo_modello("duale_localizzazione")
mu = dl.addVars(m, name="mu")
pi = dl.addVars(n, lb=-GRB.INFINITY, name="pi")
dl.setObjective(gp.quicksum(d[c] * pi[c] for c in R(n)), GRB.MAXIMIZE)
dl.addConstrs((u[l] * mu[l] <= i[l] for l in R(m)), name="rc_x")
dl.addConstrs((-mu[l] + pi[c] <= t[l][c] for l in R(m) for c in R(n)), name="rc_y")
return dl
m1, x1, y1 = modello_1(t1, u1, i1, d1)
salva_modello(m1, "fam08_1_primale")
# ---------- 2. IL RILASSAMENTO LP ----------
zlp1, zlp1r, _ = rilassamenti(m1)
# ---------- 3. IL DUALE DEL RILASSAMENTO (LOWER BOUND) ----------
d1_ = duale_1(t1, u1, i1, d1)
salva_modello(d1_, "fam08_1_duale")
mano = {f"mu[{l}]": i1[l] / u1[l] for l in R(m)}
mano.update({f"pi[{c}]": min(t1[l][c] + mano[f"mu[{l}]"] for l in R(m)) for c in R(n)})
lb1, viol = valuta(d1_, mano)
assert viol <= 1e-9, viol
print("Soluzione duale a mano: mu_l = i_l/u_l = " + ", ".join(frazione(i1[l] / u1[l]) for l in R(m))
+ "; pi_c = min_l (t_lc + mu_l) = " + ", ".join(frazione(mano[f"pi[{c}]"]) for c in R(n))
+ f" -> lb = {frazione(lb1)}")
dualita_forte(d1_, zlp1)
# ---------- 4. EURISTICA COSTRUTTIVA (UPPER BOUND) ----------
print("Euristica: si scandiscono le sedi in ordine, riempendo la domanda residua dei clienti")
print("con la capacita' residua di ciascuna sede, senza superare né l'una né l'altra.")
def euristica_1(t, u, i, d):
m, n = len(u), len(d)
y, x, rc, rd, passi = {}, [0] * m, list(u), list(d), []
for l in R(m):
for c in R(n):
if rd[c] > 0 and rc[l] > 0:
q = min(rd[c], rc[l])
y[(l, c)] = q
rd[c] -= q
rc[l] -= q
passi.append(f"Sede {l + 1}, cliente {c + 1}: si spedisce min(rd={rd[c] + q}, rc={rc[l] + q}) = {q}; "
f"rd[{c + 1}] = {rd[c]}, rc[{l + 1}] = {rc[l]}.")
if rc[l] < u[l]:
x[l] = 1
passi.append(f"La sede {l + 1} ha spedito qualcosa (rc = {rc[l]} < u = {u[l]}): si apre, x[{l + 1}] = 1.")
ok = all(v == 0 for v in rd)
return x, y, passi, ok
xe, ye, passi, ok = euristica_1(t1, u1, i1, d1)
for i, s in enumerate(passi, 1):
print(f" Passo {i}. {s}")
assert ok, "euristica non ammissibile: domanda non soddisfatta"
ub1 = sum(i1[l] * xe[l] for l in R(m)) + sum(t1[l][c] * ye.get((l, c), 0) for l in R(m) for c in R(n))
sol_eur = {f"x[{l}]": xe[l] for l in R(m)}
sol_eur.update({f"y[{l},{c}]": v for (l, c), v in ye.items()})
assert ammissibile(m1, sol_eur)
print(f" ub = {ub1}")
# ---------- 5. SOLUZIONE OTTIMA DEL MILP ----------
z1 = risolvi(m1)
print("Soluzione ottima del MILP:")
stampa_soluzione(m1, solo_non_nulle=True)
riga = registra_bound("1 localizzazione capacitata", ub1, lb1, zlp1, zlp1r, z1)
salva_dati(pd.DataFrame([riga]), "fam08_1_bound")
# ---------- 6. DOMANDE DI MODELLAZIONE AGGIUNTIVE ----------
varianti = {}
def variante(nome, mod):
z = risolvi(mod)
print(f" {nome:70s} z = {frazione(z)}")
return z
# 1a: ogni sede aperta deve spedire almeno 5 litri (lotto minimo / semicontinua)
mod, x, y = modello_1(t1, u1, i1, d1)
mod.addConstrs((gp.quicksum(y[l, c] for c in R(n)) >= 5 * x[l] for l in R(m)), name="lotto_minimo")
varianti["1a"] = variante("1a. Ogni sede aperta spedisce almeno 5 litri (sum_c y_lc >= 5 x_l)", mod)
# 1b: la sede 2 si apre solo se si apre la sede 1
mod, x, y = modello_1(t1, u1, i1, d1)
mod.addConstr(x[1] <= x[0], name="2_solo_se_1")
varianti["1b"] = variante("1b. La sede 2 si apre solo se si apre la sede 1 (x_2 <= x_1)", mod)
salva_dati(pd.DataFrame({"variante": list(varianti), "z": list(varianti.values())}), "fam08_1_varianti")
# ---------- 7. IL SANDWICH SULLA VARIANTE 1b ----------
intestazione("1b. Il sandwich sulla variante: la sede 2 si apre solo se si apre la sede 1")
def modello_1b(t, u, i, d):
mod_, xx, yy = modello_1(t, u, i, d)
mod_.addConstr(xx[1] - xx[0] <= 0, name="2_solo_se_1")
return mod_, xx, yy
def duale_1b(t, u, i, d):
"""Al duale di 8.1 si aggiunge rho <= 0 per il vincolo x_2 - x_1 <= 0:
allenta la colonna della sede 1 e stringe quella della sede 2. Il termine
noto e' zero, quindi l'obiettivo non cambia."""
mm, nn = len(u), len(d)
dl = nuovo_modello("duale_localizzazione_1b")
mu = dl.addVars(mm, name="mu")
pi = dl.addVars(nn, lb=-GRB.INFINITY, name="pi")
rho = dl.addVar(lb=-GRB.INFINITY, ub=0.0, name="rho")
dl.setObjective(gp.quicksum(d[c] * pi[c] for c in R(nn)), GRB.MAXIMIZE)
dl.addConstr(u[0] * mu[0] - rho <= i[0], name="rc_x0")
dl.addConstr(u[1] * mu[1] + rho <= i[1], name="rc_x1")
dl.addConstrs((-mu[l] + pi[c] <= t[l][c] for l in R(mm) for c in R(nn)), name="rc_y")
return dl
m1b, x1b, y1b = modello_1b(t1, u1, i1, d1)
salva_modello(m1b, "fam08_1b_primale")
# -- euristica ammissibile: la base, riparata aprendo anche la sede 1 --
print("Euristica costruttiva: si parte dalla soluzione del problema base; se apre la sede 2")
print("senza la 1, si apre anche la 1 (la domanda resta servita, cambia solo il costo fisso).")
aperte = [l for l in R(m) if xe[l]]
print(f" soluzione base: sedi aperte {[l + 1 for l in aperte]}, costo {frazione(ub1)}")
aperte_b = sorted(set(aperte) | ({0} if 1 in aperte else set()))
ub1b = sum(i1[l] for l in aperte_b) + sum(t1[l][c] * ye.get((l, c), 0)
for l in R(m) for c in R(n))
sol_1b = {f"x[{l}]": (1 if l in aperte_b else 0) for l in R(m)}
sol_1b.update({f"y[{l},{c}]": v for (l, c), v in ye.items()})
assert ammissibile(m1b, sol_1b), "la soluzione euristica della variante deve essere ammissibile"
print(f" dopo la riparazione: sedi {[l + 1 for l in aperte_b]} -> ub = {frazione(ub1b)}")
# -- certificato duale: rho sposta costo fisso dalla sede 2 alla sede 1 --
d1b_ = duale_1b(t1, u1, i1, d1)
salva_modello(d1b_, "fam08_1b_duale")
def valore_duale_1b(rho_val):
mu_v = {0: (i1[0] + rho_val) / u1[0], 1: (i1[1] - rho_val) / u1[1]}
if min(mu_v.values()) < 0:
return None, None, None
pi_v = {c: min(t1[l][c] + mu_v[l] for l in R(m)) for c in R(n)}
return sum(d1[c] * pi_v[c] for c in R(n)), mu_v, pi_v
# rho e' non positivo: si prova sulla griglia dei valori che rendono tesa una
# colonna, cioe' dove il minimo che definisce pi_c cambia sede
candidati_1b = [0.0] + [-(u1[0] * (t1[1][c] + i1[1] / u1[1] - t1[0][c]) - i1[0])
* u1[1] / (u1[0] + u1[1]) for c in R(n)]
candidati_1b = [r for r in candidati_1b if r <= 0 and valore_duale_1b(r)[0] is not None]
scelto_1b = max(candidati_1b, key=lambda r: valore_duale_1b(r)[0])
lb1b, mu_1b, pi_1b = valore_duale_1b(scelto_1b)
mano_1b = ({f"mu[{l}]": mu_1b[l] for l in R(m)}
| {f"pi[{c}]": pi_1b[c] for c in R(n)} | {"rho": scelto_1b})
lb1b_val, viol_1b = valuta(d1b_, mano_1b)
assert viol_1b <= 1e-9, viol_1b
print("Soluzione duale a mano: rho sposta costo fisso dalla sede 2 alla sede 1, cioe'")
print(" mu_1 = (i_1 + rho)/u_1 e mu_2 = (i_2 - rho)/u_2; poi pi_c = min_l (t_lc + mu_l)")
print(f" come nel problema base. Si prova rho fra i valori che fanno cambiare quel minimo:")
print(f" rho = {frazione(scelto_1b)} -> lb = {frazione(lb1b_val)}")
if abs(scelto_1b) < 1e-9:
print(" Qui il migliore e' rho = 0: la sede 1 ha gia' il costo per litro piu' basso,")
print(" e spostarle altro costo fisso abbasserebbe mu_2 piu' di quanto alzi mu_1.")
print(" Il certificato della variante coincide con quello del problema base.")
zlp1b, zlp1br, _ = due_rilassamenti(m1b, d1b_)
z1b = risolvi(m1b)
riga_1b = registra_bound("1b sede 2 solo con sede 1", ub1b, lb1b_val, zlp1b, zlp1br, z1b)
salva_dati(pd.DataFrame([riga_1b]), "fam08_1b_bound")
assert lb1b_val <= zlp1b <= z1b <= ub1b + 1e-9
# ---------- 8. FIGURE ----------
def barre_flusso(y, m, n, titolo, nome):
"""Per ogni sede, barra impilata dei litri spediti a ciascun cliente."""
fig, ax = plt.subplots(figsize=(7.2, 3.0))
for l in R(m):
inizio = 0
for c in R(n):
q = y.get((l, c), 0)
if q > 0:
ax.barh(l, q, left=inizio, color=CICLO[c % len(CICLO)], edgecolor="white")
ax.text(inizio + q / 2, l, f"c{c + 1}", ha="center", va="center", color="white",
fontsize=9, fontweight="bold")
inizio += q
ax.set_yticks(R(m))
ax.set_yticklabels([f"sede {l + 1}" for l in R(m)])
ax.set_xlabel("litri spediti")
ax.set_title(titolo)
ax.invert_yaxis()
salva_figura(fig, nome)
ott_y = {(l, c): y1[l, c].X for l in R(m) for c in R(n) if y1[l, c].X > 1e-6}
barre_flusso(ott_y, m, n, f"Localizzazione capacitata: soluzione ottima (z = {frazione(z1)})", "cap08_capacitata_ottimo")
print("Fine.")