Copertura del segnale con interferenza
Classe: BIP · Legami: se e solo se (soglia + interferenza) · Script: python/fam08_3_copertura.py
Difficoltà: ★★★★☆ · Tempo: 30–45 min
Problema 8.3
Un operatore sceglie al più \(k \in \mathbb{Z}_{\ge 1}\) sedi, fra \(m \in \mathbb{Z}_{\ge 1}\) candidate, per servire \(n \in \mathbb{Z}_{\ge 1}\) clienti. \(s_{lc} \in \mathbb{Q}_{\ge 0}\) è il segnale ricevuto dal cliente \(c\) se \(l\) è installata. Un cliente è coperto se e solo se il segnale totale è almeno \(t \in \mathbb{Q}_{>0}\) e al più una sede genera per lui un segnale \(\ge b \in \mathbb{Q}_{>0}\). \(p_c \in \mathbb{Q}_{>0}\) è il profitto se coperto. Si vuole massimizzare il profitto totale.
Il problema a parole. Decidiamo quali sedi installare (al più \(k\)). L'obiettivo: profitto totale massimo. I vincoli: un cliente è coperto se e solo se riceve segnale sufficiente e non troppa interferenza; al più \(k\) sedi installate.
Modello
Dati. \(m\), \(n\), \(s_{lc} \in \mathbb{Q}_{\ge 0}\), \(p_c \in \mathbb{Q}_{>0}\), soglia \(t\), limite di interferenza \(b\), budget \(k\). Per ogni cliente \(c\): \(\mathscr{L}_c = \{l : s_{lc} \ge b\}\).
Variabili decisionali. \(m\) binarie \(x_l\) (sede installata), \(n\) binarie \(y_c\) (cliente coperto).
- l'obiettivo massimizza il profitto totale;
- il primo vincolo lega copertura e segnale ricevuto (\(n\) vincoli);
- il secondo lega copertura e interferenza (\(n\) vincoli);
- il terzo limita a \(k\) le sedi installate (un vincolo).
Il legame: un se e solo se. Un verso — \(y_c=1 \Rightarrow\) segnale \(\ge t\) e al più una sede forte — è imposto direttamente dai due vincoli. L'altro verso — se entrambe le condizioni valgono, il cliente è coperto — non è imposto dai vincoli (che ammettono anche \(y_c=0\)), ma segue dall'ottimalità: poiché \(p_c>0\) e \(y_c\) compare solo in questi due vincoli, alzarla a \(1\) resta ammissibile e aumenta l'obiettivo. Lo stesso schema del problema 7.6.
Il modello in gurobipy
mod = gp.Model("copertura_interferenza")
x = mod.addVars(m, vtype=GRB.BINARY, name="x")
y = mod.addVars(n, vtype=GRB.BINARY, name="y")
mod.setObjective(gp.quicksum(p[c] * y[c] for c in range(n)), GRB.MAXIMIZE)
mod.addConstrs((-gp.quicksum(s[l][c] * x[l] for l in range(m)) + t * y[c] <= 0
for c in range(n)), name="soglia")
mod.addConstrs((gp.quicksum(x[l] for l in L[c]) + (m - 1) * y[c] <= m
for c in range(n)), name="interferenza")
mod.addConstr(x.sum() <= k, name="budget")
L'istanza
\(m=3\), \(n=5\), \(t=5\), \(b=4\), \(k=2\):
| \(s_{lc}\) | \(c=1\) | \(c=2\) | \(c=3\) | \(c=4\) | \(c=5\) |
|---|---|---|---|---|---|
| \(l=1\) | 6 | 0 | 5 | 3 | 1 |
| \(l=2\) | 4 | 5 | 2 | 0 | 0 |
| \(l=3\) | 0 | 7 | 5 | 4 | 2 |
| \(c=1\) | \(c=2\) | \(c=3\) | \(c=4\) | \(c=5\) | |
|---|---|---|---|---|---|
| \(p_c\) | 10 | 20 | 5 | 15 | 25 |
Con \(b=4\): \(\mathscr{L}_1=\{1,2\}\), \(\mathscr{L}_2=\{2,3\}\), \(\mathscr{L}_3=\{1,3\}\), \(\mathscr{L}_4=\{3\}\), \(\mathscr{L}_5=\emptyset\).
Il modello scritto sui dati dell'istanza:
Euristica costruttiva: il bound primale
Si aprono le prime \(k\) sedi. Cliente 1: segnale \(10\ge5\) ma 2 sedi forti (\(>1\)): non coperto. Cliente 2: segnale \(5\ge5\), 1 sede forte: coperto. Cliente 3: segnale \(7\ge5\), 1 sede forte: coperto. Clienti 4 e 5: segnale insufficiente: non coperti. Valore \(20+5=25\): \(z(\mathit{MILP}) \ge \mathit{LB} = 25\).
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\pi_c=0\), \(\bar\mu=0\) e \(\bar\lambda_c = p_c/(m-1) = p_c/2\):
di valore \(m\sum_c\bar\lambda_c = 3\cdot75/2=225/2\). Per la dualità debole (problema di massimo: l'euristica dà il lower bound, il duale l'upper bound), \(\mathit{LB}=25 \le z(\mathit{MILP}) \le z(\mathit{LP}) \le \mathit{UB}=225/2\).
Quello che dice il solver. \(z(\mathit{LP}) = 41925/646 \approx 64{,}9\), \(z(\mathit{LP}^+) = 125/2 = 62{,}5\). \(z(\mathit{MILP}) = 45\), con le sedi 1 e 3 installate e i clienti 1, 2, 4 coperti (non 3 né 5): diverso da quanto trovato dall'euristica. Gap euristica \(44{,}4\%\).
| \(LB\) | \(UB\) (duale) | \(z(\mathit{LP})\) | \(z(\mathit{LP}^+)\) | \(z(\mathit{MILP})\) | gap dell'euristica |
|---|---|---|---|---|---|
| 25 | \(225/2\) | \(41925/646\) | \(125/2\) | 45 | \(44{,}4\%\) |
Considerazioni aggiuntive
- Il cliente 5 non può mai essere coperto: segnale massimo \(1+0+2=3<5\) anche aprendo tutte le sedi.
- Per i clienti con \(|\mathscr{L}_c|\le1\) (4 e 5) il vincolo di interferenza è ridondante.
Domande di modellazione aggiuntive
8.3.1 — Installazione condizionata
La sede 1 può essere installata solo se lo è anche la sede 3. Come si modella? Qual è il nuovo ottimo?
Una variante svolta: copertura minima garantita
Per contratto, almeno \(3\) clienti devono essere coperti.
Si aggiunge il vincolo lineare
(un vincolo lineare). Sull'istanza l'ottimo del problema 8.3 copre già \(3\) clienti, quindi il vincolo non è vincolante e l'ottimo resta \(45\).
«Almeno tre clienti coperti» aggiunge al duale \(\omega \le 0\). Conviene lasciarlo a zero, e il conto dice perché: abbassarlo costringe ogni \(\lambda_c\) a salire di \(-\omega/(m-1)\), e nell'obiettivo quelle \(\lambda\) pesano più di quanto il termine noto faccia risparmiare. L'euristica, con poche sedi, prova ogni scelta di \(k\) e tiene la migliore fra quelle ammissibili.
| valore | che cos'è | |
|---|---|---|
| \(\mathit{UB}\) | \(\frac{225}{2}\) | certificato duale costruito a mano |
| \(\mathit{LB}\) | \(45\) | soluzione euristica |
| \(z(\mathit{LP})\) | \(\frac{41925}{646}\) | rilassamento senza i bound |
| \(z(\mathit{LP}^+)\) | \(\frac{125}{2}\) | rilassamento con i bound |
| \(z(\mathit{MILP})\) | \(45\) | ottimo del MILP |
Codice
Script completo —
python/fam08_3_copertura.py
(riproducibile con python3 python/fam08_3_copertura.py dalla cartella
python/). Notebook —
notebooks/fam08_3_copertura.ipynb
— che si apre in Colab dal badge in cima alla pagina.
Mostra lo script completo — python/fam08_3_copertura.py (227 righe)
"""Problema 8.3 -- Copertura del segnale con interferenza (massimo profitto).
Un «se e solo se» come nel problema 7.6: un verso (soglia+interferenza
=> coperto) è imposto da due famiglie di vincoli di link; l'altro verso
(coperto => condizioni soddisfatte) segue dall'obiettivo.
"""
from itertools import combinations as combinazioni
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 intestazione, plt, salva_dati, salva_figura
from esteso import salva_modello
R = range
# ---------- 1. MODELLO E ISTANZA ----------
intestazione("3. Copertura con interferenza: soglia di segnale e al più una sede forte")
s3 = [[6, 0, 5, 3, 1], [4, 5, 2, 0, 0], [0, 7, 5, 4, 2]] # segnale sede l -> cliente c
p3 = [10, 20, 5, 15, 25] # profitto se il cliente c è coperto
t3, b3, k3 = 5, 4, 2 # soglia di segnale, limite di interferenza, budget sedi
m, n = 3, 5
L3 = [[l for l in R(m) if s3[l][c] >= b3] for c in R(n)] # L_c: sedi "forti" per il cliente c
salva_dati(pd.DataFrame([{"sede": l + 1, "cliente": c + 1, "s": s3[l][c]}
for l in R(m) for c in R(n)]), "fam08_3_segnale")
salva_dati(pd.DataFrame({"cliente": R(1, n + 1), "p": p3}), "fam08_3_clienti")
def modello_3(s, p, t, b, k):
m, n = len(s), len(p)
L = [[l for l in R(m) if s[l][c] >= b] for c in R(n)]
mod = nuovo_modello("copertura_interferenza")
x = mod.addVars(m, vtype=GRB.BINARY, name="x")
y = mod.addVars(n, vtype=GRB.BINARY, name="y")
mod.setObjective(gp.quicksum(p[c] * y[c] for c in R(n)), GRB.MAXIMIZE)
mod.addConstrs((-gp.quicksum(s[l][c] * x[l] for l in R(m)) + t * y[c] <= 0 for c in R(n)),
name="soglia")
mod.addConstrs((gp.quicksum(x[l] for l in L[c]) + (m - 1) * y[c] <= m for c in R(n)),
name="interferenza")
mod.addConstr(x.sum() <= k, name="budget")
return mod, x, y, L
def duale_3(s, p, t, b, k):
"""min sum m lam_c + k mu; -sum_c s_lc pi_c + sum_{c in C_l} lam_c + mu >= 0;
t pi_c + (m-1) lam_c >= p_c; pi,lam,mu >= 0."""
m, n = len(s), len(p)
L = [[l for l in R(m) if s[l][c] >= b] for c in R(n)]
C = [[c for c in R(n) if l in L[c]] for l in R(m)]
dl = nuovo_modello("duale_copertura")
pi = dl.addVars(n, name="pi")
lam = dl.addVars(n, name="lam")
mu = dl.addVar(name="mu")
dl.setObjective(m * lam.sum() + k * mu, GRB.MINIMIZE)
dl.addConstrs((-gp.quicksum(s[l][c] * pi[c] for c in R(n)) + gp.quicksum(lam[c] for c in C[l]) + mu >= 0
for l in R(m)), name="rc_x")
dl.addConstrs((t * pi[c] + (m - 1) * lam[c] >= p[c] for c in R(n)), name="rc_y")
return dl
m3, x3, y3, L3m = modello_3(s3, p3, t3, b3, k3)
salva_modello(m3, "fam08_3_primale")
# ---------- 2. IL RILASSAMENTO LP ----------
zlp3, zlp3r, _ = rilassamenti(m3)
# ---------- 3. IL DUALE DEL RILASSAMENTO (UPPER BOUND: E' UN MASSIMO) ----------
d3 = duale_3(s3, p3, t3, b3, k3)
salva_modello(d3, "fam08_3_duale")
mano = {"mu": 0.0}
mano.update({f"pi[{c}]": 0.0 for c in R(n)})
mano.update({f"lam[{c}]": p3[c] / 2 for c in R(n)})
ub3, viol = valuta(d3, mano)
assert viol <= 1e-9, viol
print("Soluzione duale a mano: pi = 0, mu = 0, lam_c = p_c/2 = "
+ ", ".join(frazione(p3[c] / 2) for c in R(n)) + f" -> ub = {frazione(ub3)}")
dualita_forte(d3, zlp3)
# ---------- 4. EURISTICA COSTRUTTIVA (LOWER BOUND: E' UN MASSIMO) ----------
print("Euristica: si aprono le prime k sedi; un cliente è coperto se il segnale totale")
print("raggiunge la soglia e al più una sede forte lo raggiunge.")
def euristica_3(s, p, t, b, k):
m, n = len(s), len(p)
x = [1 if l < k else 0 for l in R(m)]
y, passi = [0] * n, []
for c in R(n):
ts = sum(s[l][c] for l in R(k))
ni = sum(1 for l in R(k) if s[l][c] >= b)
y[c] = 1 if (ts >= t and ni <= 1) else 0
passi.append(f"Cliente {c + 1}: segnale totale = {ts}, sedi forti = {ni}; "
f"{'coperto' if y[c] else 'non coperto'}.")
return x, y, passi
xe, ye, passi = euristica_3(s3, p3, t3, b3, k3)
print(f" Si aprono le prime k = {k3} sedi: x = {xe}.")
for i, s in enumerate(passi, 1):
print(f" Passo {i}. {s}")
lb3 = sum(p3[c] * ye[c] for c in R(n))
print(f" lb = {lb3}")
# ---------- 5. SOLUZIONE OTTIMA DEL MILP ----------
z3 = risolvi(m3)
print("Soluzione ottima del MILP:")
stampa_soluzione(m3, solo_non_nulle=True)
riga = registra_bound("3 copertura", ub3, lb3, zlp3, zlp3r, z3, senso="max")
salva_dati(pd.DataFrame([riga]), "fam08_3_bound")
# ---------- 6. DOMANDE DI MODELLAZIONE AGGIUNTIVE ----------
varianti = {}
def variante(nome, mod):
z = risolvi(mod)
print(f" {nome:70s} z = {frazione(z)}")
return z
# 3a: almeno 3 clienti devono essere coperti
mod, x, y, L = modello_3(s3, p3, t3, b3, k3)
mod.addConstr(y.sum() >= 3, name="copertura_minima")
varianti["3a"] = variante("3a. Almeno 3 clienti coperti (sum y_c >= 3)", mod)
# 3b: se si apre la sede 1 si apre anche la sede 3
mod, x, y, L = modello_3(s3, p3, t3, b3, k3)
mod.addConstr(x[0] <= x[2], name="1_implica_3")
varianti["3b"] = variante("3b. Se si apre la sede 1 si apre anche la 3 (x_1 <= x_3)", mod)
salva_dati(pd.DataFrame({"variante": list(varianti), "z": list(varianti.values())}), "fam08_3_varianti")
# ---------- 7. IL SANDWICH SULLA VARIANTE 3a ----------
intestazione("3a. Il sandwich sulla variante: almeno 3 clienti coperti")
def modello_3a(s, p, t, b, k, minimo=3):
mod_, xx, yy, LL = modello_3(s, p, t, b, k)
mod_.addConstr(yy.sum() >= minimo, name="copertura_minima")
return mod_, xx, yy, LL
def duale_3a(s, p, t, b, k, minimo=3):
"""Al duale di 8.3 si aggiunge omega <= 0 per il vincolo sum_c y_c >= minimo
(verso >= in un massimo): entra nell'obiettivo con il suo termine noto e
nelle colonne delle y_c."""
mm, nn = len(s), len(p)
L = [[l for l in R(mm) if s[l][c] >= b] for c in R(nn)]
C = [[c for c in R(nn) if l in L[c]] for l in R(mm)]
dl = nuovo_modello("duale_copertura_3a")
pi = dl.addVars(nn, name="pi")
lam = dl.addVars(nn, name="lam")
mu = dl.addVar(name="mu")
om = dl.addVar(lb=-GRB.INFINITY, ub=0.0, name="omega")
dl.setObjective(mm * lam.sum() + k * mu + minimo * om, GRB.MINIMIZE)
dl.addConstrs((-gp.quicksum(s[l][c] * pi[c] for c in R(nn))
+ gp.quicksum(lam[c] for c in C[l]) + mu >= 0 for l in R(mm)), name="rc_x")
dl.addConstrs((t * pi[c] + (mm - 1) * lam[c] + om >= p[c] for c in R(nn)), name="rc_y")
return dl
m3a, x3a, y3a, L3a = modello_3a(s3, p3, t3, b3, k3)
salva_modello(m3a, "fam08_3a_primale")
# -- euristica ammissibile: si prova ogni scelta di k sedi e si tiene la migliore --
print("Euristica costruttiva: le sedi sono poche, quindi si prova ogni scelta di k sedi e")
print("si tiene quella che copre almeno 3 clienti con il profitto piu' alto.")
migliore_3a = None
for scelta in combinazioni(R(m), k3):
cop = []
for c in R(n):
segnale = sum(s3[l][c] for l in scelta)
forti = sum(1 for l in scelta if s3[l][c] >= b3)
if segnale >= t3 and forti <= 1:
cop.append(c)
valore = sum(p3[c] for c in cop)
print(f" sedi {[l + 1 for l in scelta]}: clienti coperti {[c + 1 for c in cop]}, "
f"profitto {valore}" + ("" if len(cop) >= 3 else " (meno di 3: non ammissibile)"))
if len(cop) >= 3 and (migliore_3a is None or valore > migliore_3a[0]):
migliore_3a = (valore, scelta, cop)
lb3a, sedi_3a, cop_3a = migliore_3a
sol_3a = ({f"x[{l}]": (1 if l in sedi_3a else 0) for l in R(m)}
| {f"y[{c}]": (1 if c in cop_3a else 0) for c in R(n)})
assert ammissibile(m3a, sol_3a), "la soluzione euristica della variante deve essere ammissibile"
print(f" la migliore e' {[l + 1 for l in sedi_3a]} -> lb = {frazione(lb3a)}")
# -- certificato duale: omega resta a zero, e si vede dal conto --
d3a = duale_3a(s3, p3, t3, b3, k3)
salva_modello(d3a, "fam08_3a_duale")
mano_3a = {"mu": 0.0, "omega": 0.0}
mano_3a.update({f"pi[{c}]": 0.0 for c in R(n)})
mano_3a.update({f"lam[{c}]": p3[c] / (m - 1) for c in R(n)})
ub3a, viol_3a = valuta(d3a, mano_3a)
assert viol_3a <= 1e-9, viol_3a
print("Soluzione duale a mano: pi = 0, mu = 0 e lam_c = p_c/(m-1) come nel problema base.")
print(" Il nuovo omega conviene lasciarlo a zero: abbassarlo costringe lam_c a salire di")
print(" -omega/(m-1) per ciascuno degli n clienti, e nell'obiettivo quelle m lam pesano")
print(" m*n/(m-1) volte piu' di quanto il termine noto faccia risparmiare.")
print(f" -> ub = {frazione(ub3a)}")
zlp3a, zlp3ar, _ = due_rilassamenti(m3a, d3a)
z3a = risolvi(m3a)
riga_3a = registra_bound("3a almeno 3 clienti coperti", ub3a, lb3a, zlp3a, zlp3ar, z3a, senso="max")
salva_dati(pd.DataFrame([riga_3a]), "fam08_3a_bound")
assert lb3a <= z3a <= zlp3a + 1e-9 <= ub3a + 1e-9
# ---------- 8. FIGURE ----------
fig, ax = plt.subplots(figsize=(7.2, 3.2))
ott_x = [l for l in R(m) if x3[l].X > 0.5]
larghezza = 0.6
for c in R(n):
colore = "#1E8449" if y3[c].X > 0.5 else "#C0392B"
ax.bar(c, p3[c], color=colore, width=larghezza)
ax.text(c, p3[c] + 0.5, "coperto" if y3[c].X > 0.5 else "non coperto", ha="center", fontsize=8)
ax.set_xticks(R(n))
ax.set_xticklabels([f"cliente {c + 1}" for c in R(n)])
ax.set_ylabel("profitto $p_c$")
ax.set_title(f"Copertura: soluzione ottima con sedi aperte {[l + 1 for l in ott_x]} (z = {frazione(z3)})")
salva_figura(fig, "cap08_copertura_ottimo")
print("Fine.")