Localizzazione di hub con costo massimo

Classe: MILP · Legami: attivazione aggregata, variabile di massimo · Script: python/fam08_4_hub.py

Difficoltà: ★★★★☆ · Tempo: 45–60 min

Problema 8.4

\(n \in \mathbb{Z}_{\ge 1}\) terminali, ciascuno da connettere a esattamente un hub; \(m \in \mathbb{Z}_{\ge 1}\) hub, ciascuno con capacità \(k \in \mathbb{Z}_{\ge 1}\) terminali e costo di attivazione \(f_j \in \mathbb{Q}_{\ge 0}\). \(c_{ij} \in \mathbb{Q}_{\ge 0}\) è il costo di connettere il terminale \(i\) all'hub \(j\). Si minimizza la somma dei costi di attivazione e del costo di connessione massimo di ciascun hub.

Il problema a parole. Decidiamo quali hub attivare e a quale hub connettere ciascun terminale. L'obiettivo: attivazione più, per ciascun hub, il costo di connessione più alto (non la somma). I vincoli: ogni terminale a esattamente un hub; un hub non attivato non serve nessuno, uno attivato serve al più \(k\).

Modello

Variabili decisionali. \(n\,m\) binarie \(x_{ij}\), \(m\) binarie \(y_j\) (hub attivato), \(m\) continue non negative \(z_j\) (costo massimo dell'hub \(j\)).

\[ \begin{aligned} \min ~~ \sum_{j=1}^{m} f_j\, y_j + \sum_{j=1}^{m} z_j & & \\ \text{soggetto a} \quad \sum_{j=1}^{m} x_{ij} &= 1, & \forall i \in \{1, 2, \dots, n\}, \\ -\sum_{i=1}^{n} x_{ij} + k\, y_j &\ge 0, & \forall j \in \{1, 2, \dots, m\}, \\ -c_{ij}\, x_{ij} + z_j &\ge 0, & \forall i \in \{1, 2, \dots, n\},\ \forall j \in \{1, 2, \dots, m\}, \\ x_{ij} &\in \{0, 1\}, & \forall i \in \{1, 2, \dots, n\},\ \forall j \in \{1, 2, \dots, m\}, \\ y_j &\in \{0, 1\}, & \forall j \in \{1, 2, \dots, m\}, \\ z_j &\ge 0, & \forall j \in \{1, 2, \dots, m\}. \end{aligned} \]

l'obiettivo minimizza costi di attivazione più costo massimo per hub;
il primo vincolo assegna ogni terminale a un hub (\(n\) vincoli);
il secondo lega assegnamento e attivazione, in forma aggregata, e impone la capacità (\(m\) vincoli);
il terzo lega assegnamento e variabile di massimo (\(n\,m\) vincoli).

Primo legame: attivazione aggregata. Se un terminale è connesso all'hub \(j\), \(j\) deve essere attivato; dalla contronominale, un hub non attivato non serve nessuno. Entrambi imposti direttamente dal secondo vincolo. Il verso opposto — un hub attivato serve almeno un terminale — non è imposto dai vincoli: \(y_j = 1\) con tutte le \(x_{ij} = 0\) è ammissibile. Segue dall'ottimalità, con una forza che dipende dal segno di \(f_j\): se \(f_j > 0\), spegnere un hub vuoto riduce strettamente il costo, quindi in ogni ottimo nessun hub attivato resta vuoto; se \(f_j = 0\) — che l'enunciato ammette, avendo dichiarato \(f_j \in \mathbb{Q}_{\ge 0}\) — lo scambio non migliora e la conclusione corretta è la più debole «esiste un ottimo in cui gli hub vuoti sono spenti». Sull'istanza \(f = (5,6,7)\) vale la versione forte. Come nel problema 7.2.

Secondo legame: variabile di massimo. Se il terminale \(i\) è connesso a \(j\), \(z_j \ge c_{ij}\): imposto direttamente. All'ottimo, \(z_j = \max_{i:x_{ij}=1} c_{ij}\) esattamente, perché l'obiettivo minimizza \(z_j\) e nessun altro vincolo la coinvolge. Come nel problema 7.4.

Il modello in gurobipy

mod = gp.Model("hub_max")
x = mod.addVars(n, m, vtype=GRB.BINARY, name="x")
y = mod.addVars(m, vtype=GRB.BINARY, name="y")
z = mod.addVars(m, name="z")
mod.setObjective(gp.quicksum(f[j] * y[j] for j in range(m)) + z.sum(), GRB.MINIMIZE)
mod.addConstrs((gp.quicksum(x[i, j] for j in range(m)) == 1 for i in range(n)), name="assegnamento")
mod.addConstrs((-gp.quicksum(x[i, j] for i in range(n)) + k * y[j] >= 0 for j in range(m)), name="attivazione")
mod.addConstrs((-c[i][j] * x[i, j] + z[j] >= 0 for i in range(n) for j in range(m)), name="massimo")

L'istanza

\(n=3\) terminali, \(m=3\) hub, \(k=2\):

\(c_{ij}\)	\(j=1\)	\(j=2\)	\(j=3\)
\(i=1\)	5	10	2
\(i=2\)	5	4	6
\(i=3\)	5	4	6

	\(j=1\)	\(j=2\)	\(j=3\)
\(f_j\)	5	6	7

Il modello scritto sui dati dell'istanza:

\[ \begin{array}{rrrrrrrrrrrrrrrr c l} \min & & & & & & & & & & 5y_1 & +6y_2 & +7y_3 & +z_1 & +z_2 & +z_3 & & \\ \text{soggetto a} & x_{11} & +x_{12} & +x_{13} & & & & & & & & & & & & & = & 1\\ & & & & x_{21} & +x_{22} & +x_{23} & & & & & & & & & & = & 1\\ & & & & & & & x_{31} & +x_{32} & +x_{33} & & & & & & & = & 1\\ & -x_{11} & & & -x_{21} & & & -x_{31} & & & +2y_1 & & & & & & \ge & 0\\ & & -x_{12} & & & -x_{22} & & & -x_{32} & & & +2y_2 & & & & & \ge & 0\\ & & & -x_{13} & & & -x_{23} & & & -x_{33} & & & +2y_3 & & & & \ge & 0\\ & -5x_{11} & & & & & & & & & & & & +z_1 & & & \ge & 0\\ & & -10x_{12} & & & & & & & & & & & & +z_2 & & \ge & 0\\ & & & -2x_{13} & & & & & & & & & & & & +z_3 & \ge & 0\\ & & & & -5x_{21} & & & & & & & & & +z_1 & & & \ge & 0\\ & & & & & -4x_{22} & & & & & & & & & +z_2 & & \ge & 0\\ & & & & & & -6x_{23} & & & & & & & & & +z_3 & \ge & 0\\ & & & & & & & -5x_{31} & & & & & & +z_1 & & & \ge & 0\\ & & & & & & & & -4x_{32} & & & & & & +z_2 & & \ge & 0\\ & & & & & & & & & -6x_{33} & & & & & & +z_3 & \ge & 0\\ & x_{11}, & x_{12}, & x_{13}, & x_{21}, & x_{22}, & x_{23}, & x_{31}, & x_{32}, & x_{33} & & & & & & & \in & \{0, 1\}\\ & & & & & & & & & & y_1, & y_2, & y_3 & & & & \in & \{0, 1\}\\ & & & & & & & & & & & & & z_1, & z_2, & z_3 & \ge & 0 \end{array} \]

Euristica costruttiva: il bound primale

Un next-fit (bin packing): un hub alla volta, fino a \(k\) terminali — la stessa euristica generica dello scheduling, riusata da euristiche.py. Terminale 1 e 2 sull'hub 1 (pieno), terminale 3 sull'hub 2. Costi massimi: \(z_1=\max(5,5)=5\), \(z_2=4\). Valore \(5+6+5+4=20\): \(z(\mathit{MILP}) \le \mathit{UB} = 20\).

Rilassamento LP e duale: il bound duale

Il duale del rilassamento lineare, una variabile per vincolo del primale:

\[ \begin{aligned} \max ~~ \sum_{i=1}^{n} \alpha_i & & \\ \text{soggetto a} \quad \alpha_i - \beta_j - c_{ij}\, \gamma_{ij} &\le 0, & \forall i \in \{1, 2, \dots, n\},\ \forall j \in \{1, 2, \dots, m\}, \\ k\, \beta_j &\le f_j, & \forall j \in \{1, 2, \dots, m\}, \\ \sum_{i=1}^{n} \gamma_{ij} &\le 1, & \forall j \in \{1, 2, \dots, m\}, \\ \alpha_i &\gtreqless 0, & \forall i \in \{1, 2, \dots, n\}, \\ \beta_j &\ge 0, & \forall j \in \{1, 2, \dots, m\}, \\ \gamma_{ij} &\ge 0, & \forall i \in \{1, 2, \dots, n\},\ \forall j \in \{1, 2, \dots, m\}. \end{aligned} \]

Lo stesso duale, scritto sui dati dell'istanza:

\[ \begin{array}{rrrrrrrrrrrrrrrr c l} \max & \alpha_1 & +\alpha_2 & +\alpha_3 & & & & & & & & & & & & & & \\ \text{soggetto a} & \alpha_1 & & & -\beta_1 & & & -5\gamma_{11} & & & & & & & & & \le & 0\\ & \alpha_1 & & & & -\beta_2 & & & -10\gamma_{12} & & & & & & & & \le & 0\\ & \alpha_1 & & & & & -\beta_3 & & & -2\gamma_{13} & & & & & & & \le & 0\\ & & \alpha_2 & & -\beta_1 & & & & & & -5\gamma_{21} & & & & & & \le & 0\\ & & \alpha_2 & & & -\beta_2 & & & & & & -4\gamma_{22} & & & & & \le & 0\\ & & \alpha_2 & & & & -\beta_3 & & & & & & -6\gamma_{23} & & & & \le & 0\\ & & & \alpha_3 & -\beta_1 & & & & & & & & & -5\gamma_{31} & & & \le & 0\\ & & & \alpha_3 & & -\beta_2 & & & & & & & & & -4\gamma_{32} & & \le & 0\\ & & & \alpha_3 & & & -\beta_3 & & & & & & & & & -6\gamma_{33} & \le & 0\\ & & & & 2\beta_1 & & & & & & & & & & & & \le & 5\\ & & & & & 2\beta_2 & & & & & & & & & & & \le & 6\\ & & & & & & 2\beta_3 & & & & & & & & & & \le & 7\\ & & & & & & & \gamma_{11} & & & +\gamma_{21} & & & +\gamma_{31} & & & \le & 1\\ & & & & & & & & \gamma_{12} & & & +\gamma_{22} & & & +\gamma_{32} & & \le & 1\\ & & & & & & & & & \gamma_{13} & & & +\gamma_{23} & & & +\gamma_{33} & \le & 1\\ & \alpha_1, & \alpha_2, & \alpha_3 & & & & & & & & & & & & & \gtreqless & 0\\ & & & & \beta_1, & \beta_2, & \beta_3 & & & & & & & & & & \ge & 0\\ & & & & & & & \gamma_{11}, & \gamma_{12}, & \gamma_{13}, & \gamma_{21}, & \gamma_{22}, & \gamma_{23}, & \gamma_{31}, & \gamma_{32}, & \gamma_{33} & \ge & 0 \end{array} \]

Con \(\bar\gamma_{ij}=0\) e \(\bar\beta_j = f_j/k\) (il massimo ammesso), il vincolo su \(\alpha_i\) vale per ogni hub \(j\), non solo il più conveniente: \(\bar\alpha_i = \min_j \bar\beta_j\).

\[ \bar\beta = (5/2,\ 3,\ 7/2),\qquad \bar\alpha_i = 5/2\ \ \forall i, \]

di valore \(3\cdot5/2=15/2\). Per la dualità debole, \(\mathit{LB}=15/2 \le z(\mathit{LP}) \le z(\mathit{MILP}) \le \mathit{UB}=20\).

Un tranello frequente

Il vincolo su \(\alpha_i\) vale per ogni hub \(j\): fissare \(\bar\gamma_{ij}=0\) solo per gli hub «non convenienti» non basta a liberare \(\alpha_i\) da quel vincolo. \(\alpha_i\) resta limitato dal minimo su tutti gli hub, non da uno solo.

Quello che dice il solver. \(z(\mathit{LP})=25/2\), \(z(\mathit{LP}^+)=1015/78\approx13{,}0\). \(z(\mathit{MILP})=19\), con gli hub 1 e 3 attivati (non 1 e 2): il terminale 1 da solo sull'hub 3 (il più economico per lui), i terminali 2 e 3 sull'hub 1. Gap euristica \(5{,}3\%\).

\(UB\)	\(LB\) (duale)	\(z(\mathit{LP})\)	\(z(\mathit{LP}^+)\)	\(z(\mathit{MILP})\)	gap dell'euristica
20	\(15/2\)	\(25/2\)	\(1015/78\)	19	\(5{,}3\%\)

Soluzione ottima

Considerazioni aggiuntive

\(x_{ij} \le y_j\) (disaggregato) è implicato dal vincolo aggregato di attivazione sui punti interi, non nel rilassamento: aggiungerlo non cambia \(z(\mathit{MILP})\) e alza \(z(\mathit{LP}^+)\) da \(1015/78\) a \(79/6\) (domanda 8.4.1).
Con \(M_j=\max_i c_{ij}\), \(z_j \le M_j y_j\) non è una disuguaglianza valida (il modello ammette \(z_j>0\) con \(y_j=0\)), ma è un vincolo che preserva l'ottimalità: minimizzando \(z_j\), l'ottimo la annulla comunque quando \(y_j=0\).

Domande di modellazione aggiuntive

8.4.1 — Link di attivazione disaggregato

Si aggiungano al modello i link disaggregati \(x_{ij} \le y_j\). Cambia l'ottimo? Cambia il rilassamento? E che cosa succede se, invece di aggiungerli, si sostituisce con essi il vincolo aggregato?

Una variante svolta: connessione vietata

Il terminale 1 non può essere connesso all'hub 2 (un vincolo di sicurezza).

Si fissa la variabile a zero con il vincolo lineare

\[ x_{12} = 0 \]

(un vincolo lineare). Sull'istanza l'ottimo del problema 8.4 non usa già \(x_{12}\) (il terminale 1 è connesso all'hub 3), quindi il vincolo aggiuntivo non è vincolante e l'ottimo resta \(19\).

Vietare una connessione significa togliere una colonna dal primale, quindi togliere un vincolo dal duale: \(\alpha_1\) non deve più reggere il confronto con l'hub proibito e può salire. Gli si dà un \(\gamma\) su ciascuno degli hub che gli restano — il budget di ogni hub è \(1\) e nessun altro terminale lo usa — e il bound cresce. È il caso in cui un vincolo in più nel primale migliora il certificato, invece di lasciarlo fermo.

	valore	che cos'è
\(\mathit{UB}\)	\(20\)	soluzione euristica
\(\mathit{LB}\)	\(\frac{21}{2}\)	certificato duale costruito a mano
\(z(\mathit{LP})\)	\(\frac{25}{2}\)	rilassamento senza i bound
\(z(\mathit{LP}^+)\)	\(\frac{79}{6}\)	rilassamento con i bound
\(z(\mathit{MILP})\)	\(19\)	ottimo del MILP

Codice

Script completo — python/fam08_4_hub.py (riproducibile con python3 python/fam08_4_hub.py dalla cartella python/, richiama next_fit da euristiche.py). Notebook — notebooks/fam08_4_hub.ipynb — che si apre in Colab dal badge in cima alla pagina.

Mostra lo script completo — python/fam08_4_hub.py (249 righe)

"""Problema 8.4 -- Localizzazione di hub con costo di connessione massimo.

Due link: attivazione (aggregata, come nello scheduling 7.2) e variabile di
massimo z_j = max_i {c_ij : x_ij = 1} (stesso schema del tempo di lavorazione 7.4).
L'euristica
next-fit è quella generica di euristiche.py: gli hub sono le "macchine" (capacità
k) e i terminali i "lavori" (tempo unitario, indipendente dalla macchina).
"""
import gurobipy as gp
import pandas as pd
from gurobipy import GRB

from euristiche import matrice, next_fit
from mip import (ammissibile, dualita_forte, due_rilassamenti, frazione,
                 nuovo_modello, registra_bound, rilassamenti, rilassamento, risolvi,
                 stampa_soluzione, valuta)
from stile import intestazione, plt, salva_dati, salva_figura
from esteso import salva_modello

R = range

# ---------- 1. MODELLO E ISTANZA ----------

intestazione("4. Localizzazione di hub: attivazione e costo di connessione massimo")
c4 = [[5, 10, 2], [5, 4, 6], [5, 4, 6]]   # costo di connessione terminale i -> hub j
f4 = [5, 6, 7]                             # costo di attivazione hub j
k4 = 2                                     # capacità di ciascun hub
n, m = 3, 3
salva_dati(pd.DataFrame([{"terminale": i + 1, "hub": j + 1, "c": c4[i][j]}
                         for i in R(n) for j in R(m)]), "fam08_4_costi")
salva_dati(pd.DataFrame({"hub": R(1, m + 1), "f": f4}), "fam08_4_attivazione")


def modello_4(c, f, k):
    n, m = len(c), len(f)
    mod = nuovo_modello("hub_max")
    x = mod.addVars(n, m, vtype=GRB.BINARY, name="x")
    y = mod.addVars(m, vtype=GRB.BINARY, name="y")
    z = mod.addVars(m, name="z")
    mod.setObjective(gp.quicksum(f[j] * y[j] for j in R(m)) + z.sum(), GRB.MINIMIZE)
    mod.addConstrs((gp.quicksum(x[i, j] for j in R(m)) == 1 for i in R(n)), name="assegnamento")
    mod.addConstrs((-gp.quicksum(x[i, j] for i in R(n)) + k * y[j] >= 0 for j in R(m)), name="attivazione")
    mod.addConstrs((-c[i][j] * x[i, j] + z[j] >= 0 for i in R(n) for j in R(m)), name="massimo")
    return mod, x, y, z


def duale_4(c, f, k):
    """max sum_i alpha_i;  alpha_i - beta_j - c_ij gamma_ij <= 0;  k beta_j <= f_j;
    sum_i gamma_ij <= 1;  alpha libero, beta,gamma >= 0."""
    n, m = len(c), len(f)
    dl = nuovo_modello("duale_hub")
    alpha = dl.addVars(n, lb=-GRB.INFINITY, name="alpha")
    beta = dl.addVars(m, name="beta")
    gamma = dl.addVars(n, m, name="gamma")
    dl.setObjective(alpha.sum(), GRB.MAXIMIZE)
    dl.addConstrs((alpha[i] - beta[j] - c[i][j] * gamma[i, j] <= 0 for i in R(n) for j in R(m)), name="rc_x")
    dl.addConstrs((k * beta[j] <= f[j] for j in R(m)), name="rc_y")
    dl.addConstrs((gp.quicksum(gamma[i, j] for i in R(n)) <= 1 for j in R(m)), name="rc_z")
    return dl


m4, x4, y4, z4 = modello_4(c4, f4, k4)
salva_modello(m4, "fam08_4_primale")

# ---------- 2. IL RILASSAMENTO LP ----------
zlp4, zlp4r, _ = rilassamenti(m4)

# ---------- 3. IL DUALE DEL RILASSAMENTO (LOWER BOUND) ----------

d4 = duale_4(c4, f4, k4)
salva_modello(d4, "fam08_4_duale")
beta_mano = [f4[j] / k4 for j in R(m)]     # il massimo ammesso da k*beta_j <= f_j
alpha_mano = min(beta_mano)                # deve reggere per OGNI hub j, non solo il più conveniente
mano = {f"gamma[{i},{j}]": 0.0 for i in R(n) for j in R(m)}
mano.update({f"beta[{j}]": beta_mano[j] for j in R(m)})
mano.update({f"alpha[{i}]": alpha_mano for i in R(n)})
lb4, viol = valuta(d4, mano)
assert viol <= 1e-9, viol
print(f"Soluzione duale a mano: gamma = 0, beta_j = f_j/k = {[frazione(b) for b in beta_mano]}, "
      f"alpha_i = min_j beta_j = {frazione(alpha_mano)}  ->  lb = {frazione(lb4)}")
dualita_forte(d4, zlp4)

# ---------- 4. EURISTICA COSTRUTTIVA (UPPER BOUND) ----------

print("Euristica next-fit: si riempiono gli hub uno alla volta fino a k terminali,")
print("poi si passa al successivo (stessa euristica generica dei problemi di scheduling).")
t4 = matrice([1] * n, m)   # tempo unitario per ogni terminale, indipendente dall'hub
a4 = [k4] * m                    # capacità residua di ciascun hub
esito4 = next_fit(t4, a4)
esito4.traccia.stampa()
assert esito4.ok
ye = esito4.y
ze = [0.0] * m
for j in R(m):
    if ye[j]:
        ze[j] = max(c4[i][j] for i in R(n) if esito4.x.get((i, j)) == 1)
ub4 = sum(f4[j] * ye[j] for j in R(m)) + sum(ze)
print(f"  y = {ye}, z = {ze}  ->  ub = {frazione(ub4)}")

# ---------- 5. SOLUZIONE OTTIMA DEL MILP ----------

z4v = risolvi(m4)
print("Soluzione ottima del MILP:")
stampa_soluzione(m4, solo_non_nulle=True)
riga = registra_bound("4 hub", ub4, lb4, zlp4, zlp4r, z4v, senso="min")
salva_dati(pd.DataFrame([riga]), "fam08_4_bound")

# ---------- 6. DOMANDE DI MODELLAZIONE AGGIUNTIVE ----------

varianti = {}


def variante(nome, mod):
    z = risolvi(mod)
    print(f"  {nome:70s} z = {frazione(z)}")
    return z


# 4a: si AGGIUNGONO i link disaggregati x_ij <= y_j al vincolo aggregato
mod, x, y, z = modello_4(c4, f4, k4)
mod.addConstrs((x[i, j] <= y[j] for i in R(n) for j in R(m)), name="attivazione_disaggregata")
varianti["4a"] = variante("4a. Link disaggregati AGGIUNTI a quello aggregato (x_ij <= y_j)", mod)
zlp_4a, _, _ = rilassamento(mod, rafforzato=True)
zlp_base, _, _ = rilassamento(modello_4(c4, f4, k4)[0], rafforzato=True)
print(f"      rilassamento: z(LP+) passa da {frazione(zlp_base)} a {frazione(zlp_4a)}: i link")
print("      disaggregati sono disuguaglianze valide implicate dal vincolo aggregato sui")
print("      punti interi, ma non dal rilassamento, e lo rafforzano.")

# 4a-bis: il tranello. Se si SOSTITUISCE il vincolo aggregato con i soli link
# disaggregati, si perde anche la capacita' k: il modello non e' piu' quello del
# problema. Va tenuta esplicitamente, oppure si parla di aggiunta e non di sostituzione.
mod, x, y, z = modello_4(c4, f4, k4)
mod.update()
mod.remove([cc for cc in mod.getConstrs() if cc.ConstrName.startswith("attivazione")])
mod.update()
mod.addConstrs((x[i, j] <= y[j] for i in R(n) for j in R(m)), name="solo_disaggregato")
varianti["4a_senza_capacita"] = variante(
    "4a'. SOSTITUENDO l'aggregato con i soli link disaggregati (capacita' persa)", mod)
mod, x, y, z = modello_4(c4, f4, k4)
mod.update()
mod.remove([cc for cc in mod.getConstrs() if cc.ConstrName.startswith("attivazione")])
mod.update()
mod.addConstrs((x[i, j] <= y[j] for i in R(n) for j in R(m)), name="solo_disaggregato")
mod.addConstrs((gp.quicksum(x[i, j] for i in R(n)) <= k4 for j in R(m)), name="capacita")
varianti["4a_con_capacita"] = variante(
    "4a''. Sostituzione corretta: link disaggregati + capacita' separata", mod)
assert varianti["4a_senza_capacita"] < varianti["4a"], "senza capacita' l'ottimo scende"
assert varianti["4a_con_capacita"] == varianti["4a"], "con la capacita' l'ottimo non cambia"
# 4b: il terminale 1 non può essere connesso all'hub 2
mod, x, y, z = modello_4(c4, f4, k4)
mod.addConstr(x[0, 1] == 0, name="terminale1_non_hub2")
varianti["4b"] = variante("4b. Il terminale 1 non può connettersi all'hub 2 (x_12 = 0)", mod)
salva_dati(pd.DataFrame({"variante": list(varianti), "z": list(varianti.values())}), "fam08_4_varianti")

# ---------- 7. IL SANDWICH SULLA VARIANTE 4b ----------
intestazione("4b. Il sandwich sulla variante: il terminale 1 non puo' usare l'hub 2")

VIETATA = (0, 1)      # (terminale, hub) proibito


def modello_4b(c, f, k, vietata=VIETATA):
    mod_, xx, yy, zz = modello_4(c, f, k)
    mod_.addConstr(xx[vietata] == 0, name="connessione_vietata")
    return mod_, xx, yy, zz


def duale_4b(c, f, k, vietata=VIETATA):
    """Vietare una connessione vuol dire togliere una colonna dal primale,
    quindi togliere il vincolo corrispondente dal duale: alpha_1 non deve piu'
    reggere il confronto con l'hub 2, e puo' salire."""
    nn, mm = len(c), len(f)
    dl = nuovo_modello("duale_hub_4b")
    alpha = dl.addVars(nn, lb=-GRB.INFINITY, name="alpha")
    beta = dl.addVars(mm, name="beta")
    gamma = dl.addVars(nn, mm, name="gamma")
    dl.setObjective(alpha.sum(), GRB.MAXIMIZE)
    dl.addConstrs((alpha[i] - beta[j] - c[i][j] * gamma[i, j] <= 0
                   for i in R(nn) for j in R(mm) if (i, j) != vietata), name="rc_x")
    dl.addConstrs((k * beta[j] <= f[j] for j in R(mm)), name="rc_y")
    dl.addConstrs((gp.quicksum(gamma[i, j] for i in R(nn)) <= 1 for j in R(mm)), name="rc_z")
    return dl


m4b, x4b, y4b, z4b = modello_4b(c4, f4, k4)
salva_modello(m4b, "fam08_4b_primale")

# -- euristica ammissibile: la base, spostando il terminale vietato --
print("Euristica costruttiva: si parte dalla soluzione del problema base e, se il terminale")
print("vietato sta sull'hub proibito, lo si sposta sull'hub capiente di costo piu' basso.")
ass_4b = {i: j for (i, j), v in esito4.x.items() if v == 1}
ti, hj = VIETATA
if ass_4b.get(ti) == hj:
    carichi = {j: sum(1 for q, jj in ass_4b.items() if jj == j and q != ti) for j in R(m)}
    nuovo_hub = min((j for j in R(m) if j != hj and carichi[j] < k4), key=lambda j: c4[ti][j])
    print(f"  il terminale {ti + 1} era sull'hub {hj + 1}: si sposta sull'hub {nuovo_hub + 1}")
    ass_4b[ti] = nuovo_hub
else:
    print(f"  il terminale {ti + 1} non era sull'hub {hj + 1}: niente da riparare")
y_4b = [1 if any(j == q for q in ass_4b.values()) else 0 for j in R(m)]
z_4b = [max([c4[i][j] for i, q in ass_4b.items() if q == j] + [0.0]) for j in R(m)]
ub4b = sum(f4[j] * y_4b[j] for j in R(m)) + sum(z_4b)
sol_4b = ({f"x[{i},{j}]": (1 if ass_4b[i] == j else 0) for i in R(n) for j in R(m)}
          | {f"y[{j}]": y_4b[j] for j in R(m)} | {f"z[{j}]": z_4b[j] for j in R(m)})
assert ammissibile(m4b, sol_4b), "la soluzione euristica della variante deve essere ammissibile"
print(f"  ub = {frazione(ub4b)}")

# -- certificato duale: senza quella colonna, alpha_1 sale --
d4b = duale_4b(c4, f4, k4)
salva_modello(d4b, "fam08_4b_duale")
beta_4b = {j: f4[j] / k4 for j in R(m)}
# il terminale vietato puo' spendere un gamma su ogni hub che gli resta: ogni hub
# ha budget 1 e nessun altro terminale lo usa nella ricetta
mano_4b = {f"beta[{j}]": beta_4b[j] for j in R(m)}
mano_4b.update({f"gamma[{i},{j}]": 0.0 for i in R(n) for j in R(m)})
mano_4b.update({f"gamma[{ti},{j}]": 1.0 for j in R(m) if j != hj})
alpha_vietato = min(beta_4b[j] + c4[ti][j] for j in R(m) if j != hj)
alpha_altri = min(beta_4b.values())
mano_4b.update({f"alpha[{i}]": (alpha_vietato if i == ti else alpha_altri) for i in R(n)})
lb4b, viol_4b = valuta(d4b, mano_4b)
assert viol_4b <= 1e-9, viol_4b
print("Soluzione duale a mano: beta_j = f_j/k come nel problema base. Per i terminali")
print("  liberi alpha_i = min_j beta_j. Il terminale vietato, invece, non deve piu' reggere")
print("  il confronto con l'hub proibito: gli si da' un gamma su ciascuno degli hub che gli")
print("  restano (il budget di ogni hub e' 1 e nessun altro lo usa), e allora")
print(f"  alpha_{ti + 1} = min_(j != {hj + 1}) (beta_j + c_{ti + 1}j) = {frazione(alpha_vietato)}")
print(f"  invece di {frazione(alpha_altri)}.")
print(f"  ->  lb = {frazione(lb4b)}  (la ricetta del problema base darebbe "
      f"{frazione(n * alpha_altri)})")
zlp4b, zlp4br, _ = due_rilassamenti(m4b, d4b)
z4b_val = risolvi(m4b)
riga_4b = registra_bound("4b connessione vietata", ub4b, lb4b, zlp4b, zlp4br, z4b_val)
salva_dati(pd.DataFrame([riga_4b]), "fam08_4b_bound")
assert lb4b <= zlp4b <= z4b_val <= ub4b + 1e-9

# ---------- 8. FIGURE ----------

fig, ax = plt.subplots(figsize=(6.4, 3.2))
colori = ["#16324A", "#0E7490", "#CA6F1E"]
for j in R(m):
    if y4[j].X > 0.5:
        assegnati = [i + 1 for i in R(n) if x4[i, j].X > 0.5]
        ax.barh(j, z4[j].X, color=colori[j % 3], label=f"hub {j + 1}: terminali {assegnati}")
ax.set_yticks(R(m))
ax.set_yticklabels([f"hub {j + 1}" for j in R(m)])
ax.set_xlabel("costo di connessione massimo $z_j$")
ax.set_title(f"Soluzione ottima (z = {frazione(z4v)})")
ax.legend(fontsize=7, loc="lower right")
salva_figura(fig, "cap08_hub_ottimo")
print("Fine.")