Veicoli: lotto minimo e premio per la varietà
Classe: MILP · Legami: lotto minimo (semi-intera), contare i tipi, se e solo se · Script: python/fam09_3_veicoli.py
Difficoltà: ★★★★☆ · Tempo: 45–60 min
Problema 9.3
Un'azienda produce \(s \in \mathbb{Z}_{\ge 1}\) tipi di veicolo usando \(k \in \mathbb{Z}_{\ge 1}\) risorse. Per ogni risorsa \(i \in \{1, 2, \dots, k\}\) e ogni tipo \(j \in \{1, 2, \dots, s\}\), il valore \(a_{ij} \in \mathbb{Q}_{\ge 0}\) è la quantità di risorsa \(i\) necessaria per una unità del tipo \(j\), e \(b_i \in \mathbb{Q}_{>0}\) è la disponibilità della risorsa \(i\). Per ogni tipo \(j\), il valore \(\bar p_j \in \mathbb{Q}_{>0}\) è il profitto di una unità e \(\bar q_j \in \mathbb{Z}_{\ge 1}\) la quantità minima da produrre se si decide di produrre quel tipo. Se la produzione comprende almeno due tipi diversi, l'azienda incassa un premio \(\bar r \in \mathbb{Q}_{>0}\) (un contributo per la diversificazione). L'azienda vuole massimizzare il profitto totale.
Il problema a parole. Decidiamo quante unità produrre di ciascun tipo. L'obiettivo: profitto massimo, premio compreso. I vincoli: le risorse non si superano; e di un tipo si producono zero unità oppure almeno \(\bar q_j\).
Modello
Variabili. \(x_j \in \mathbb{Z}_{\ge 0}\) unità prodotte del tipo \(j\); \(y_j \in \{0,1\}\) vale \(1\) se il tipo \(j\) viene prodotto; \(z \in \{0,1\}\) vale \(1\) se si incassa il premio. Il dato \(M_j = \min \{ \lfloor b_i / a_{ij} \rfloor : a_{ij} > 0 \}\) è il massimo producibile del solo tipo \(j\): il minimo si prende sulle sole risorse che il tipo consuma, perché una risorsa con \(a_{ij} = 0\) non lo limita.
Descrizione. L'obiettivo somma i profitti dei veicoli prodotti e il premio per la varietà. I vincoli di risorsa, uno per risorsa, sono le disponibilità. I due vincoli di lotto minimo e di attivazione, uno per tipo ciascuno, rendono \(x_j\) semi-intera: o zero, o almeno \(\bar q_j\) e al più \(M_j\). Il vincolo di premio, uno solo, dice che il premio si incassa solo se i tipi attivi sono almeno due.
Il legame fra le variabili: la semi-interezza
I due vincoli insieme dicono
Se \(y_j = 0\) entrambi danno \(x_j = 0\); se \(y_j = 1\) si ottiene \(\bar q_j \le x_j \le M_j\). Nessuno dei due basta da solo: senza l'attivazione un tipo potrebbe essere prodotto con \(y_j = 0\); senza il lotto minimo la soglia sarebbe vuota.
| \(y_j\) | che cosa impongono i due vincoli | \(x_j\) |
|---|---|---|
| \(0\) | \(0 \le x_j \le 0\) | \(x_j = 0\) |
| \(1\) | \(\bar q_j \le x_j \le M_j\) | lotto ammesso |
Il premio si incassa solo con almeno due tipi
Il vincolo si legge \(2 z \le \sum_{j=1}^{s} y_j\). Se \(z = 1\) allora \(\sum_j y_j \ge 2\): almeno due tipi sono attivati e, per la semi-interezza, effettivamente prodotti. Il verso opposto — se ci sono due tipi attivi allora \(z = 1\) — non è imposto da alcun vincolo, ma segue dall'ottimalità: porre \(z = 1\) resta ammissibile e aumenta l'obiettivo di \(\bar r > 0\). Con \(\bar r = 0\) l'argomento cade, e \(z\) smette di essere un indicatore fedele (domanda 9.3.2).
Il modello in gurobipy
mm = gp.Model("veicoli")
x = mm.addVars(n, vtype=GRB.INTEGER, name="x")
y = mm.addVars(n, vtype=GRB.BINARY, name="y")
z = mm.addVar(vtype=GRB.BINARY, name="z")
mm.setObjective(gp.quicksum(p[j] * x[j] for j in range(n)) + r * z, GRB.MAXIMIZE)
mm.addConstrs((gp.quicksum(a[i][j] * x[j] for j in range(n)) <= b[i]
for i in range(m)), name="risorsa")
mm.addConstrs((x[j] - q[j] * y[j] >= 0 for j in range(n)), name="lotto_minimo")
mm.addConstrs((x[j] - M[j] * y[j] <= 0 for j in range(n)), name="attiva")
mm.addConstr(-gp.quicksum(y[j] for j in range(n)) + 2 * z <= 0, name="premio")
L'istanza
\(s = 3\) tipi, \(k = 2\) risorse (acciaio in tonnellate e ore di manodopera), \(\bar q_j = 10\) per ogni tipo, \(\bar r = 500\).
| \(a_{ij}\) | \(j=1\) | \(j=2\) | \(j=3\) | \(b_i\) | |
|---|---|---|---|---|---|
| \(i=1\) (acciaio) | 2 | 3 | 5 | 100 | |
| \(i=2\) (ore) | 30 | 25 | 40 | 1200 |
| \(j=1\) | \(j=2\) | \(j=3\) | |
|---|---|---|---|
| \(\bar p_j\) | 200 | 250 | 300 |
| \(\bar q_j\) | 10 | 10 | 10 |
| \(M_j\) | 40 | 33 | 20 |
I big-M sono calcolati dai dati: \(M_1 = \min(\lfloor 100/2 \rfloor, \lfloor 1200/30 \rfloor) = 40\), e analogamente \(M_2 = 33\), \(M_3 = 20\).
Il modello scritto sui dati dell'istanza:
Euristica costruttiva: il bound primale
Il problema è di massimo, quindi l'euristica dà il bound primale, che sta sotto l'ottimo. Si ordinano i tipi per profitto rapportato al consumo della risorsa più stretta, si attivano i due migliori al lotto minimo (così il premio è garantito) e poi si riempie con il più redditizio fra quelli attivati.
Sull'istanza si attivano i tipi \(2\) e \(1\) al lotto minimo, poi si riempie: la produzione è \((11, 26, 0)\), l'acciaio si esaurisce e restano \(220\) ore. Il profitto è \(8700\) più il premio di \(500\):
Rilassamento LP e duale: il bound duale
Il primale è di massimo con vincoli \(\le\) e \(\ge\): seguendo la tabella di conversione, si associano \(\pi_i \ge 0\) alle risorse, \(\ell_j \ge 0\) al lotto minimo (scritto come \(-\lambda_j\)), \(\beta_j \ge 0\) all'attivazione e \(\gamma \ge 0\) al premio.
Lo stesso duale, scritto sui dati dell'istanza:
Descrizione. \(\pi_i\) è il prezzo di una unità della risorsa \(i\); \(\ell_j\) e \(\beta_j\) sono i prezzi dei due vincoli di semi-interezza del tipo \(j\), e \(\gamma\) quello del premio. L'obiettivo valuta a quei prezzi tutte le risorse disponibili. Il primo gruppo sono le colonne delle \(x_j\): le risorse che una unità del tipo \(j\) consuma, corrette dai due vincoli di semi-interezza, devono coprire il profitto \(\bar p_j\). Il secondo sono le colonne delle \(y_j\): accendere il tipo \(j\) obbliga a produrne almeno \(\bar q_j\) e ne concede al più \(M_j\), e il saldo deve coprire il premio \(\gamma\). L'ultimo è la colonna di \(z\): il premio si incassa solo con due tipi attivi, e infatti va coperto da \(2\gamma\).
Ricetta, in tre passi.
- Il vincolo su \(z\) obbliga \(\gamma \ge \bar r/2\): si prende il minimo, \(\bar\gamma = 250\).
- Si pone \(\bar\beta_j = 0\) e si ricava il più piccolo \(\ell_j\) ammissibile, \(\bar\ell_j = \bar\gamma / \bar q_j = 25\) per ogni \(j\): ogni tipo attivato «porta con sé» la sua quota di premio.
- Resta \(\sum_i a_{ij}\, \pi_i \ge \bar p_j + \bar\ell_j\). Si valuta una sola risorsa, al prezzo che copre tutti i tipi, e si tiene quella che dà il bound più piccolo.
Il bound migliore è quello dell'acciaio: \(z(\mathit{MILP}) \le \mathit{UB} = 11\,250\).
Soluzione ottima
La produzione ottima è \((26, 16, 0)\): si attivano i tipi \(1\) e \(2\), si incassa il premio, si consumano tutte e \(100\) le tonnellate di acciaio e \(1180\) ore sulle \(1200\) disponibili.
| \(LB\) (euristica) | \(z(\mathit{MILP})\) | \(z(\mathit{LP}^+)\) | \(z(\mathit{LP})\) | \(UB\) (duale) | gap dell'euristica |
|---|---|---|---|---|---|
| 9200 | 9700 | 9750 | \(20625/2\) | 11250 | \(5{,}2\%\) |
Qui il rilassamento con i bound batte il duale a mano
È l'unico problema del capitolo in cui \(z(\mathit{LP}^+) = 9750\) è migliore del bound duale costruito a mano (\(11\,250\)), e di pochissimo peggiore dell'ottimo (\(9700\)). La ragione è che il rilassamento senza i bound lascia \(y_j\) e \(z\) crescere sopra \(1\), e con esse il premio: \(z(\mathit{LP}) = 20625/2 \approx 10\,312\). Aggiungere \(y_j \le 1\) e \(z \le 1\) toglie proprio quella libertà. Quando il modello contiene indicatori premiati, il rilassamento con i bound non è un dettaglio.
Considerazioni aggiuntive
- Il big-M \(M_j\) è il massimo producibile del solo tipo \(j\), non del piano complessivo: è già molto più stretto di una costante arbitraria.
- Se per qualche tipo fosse \(\bar q_j > M_j\), quel tipo sarebbe impossibile da produrre e si potrebbe eliminare dal modello con un controllo sui dati.
- Il premio è modellato con una binaria e un vincolo. Con una soglia diversa (almeno \(f\) tipi) basta sostituire il coefficiente \(2\) con \(f\); con più premi a scaglioni servirebbero una variabile per scaglione e una catena di vincoli, come nelle funzioni a tratti.
Domande di modellazione aggiuntive
9.3.1 — Premio nullo
Il contributo per la diversificazione viene abolito, cioè \(\bar r = 0\). Che cosa succede alla variabile \(z\)?
Una variante svolta: premio più esigente
Il premio si incassa solo se si producono almeno tre tipi diversi.
Basta cambiare un coefficiente: il vincolo diventa
Sull'istanza il nuovo ottimo è \(9200\), cioè \(500\) in meno. Attivare anche il terzo tipo costerebbe più del premio: conviene rinunciare al premio, e infatti all'ottimo \(z = 0\). Il valore coincide, per caso, con quello della soluzione euristica del problema di partenza.
Il vincolo del premio tocca una sola colonna del duale, quella di \(z\), che adesso vede \(\bar r / 3\) invece di \(\bar r / 2\); tutto il resto del duale è identico. La ricetta è quella del problema base con \(\gamma = \bar r / 3 = 500/3\), e il bound scende da \(11250\) a \(32500/3\): chiedere tre tipi invece di due spalma il premio su più attivazioni, quindi ogni tipo ne porta una quota minore e i prezzi duali si abbassano. L'euristica parte dal piano base, prova ad accendere un terzo tipo al suo lotto minimo e con queste risorse non ci riesce: resta a \(8700\), e l'ottimo scende da \(9700\) a \(9200\).
| valore | che cos'è | |
|---|---|---|
| \(\mathit{UB}\) | \(\frac{32500}{3}\) | certificato duale costruito a mano |
| \(\mathit{LB}\) | \(8700\) | soluzione euristica |
| \(z(\mathit{LP})\) | \(\frac{29875}{3}\) | rilassamento senza i bound |
| \(z(\mathit{LP}^+)\) | \(\frac{28750}{3}\) | rilassamento con i bound |
| \(z(\mathit{MILP})\) | \(9200\) | ottimo del MILP |
Codice
Script completo —
python/fam09_3_veicoli.py
(riproducibile con python3 python/fam09_3_veicoli.py dalla cartella
python/). Notebook —
notebooks/fam09_3_veicoli.ipynb
— che si apre in Colab dal badge in cima alla pagina.
Mostra lo script completo — python/fam09_3_veicoli.py (287 righe)
"""Problema 9.3 -- Veicoli: lotto minimo e premio per la varieta'.
Tre tecniche insieme: la variabile semi-intera del lotto minimo (3.3), il
conteggio dei tipi attivi (3.11) e un premio «se e solo se» si producono almeno
due tipi (3.10). Il premio si incassa solo se il conteggio arriva a due: il
verso mancante segue dall'ottimalita' perche' il premio e' positivo.
"""
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, valuta)
from stile import ARANCIO, BLU, ROSSO, TEAL, VERDE, intestazione, plt, salva_dati, salva_figura
from esteso import salva_modello
R = range
# ---------- 1. MODELLO E ISTANZA ----------
intestazione("9.3 Veicoli: lotto minimo per tipo e premio se si producono almeno due tipi")
a3 = [[2, 3, 5], # acciaio (tonnellate) per unita' dei tre tipi
[30, 25, 40]] # ore di manodopera per unita'
b3 = [100, 1200] # acciaio e ore disponibili
p3 = [200, 250, 300] # profitto per unita'
q3 = [10, 10, 10] # quantita' minima se il tipo si produce
r3 = 500 # premio se si producono almeno due tipi
n3, m3 = 3, 2
# il piu' piccolo big-M valido per tipo: quante unita' al massimo consentono le risorse
M3 = [min(b3[i] // a3[i][j] for i in R(m3)) for j in R(n3)]
salva_dati(pd.DataFrame({"tipo": R(1, n3 + 1), "acciaio": a3[0], "ore": a3[1],
"profitto": p3, "minimo": q3, "M": M3}), "fam09_3_dati")
print(f" Risorse: {b3[0]} t di acciaio, {b3[1]} ore. Big-M per tipo (dai soli dati): {M3}")
def modello_3(a, b, p, q, r):
n, m = len(p), len(b)
M = [min(b[i] // a[i][j] for i in R(m)) for j in R(n)]
mm = nuovo_modello("veicoli")
x = mm.addVars(n, vtype=GRB.INTEGER, name="x") # unita' prodotte
y = mm.addVars(n, vtype=GRB.BINARY, name="y") # tipo attivato
z = mm.addVar(vtype=GRB.BINARY, name="z") # premio per la varieta'
mm.setObjective(gp.quicksum(p[j] * x[j] for j in R(n)) + r * z, GRB.MAXIMIZE)
mm.addConstrs((gp.quicksum(a[i][j] * x[j] for j in R(n)) <= b[i] for i in R(m)),
name="risorsa")
mm.addConstrs((x[j] - q[j] * y[j] >= 0 for j in R(n)), name="lotto_minimo")
mm.addConstrs((x[j] - M[j] * y[j] <= 0 for j in R(n)), name="attiva")
mm.addConstr(-gp.quicksum(y[j] for j in R(n)) + 2 * z <= 0, name="premio")
return mm, x, y, z
def duale_3(a, b, p, q, r):
"""min sum_i b_i pi_i; sum_i a_ij pi_i - alpha_j + beta_j >= p_j;
q_j alpha_j - M_j beta_j + gamma >= 0; -2 gamma >= r; pi, alpha, beta >= 0, gamma <= 0.
(scritto con i segni della tabella di conversione per un primale di massimo)"""
n, m = len(p), len(b)
M = [min(b[i] // a[i][j] for i in R(m)) for j in R(n)]
dl = nuovo_modello("duale_veicoli")
pi = dl.addVars(m, name="pi") # risorse (<= in un max)
alpha = dl.addVars(n, lb=-GRB.INFINITY, ub=0.0, name="alpha") # lotto minimo (>= in un max)
beta = dl.addVars(n, name="beta") # attivazione (<=)
gamma = dl.addVar(name="gamma") # premio (<=)
dl.setObjective(gp.quicksum(b[i] * pi[i] for i in R(m)), GRB.MINIMIZE)
dl.addConstrs((gp.quicksum(a[i][j] * pi[i] for i in R(m)) + alpha[j] + beta[j] >= p[j]
for j in R(n)), name="rc_x")
dl.addConstrs((-q[j] * alpha[j] - M[j] * beta[j] - gamma >= 0 for j in R(n)), name="rc_y")
dl.addConstr(2 * gamma >= r, name="rc_z")
return dl
m3m, x3, y3, z3 = modello_3(a3, b3, p3, q3, r3)
salva_modello(m3m, "fam09_3_primale")
# ---------- 2. IL RILASSAMENTO LP ----------
zlp3, zlp3r, _ = rilassamenti(m3m)
# ---------- 3. IL DUALE DEL RILASSAMENTO (UPPER BOUND: E' UN MASSIMO) ----------
dl3 = duale_3(a3, b3, p3, q3, r3)
salva_modello(dl3, "fam09_3_duale")
# ricetta: gamma = r/2 (il minimo ammesso dal vincolo 2 gamma >= r), beta = 0, e
# lambda_j = gamma / q_j (ogni tipo attivato "porta" la sua quota di premio); poi si
# valuta una sola risorsa al prezzo che copre tutti i tipi, e si sceglie la migliore
gamma = r3 / 2
lam = [gamma / q3[j] for j in R(n3)]
bound = {}
for i in R(m3):
prezzo = max((p3[j] + lam[j]) / a3[i][j] for j in R(n3))
bound[i] = b3[i] * prezzo
critica = min(bound, key=bound.get)
prezzo = max((p3[j] + lam[j]) / a3[critica][j] for j in R(n3))
mano = {"gamma": gamma} | {f"pi[{i}]": 0.0 for i in R(m3)} \
| {f"alpha[{j}]": -lam[j] for j in R(n3)} | {f"beta[{j}]": 0.0 for j in R(n3)}
mano[f"pi[{critica}]"] = prezzo
ub3, viol = valuta(dl3, mano)
assert viol <= 1e-9, (viol, mano)
print(f" Duale a mano: gamma = r/2 = {frazione(gamma)} (il minimo che soddisfa 2 gamma >= r),")
print(f" beta = 0 e lambda_j = gamma / q_j = " + ", ".join(frazione(v) for v in lam)
+ ": ogni tipo")
print(" attivato porta la sua quota di premio. Poi si valuta una sola risorsa al prezzo")
print(" che copre tutti i tipi, max_j (p_j + lambda_j) / a_ij, e si tiene la piu' stretta:")
for i in R(m3):
print(f" risorsa {i + 1}: prezzo {frazione(max((p3[j] + lam[j]) / a3[i][j] for j in R(n3)))}"
f" -> b_i * prezzo = {frazione(bound[i])}")
print(f" Il minimo e' la risorsa {critica + 1}: ub = {frazione(ub3)}")
dualita_forte(dl3, zlp3)
# ---------- 4. EURISTICA COSTRUTTIVA (LOWER BOUND: E' UN MASSIMO) ----------
# euristica costruttiva: si attivano due tipi (per incassare il premio) partendo dai profitti per
# unita' di risorsa piu' scarsa, poi si riempie con il tipo migliore
def euristica(a, b, p, q, r):
n, m = len(p), len(b)
# rapporto profitto / consumo della risorsa piu' stretta
ordine = sorted(R(n), key=lambda j: -p[j] / max(a[i][j] / b[i] for i in R(m)))
x = [0] * n
res = list(b)
attivi = []
for j in ordine: # prima il lotto minimo dei due tipi migliori
if len(attivi) < 2 and all(res[i] >= a[i][j] * q[j] for i in R(m)):
x[j] = q[j]
for i in R(m):
res[i] -= a[i][j] * q[j]
attivi.append(j)
for j in ordine: # poi si riempie con il tipo piu' redditizio
if x[j] == 0:
continue
extra = min(res[i] // a[i][j] for i in R(m))
x[j] += extra
for i in R(m):
res[i] -= a[i][j] * extra
return x, attivi, res
x_eur, attivi, res = euristica(a3, b3, p3, q3, r3)
lb3 = sum(p3[j] * x_eur[j] for j in R(n3)) + (r3 if len(attivi) >= 2 else 0)
sol_eur = {f"x[{j}]": x_eur[j] for j in R(n3)} \
| {f"y[{j}]": 1 if x_eur[j] > 0 else 0 for j in R(n3)} | {"z": 1 if len(attivi) >= 2 else 0}
assert ammissibile(m3m, sol_eur)
print(f" Euristica: si attivano i tipi {[j + 1 for j in attivi]} al lotto minimo, poi si")
print(f" riempie col piu' redditizio; produzione {x_eur}, risorse residue {res}")
print(f" lb = {sum(p3[j] * x_eur[j] for j in R(n3))} + {r3} di premio = {frazione(lb3)}")
# ---------- 5. OTTIMO DEL MILP ----------
z3v = risolvi(m3m)
print(" Soluzione ottima: produzione " + ", ".join(str(round(x3[j].X)) for j in R(n3))
+ f"; tipi attivi {[j + 1 for j in R(n3) if y3[j].X > 0.5]}; premio incassato: "
+ ("si" if z3.X > 0.5 else "no"))
print(" Risorse usate: " + ", ".join(
f"{frazione(sum(a3[i][j] * round(x3[j].X) for j in R(n3)))} su {b3[i]}" for i in R(m3)))
riga = registra_bound("3 veicoli", ub3, lb3, zlp3, zlp3r, z3v, senso="max")
salva_dati(pd.DataFrame([riga]), "fam09_3_bound")
assert lb3 <= z3v <= zlp3 + 1e-6 <= ub3 + 1e-6
# ---------- 6. DOMANDE DI MODELLAZIONE AGGIUNTIVE ----------
varianti = {}
def variante(nome, m):
z = risolvi(m)
print(f" {nome:70s} z = {frazione(z)}")
return z
# 3a: il premio richiede almeno tre tipi diversi
m, x, y, z = modello_3(a3, b3, p3, q3, r3)
m.update()
m.remove([c for c in m.getConstrs() if c.ConstrName == "premio"])
m.addConstr(-gp.quicksum(y[j] for j in R(n3)) + 3 * z <= 0, name="premio3")
varianti["3a"] = variante("3a. Il premio si incassa solo con almeno tre tipi diversi", m)
# 3b: il premio e' nullo -- che cosa succede al legame "se e solo se"?
m, x, y, z = modello_3(a3, b3, p3, q3, 0)
zz = risolvi(m)
print(f" {'3b. Il premio vale 0: z non e piu un indicatore fedele':70s} z = {frazione(zz)}")
print(f" tipi attivi {[j + 1 for j in R(n3) if y[j].X > 0.5]}, ma z = {round(z.X)}: con")
print(" premio nullo l'ottimo non ha alcun motivo di alzare z, e il vincolo da solo")
print(" non lo impone. Per farne un indicatore fedele serve anche il verso opposto.")
varianti["3b"] = zz
salva_dati(pd.DataFrame({"variante": list(varianti), "z": list(varianti.values())}),
"fam09_3_varianti")
# ---------- 7. IL SANDWICH SULLA VARIANTE 3a ----------
intestazione("9.3a Il sandwich sulla variante: il premio richiede almeno tre tipi")
SOGLIA = 3
def modello_3a(a, b, p, q, r, soglia=SOGLIA):
mm, xx, yy, zz = modello_3(a, b, p, q, r)
mm.update() # i nomi dei vincoli esistono dopo l'update
mm.remove(mm.getConstrByName("premio"))
mm.addConstr(-gp.quicksum(yy[j] for j in R(len(p))) + soglia * zz <= 0, name="premio")
return mm, xx, yy, zz
def duale_3a(a, b, p, q, r, soglia=SOGLIA):
"""Rispetto al duale di 9.3 cambia una sola riga: la colonna di z, che
adesso vede il premio diviso per `soglia` invece che per due. Il resto del
duale e' identico --- il vincolo del premio tocca solo quella colonna."""
nn, mm_ = len(p), len(b)
MM = [min(b[i] // a[i][j] for i in R(mm_)) for j in R(nn)]
dl = nuovo_modello("duale_veicoli_3a")
pi = dl.addVars(mm_, name="pi")
alpha = dl.addVars(nn, lb=-GRB.INFINITY, ub=0.0, name="alpha")
beta = dl.addVars(nn, name="beta")
gm = dl.addVar(name="gamma")
dl.setObjective(gp.quicksum(b[i] * pi[i] for i in R(mm_)), GRB.MINIMIZE)
dl.addConstrs((gp.quicksum(a[i][j] * pi[i] for i in R(mm_)) + alpha[j] + beta[j] >= p[j]
for j in R(nn)), name="rc_x")
dl.addConstrs((-q[j] * alpha[j] - MM[j] * beta[j] - gm >= 0 for j in R(nn)), name="rc_y")
dl.addConstr(soglia * gm >= r, name="rc_z")
return dl
m3a, x3a, y3a, z3a = modello_3a(a3, b3, p3, q3, r3)
salva_modello(m3a, "fam09_3a_primale")
# -- euristica ammissibile: la base, poi si accende il terzo tipo se il premio conviene --
print("Euristica costruttiva: si parte dalla soluzione del problema base; se i tipi attivi")
print("sono meno di tre, si prova ad accenderne un terzo al suo lotto minimo, togliendo")
print("unita' al tipo meno redditizio finche' le risorse bastano. Si tiene il piano migliore.")
x_3a = list(x_eur)
attivi_3a = [j for j in R(n3) if x_3a[j] > 0]
base_3a = sum(p3[j] * x_3a[j] for j in R(n3)) + (r3 if len(attivi_3a) >= SOGLIA else 0)
print(f" soluzione base: tipi attivi {[j + 1 for j in attivi_3a]}, valore {frazione(base_3a)}")
migliore_3a = (base_3a, list(x_3a))
if len(attivi_3a) < SOGLIA:
for spento in (j for j in R(n3) if x_3a[j] == 0):
prova = list(x_3a)
prova[spento] = q3[spento]
peggiore = min(attivi_3a, key=lambda j: p3[j])
while any(sum(a3[i][j] * prova[j] for j in R(n3)) > b3[i] for i in R(m3)):
if prova[peggiore] == 0:
break
prova[peggiore] -= 1
if all(sum(a3[i][j] * prova[j] for j in R(n3)) <= b3[i] for i in R(m3)) \
and sum(1 for j in R(n3) if prova[j] > 0) >= SOGLIA:
valore = sum(p3[j] * prova[j] for j in R(n3)) + r3
print(f" accendendo il tipo {spento + 1} al minimo {q3[spento]}: {frazione(valore)}")
if valore > migliore_3a[0]:
migliore_3a = (valore, prova)
lb3a, x_mig = migliore_3a
att_mig = [j for j in R(n3) if x_mig[j] > 0]
sol_3a = ({f"x[{j}]": x_mig[j] for j in R(n3)}
| {f"y[{j}]": (1 if x_mig[j] > 0 else 0) for j in R(n3)}
| {"z": 1 if len(att_mig) >= SOGLIA else 0})
assert ammissibile(m3a, sol_3a), "la soluzione euristica della variante deve essere ammissibile"
print(f" piano scelto {x_mig} -> lb = {frazione(lb3a)}")
# -- certificato duale: il premio si divide per tre invece che per due --
dl3a = duale_3a(a3, b3, p3, q3, r3)
salva_modello(dl3a, "fam09_3a_duale")
gamma_3a = r3 / SOGLIA
lam_3a = [gamma_3a / q3[j] for j in R(n3)]
bound_3a = {i: b3[i] * max((p3[j] + lam_3a[j]) / a3[i][j] for j in R(n3)) for i in R(m3)}
critica_3a = min(bound_3a, key=bound_3a.get)
prezzo_3a = max((p3[j] + lam_3a[j]) / a3[critica_3a][j] for j in R(n3))
mano_3a = ({"gamma": gamma_3a} | {f"pi[{i}]": 0.0 for i in R(m3)}
| {f"alpha[{j}]": -lam_3a[j] for j in R(n3)} | {f"beta[{j}]": 0.0 for j in R(n3)})
mano_3a[f"pi[{critica_3a}]"] = prezzo_3a
ub3a, viol_3a = valuta(dl3a, mano_3a)
assert viol_3a <= 1e-9, (viol_3a, mano_3a)
print(f"Soluzione duale a mano: la stessa ricetta del problema base, con gamma = r/{SOGLIA} =")
print(f" {frazione(gamma_3a)} invece di r/2: chiedere tre tipi invece di due divide il premio")
print(" su piu' attivazioni, quindi ogni tipo ne porta una quota minore e i prezzi scendono.")
print(f" lambda_j = gamma/q_j = {[frazione(v) for v in lam_3a]}; risorsa critica: "
f"la {critica_3a + 1}")
print(f" -> ub = {frazione(ub3a)} (con la soglia a due era {frazione(ub3)})")
zlp3a, zlp3ar, _ = due_rilassamenti(m3a, dl3a)
z3a_val = risolvi(m3a)
riga_3a = registra_bound("3a premio con tre tipi", ub3a, lb3a, zlp3a, zlp3ar, z3a_val, senso="max")
salva_dati(pd.DataFrame([riga_3a]), "fam09_3a_bound")
assert lb3a <= z3a_val <= zlp3a + 1e-9 <= ub3a + 1e-9
# ---------- 8. FIGURA ----------
fig, ax = plt.subplots(figsize=(6.8, 3.0))
tipi = list(R(1, n3 + 1))
colori = [TEAL if y3[j].X > 0.5 else "#F4F6F7" for j in R(n3)]
ax.bar(tipi, [x3[j].X for j in R(n3)], color=colori, edgecolor="#7F8C8D", width=0.55)
for j in R(n3):
ax.plot([j + 0.72, j + 1.28], [q3[j], q3[j]], color=ROSSO, lw=2)
ax.plot([], [], color=ROSSO, lw=2, label="lotto minimo $q_j$")
for j in R(n3):
ax.annotate(str(round(x3[j].X)), (j + 1, x3[j].X), ha="center", va="bottom", fontsize=9)
ax.set_xticks(tipi)
ax.set_xticklabels([f"tipo {j}" for j in tipi])
ax.set_ylabel("unita' prodotte")
ax.set_title(f"9.3: piano ottimo (z = {frazione(z3v)}, premio incassato)")
ax.legend(fontsize=8)
salva_figura(fig, "cap09_veicoli_ottimo")
print("Fine.")