Vai al contenuto

Produzione e manodopera: due formulazioni equivalenti

Classe: MILP · Legami: conteggi interi, bilancio dell'organico · Script: python/fam09_2_manodopera.py
Difficoltà: ★★★★★ · Tempo: 45–60 min

Apri in Colab

Problema 9.2

Un calzaturificio deve pianificare la produzione su \(n \in \mathbb{Z}_{\ge 1}\) mesi. Per ogni mese \(t \in \{1, 2, \dots, n\}\), il valore \(d_t \in \mathbb{Z}_{\ge 0}\) è la domanda in paia e \(p_t \in \mathbb{Q}_{>0}\) il costo delle materie prime per un paio. Per ogni mese \(t \in \{1, 2, \dots, n-1\}\), il valore \(h_t \in \mathbb{Q}_{\ge 0}\) è il costo di tenere un paio in magazzino a fine mese. All'inizio dell'orizzonte l'azienda ha \(m_0 \in \mathbb{Z}_{\ge 0}\) operai; ogni operaio lavora \(r \in \mathbb{Q}_{>0}\) ore al mese e costa \(w \in \mathbb{Q}_{>0}\) euro al mese, e produrre un paio richiede \(g \in \mathbb{Q}_{>0}\) ore di lavoro. All'inizio di ogni mese si possono assumere operai, al costo di \(u \in \mathbb{Q}_{\ge 0}\) euro ciascuno; non si licenzia nessuno. L'azienda vuole decidere quanto produrre e quanti operai avere in ogni mese, al costo totale minimo.

Il problema a parole. Decidiamo quanto produrre, quanto tenere in magazzino e quanti operai assumere. L'obiettivo: costo totale (materie prime, magazzino, salari e assunzioni) minimo. I vincoli: la domanda va soddisfatta esattamente; la produzione di un mese non può richiedere più ore di quelle che gli operai in servizio riescono a fare.

Due formulazioni

La stessa decisione si può scrivere in due modi, e vale la pena metterli uno accanto all'altro: è il primo caso del corso in cui due modelli apparentemente diversi descrivono lo stesso insieme di piani.

Formulazione A: le assunzioni. Le variabili di personale sono

\[ z_t = \text{operai assunti all'inizio del mese } t, \qquad \forall t \in \{1, 2, \dots, n\}, \]

intere e non negative. Un operaio assunto al mese \(t\) resta in servizio fino alla fine dell'orizzonte, quindi costa \(u\) una volta più \(w\) per ciascuno dei \(n - t + 1\) mesi restanti. Il salario degli \(m_0\) operai iniziali, \(m_0\, w\, n\), è un termine costante: si lascia fuori dal modello e si somma al valore finale.

Formulazione B: l'organico. Le variabili di personale sono

\[ y_t = \text{operai in servizio nel mese } t, \qquad \forall t \in \{1, 2, \dots, n\}, \]

intere e non negative, con \(y_t \ge y_{t-1}\) perché non si licenzia. Il costo è \(w\, y_t\) ogni mese, più \(u\) per ogni assunzione, cioè \(u\,(y_n - m_0)\) in totale, perché le assunzioni sono gli incrementi dell'organico e la somma telescopica lascia solo gli estremi.

Modello

Variabili. \(x_t \in \mathbb{Z}_{\ge 0}\) paia prodotte nel mese \(t\); \(s_t \in \mathbb{Z}_{\ge 0}\) scorta a fine mese \(t\) (per \(t \le n-1\)); \(y_t \in \mathbb{Z}_{\ge 0}\) operai in servizio nel mese \(t\); \(z_t \in \mathbb{Z}_{\ge 0}\) operai assunti all'inizio del mese \(t\).

Modello 9.2A — con le assunzioni.

\[ \begin{aligned} \min ~~ \sum_{t=1}^{n} p_t\, x_t + \sum_{t=1}^{n-1} h_t\, s_t + \sum_{t=1}^{n} \bigl(u + w\,(n - t + 1)\bigr) z_t & & \\ \text{soggetto a} \quad x_1 - s_1 &= d_1, & \\ x_t + s_{t-1} - s_t &= d_t, & \forall t \in \{2, 3, \dots, n-1\}, \\ x_n + s_{n-1} &= d_n, & \\ -g\, x_t + r \sum_{j=1}^{t} z_j &\ge -r\, m_0, & \forall t \in \{1, 2, \dots, n\}, \\ x_t &\in \Z_{\ge 0}, & \forall t \in \{1, 2, \dots, n\}, \\ s_t &\in \Z_{\ge 0}, & \forall t \in \{1, 2, \dots, n-1\}, \\ z_t &\in \Z_{\ge 0}, & \forall t \in \{1, 2, \dots, n\}. \end{aligned} \]

Modello 9.2B — con l'organico.

\[ \begin{aligned} \min ~~ \sum_{t=1}^{n} p_t\, x_t + \sum_{t=1}^{n-1} h_t\, s_t + \sum_{t=1}^{n} w\, y_t + u\,(y_n - m_0) & & \\ \text{soggetto a} \quad x_1 - s_1 &= d_1, & \\ x_t + s_{t-1} - s_t &= d_t, & \forall t \in \{2, 3, \dots, n-1\}, \\ x_n + s_{n-1} &= d_n, & \\ -g\, x_t + r\, y_t &\ge 0, & \forall t \in \{1, 2, \dots, n\}, \\ y_1 &\ge m_0, & \\ -y_{t-1} + y_t &\ge 0, & \forall t \in \{2, 3, \dots, n\}, \\ x_t &\in \Z_{\ge 0}, & \forall t \in \{1, 2, \dots, n\}, \\ s_t &\in \Z_{\ge 0}, & \forall t \in \{1, 2, \dots, n-1\}, \\ y_t &\in \Z_{\ge 0}, & \forall t \in \{1, 2, \dots, n\}. \end{aligned} \]

Descrizione. Le due formulazioni condividono i bilanci delle scorte, uno per mese: quanto si produce più quanto si ha in magazzino copre esattamente la domanda. I vincoli delle ore, uno per mese, dicono che la produzione del mese non può richiedere più ore di quante ne facciano gli operai in servizio. Nella formulazione \(B\) ci sono in più il vincolo di organico iniziale, che parte dagli \(m_0\) operai già assunti, e i vincoli di monotonia, uno per mese a partire dal secondo, che vietano i licenziamenti. Nella \(A\) quegli stessi fatti sono nascosti nel dominio delle \(z_t\), che sono non negative.

Le due formulazioni descrivono lo stesso problema

La corrispondenza è

\[ y_t = m_0 + \sum_{j=1}^{t} z_j, \qquad\text{cioè}\qquad z_t = y_t - y_{t-1} \quad (\text{con } y_0 = m_0). \]

Ammissibilità. Con questa sostituzione il vincolo delle ore di \(A\) diventa esattamente quello di \(B\); le condizioni \(z_t \ge 0\) diventano \(y_t \ge y_{t-1}\), e \(y_1 \ge m_0\). I bilanci non contengono variabili di personale e restano identici: la corrispondenza è una biiezione fra i piani ammissibili dei due modelli.

Costo. Il costo del personale in \(B\) vale

\[ \sum_{t=1}^{n} w\, y_t + u\,(y_n - m_0) = m_0\, w\, n + \sum_{t=1}^{n} \bigl(u + w\,(n - t + 1)\bigr) z_t , \]

perché \(z_j\) compare in tutti i mesi da \(j\) in poi, cioè \(n - j + 1\) volte. Le due funzioni obiettivo differiscono per la sola costante \(m_0\, w\, n\), e i due modelli hanno gli stessi ottimi.

Perché tenerle entrambe

La formulazione \(A\) ha meno vincoli (niente monotonia) ma coefficienti di costo che dipendono dal periodo; la \(B\) ha coefficienti uniformi e si estende meglio se si aggiungono i licenziamenti (basta una seconda famiglia \(\ell_t \ge 0\) e il bilancio \(y_t = y_{t-1} + z_t - \ell_t\)). Sull'istanza anche i rilassamenti coincidono: \(z(\mathit{LP}) = 15\,960\) per entrambe. Due formulazioni equivalenti sull'intero non lo sono sempre sul rilassamento; qui lo sono, e la verifica va fatta, non data per scontata.

Il modello in gurobipy

m = gp.Model("manodopera_B")
x = m.addVars(n, vtype=GRB.INTEGER, name="x")
s = m.addVars(n - 1, vtype=GRB.INTEGER, name="s")
y = m.addVars(n, vtype=GRB.INTEGER, name="y")
m.setObjective(gp.quicksum(p[t] * x[t] for t in range(n))
               + gp.quicksum(h[t] * s[t] for t in range(n - 1))
               + gp.quicksum(w * y[t] for t in range(n)) + u * (y[n - 1] - m0), GRB.MINIMIZE)
m.addConstr(x[0] - s[0] == d[0], name="bilancio[0]")
m.addConstrs((x[t] + s[t - 1] - s[t] == d[t] for t in range(1, n - 1)), name="bilancio")
m.addConstr(x[n - 1] + s[n - 2] == d[n - 1], name=f"bilancio[{n - 1}]")
m.addConstrs((-g * x[t] + r * y[t] >= 0 for t in range(n)), name="ore")
m.addConstr(y[0] >= m0, name="organico_iniziale")
m.addConstrs((-y[t - 1] + y[t] >= 0 for t in range(1, n)), name="mai_licenziamenti")

L'istanza

\(n = 3\) mesi, \(m_0 = 2\) operai, \(w = 1500\), \(u = 100\), \(r = 160\) ore, \(g = 4\) ore per paio, \(h_t = 3\).

\(t=1\) \(t=2\) \(t=3\)
\(d_t\) 60 100 140
\(p_t\) 15 15 15

Con due operai la capacità iniziale è \(2 \cdot 160 / 4 = 80\) paia al mese: basta per il primo mese, non per gli altri due.

Il modello scritto sui dati dell'istanza:

\[ \begin{array}{rrrrrrrrr c l} \min & 15x_1 & +15x_2 & +15x_3 & +3s_1 & +3s_2 & +4600z_1 & +3100z_2 & +1600z_3 & & \\ \text{soggetto a} & x_1 & & & -s_1 & & & & & = & 60\\ & & x_2 & & +s_1 & -s_2 & & & & = & 100\\ & & & x_3 & & +s_2 & & & & = & 140\\ & -4x_1 & & & & & +160z_1 & & & \ge & -320\\ & & -4x_2 & & & & +160z_1 & +160z_2 & & \ge & -320\\ & & & -4x_3 & & & +160z_1 & +160z_2 & +160z_3 & \ge & -320\\ & x_1, & x_2, & x_3 & & & & & & \in & \Z_{\ge 0}\\ & & & & s_1, & s_2 & & & & \in & \Z_{\ge 0}\\ & & & & & & z_1, & z_2, & z_3 & \in & \Z_{\ge 0} \end{array} \]

Euristica costruttiva: il bound primale

Produzione «just in time»: ogni mese si produce esattamente la domanda, e si assume il numero minimo di operai che serve. È un'euristica costruttiva: si costruisce una sola soluzione, un elemento per volta, senza mai tornare indietro.

  • mese 1: \(\lceil 4 \cdot 60/160 \rceil = 2\) operai, nessuna assunzione;
  • mese 2: \(\lceil 4 \cdot 100/160 \rceil = 3\) operai, una assunzione;
  • mese 3: \(\lceil 4 \cdot 140/160 \rceil = 4\) operai, un'altra assunzione.

Il costo, incluso il termine costante \(m_0\, w\, n = 9000\), è \(z(\mathit{MILP}) \le \mathit{UB} = 18\,200\).

Rilassamento LP e duale: il bound duale

Sulla formulazione \(A\), con \(\mu_t\) libera su ciascun bilancio e \(\nu_t \ge 0\) su ciascun vincolo di ore:

\[ \begin{aligned} \max ~~ \sum_{t=1}^{n} d_t\, \mu_t - r\, m_0 \sum_{t=1}^{n} \nu_t & & \\ \text{soggetto a} \quad \mu_t - g\, \nu_t &\le p_t, & \forall t \in \{1, 2, \dots, n\}, \\ -\mu_t + \mu_{t+1} &\le h_t, & \forall t \in \{1, 2, \dots, n-1\}, \\ r \sum_{t=j}^{n} \nu_t &\le u + w\,(n - j + 1), & \forall j \in \{1, 2, \dots, n\}, \\ \mu_t &\gtreqless 0, & \forall t \in \{1, 2, \dots, n\}, \\ \nu_t &\ge 0, & \forall t \in \{1, 2, \dots, n\}. \end{aligned} \]

Lo stesso duale, scritto sui dati dell'istanza:

\[ \begin{array}{rrrrrrr c l} \max & 60\mu_1 & +100\mu_2 & +140\mu_3 & -320\nu_1 & -320\nu_2 & -320\nu_3 & & \\ \text{soggetto a} & \mu_1 & & & -4\nu_1 & & & \le & 15\\ & & \mu_2 & & & -4\nu_2 & & \le & 15\\ & & & \mu_3 & & & -4\nu_3 & \le & 15\\ & -\mu_1 & +\mu_2 & & & & & \le & 3\\ & & -\mu_2 & +\mu_3 & & & & \le & 3\\ & & & & 160\nu_1 & +160\nu_2 & +160\nu_3 & \le & 4600\\ & & & & & 160\nu_2 & +160\nu_3 & \le & 3100\\ & & & & & & 160\nu_3 & \le & 1600\\ & \mu_1, & \mu_2, & \mu_3 & & & & \gtreqless & 0\\ & & & & \nu_1, & \nu_2, & \nu_3 & \ge & 0 \end{array} \]

Descrizione. \(\mu_t\) è il valore di un paio disponibile nel mese \(t\) e \(\nu_t\) il prezzo di un'ora di lavoro. L'obiettivo valuta la domanda a quei prezzi e sottrae \(r\, m_0 \sum_t \nu_t\), cioè le ore che i due operai iniziali offrono gratuitamente. Il primo gruppo di vincoli sono le colonne delle \(x_t\): produrre un paio vale \(\mu_t\) e consuma \(g\) ore al prezzo \(\nu_t\), e il saldo non supera il costo \(p_t\) delle materie prime. Il secondo sono le colonne delle \(s_t\): tenere un paio in magazzino fa guadagnare \(\mu_{t+1} - \mu_t\), che non può superare \(h_t\). Il terzo sono le colonne delle \(z_j\): un operaio assunto al mese \(j\) mette a disposizione \(r\) ore in ciascuno dei mesi da \(j\) in poi, e il loro valore non può superare quanto quell'assunzione costa.

Ricetta. \(\bar\nu_t = 0\): le ore di lavoro non si valutano, i vincoli sulle assunzioni sono soddisfatti perché il membro destro è positivo, e restano \(\mu_t \le p_t\) e \(\mu_{t+1} \le \mu_t + h_t\). Il valore più grande ammissibile si costruisce in avanti,

\[\bar\mu_1 = p_1, \qquad \bar\mu_t = \min(\bar\mu_{t-1} + h_{t-1},\ p_t).\]

Sull'istanza \(\bar\mu = (15, 15, 15)\) e \(\sum_t d_t\, \bar\mu_t = 15 \cdot 300 = 4500\), a cui va sommato il termine costante \(9000\):

\[\mathit{LB} = 4500 + 9000 = 13\,500 .\]

È il costo delle sole materie prime più il salario degli operai già in servizio: la manodopera aggiuntiva è regalata.

Soluzione ottima

\(t=1\) \(t=2\) \(t=3\)
produzione \(x_t\) 60 120 120
organico \(y_t\) 2 3 3
assunzioni \(z_t\) 0 1 0
scorta \(s_t\) 0 20 —
\(UB\) \(LB\) (duale) \(z(\mathit{LP})\) \(z(\mathit{LP}^+)\) \(z(\mathit{MILP})\) gap dell'euristica
18200 13500 15960 15960 16660 \(9{,}2\%\)

Piano ottimo

Il piano ottimo assume un operaio invece di due e anticipa venti paia dal terzo mese al secondo: pagare \(3\) euro di magazzino per venti paia costa \(60\) euro contro i \(1600\) di un'assunzione al terzo mese.

Considerazioni aggiuntive

  • Il vincolo di monotonia è ciò che lega i mesi dal lato dell'organico: se si potesse licenziare a costo zero, ogni mese sceglierebbe da solo il più piccolo organico che basta a produrre \(x_t\), e quella famiglia sparirebbe. Il problema non si spezzerebbe però in \(n\) problemi indipendenti: restano le scorte \(s_t\) a legare un mese al successivo.
  • Le variabili \(x_t\) e \(s_t\) sono dichiarate intere perché le paia di scarpe non si spezzano. Qui si potrebbero lasciare continue senza cambiare l'ottimo --- ma non perché la matrice sia totalmente unimodulare: non lo è, e si vede dal rilassamento, che vale \(15\,960\) contro un ottimo intero di \(16\,660\). La ragione è più locale: fissato l'organico \(y\) a valori interi, quello che resta in \((x, s)\) è un problema di bilancio con matrice d'intervallo, totalmente unimodulare, e termini noti interi (le domande, e le capacità \(r\, y_t / g = 40\, y_t\)); un suo vertice ottimo è quindi intero. È un argomento sul sottoproblema, non sul modello intero.
  • Il termine costante \(m_0\, w\, n\) va ricordato in ogni confronto: dimenticarlo fa apparire la formulazione \(A\) molto più economica della \(B\).

Domande di modellazione aggiuntive

9.2.1 — Straordinari

Ogni operaio può fare fino a \(40\) ore di straordinario al mese, pagate \(25\) euro l'ora. Come cambia il modello? Conviene usarli?

Una variante svolta: assunzioni molto costose

Il costo di assunzione sale da \(100\) a \(3000\) euro (selezione e formazione).

Cambia solo il dato \(u\). Il nuovo ottimo è \(19\,560\), cioè \(2900\) euro in più. Il piano non cambia: si assume comunque un operaio, perché senza quell'assunzione il problema è inammissibile (\(4 \cdot 120 / 160 = 3\) operai servono comunque). La differenza è esattamente \(3000 - 100 = 2900\): il costo aggiuntivo di un'assunzione inevitabile. È un esempio utile di analisi di sensitività fatta risolvendo di nuovo il modello.

La variante non aggiunge vincoli: cambia un dato, \(u = 3000\) invece di \(100\). Modello, duale e ricetta restano quelli del problema base, e il certificato non si accorge del cambiamento, perché \(\mu_t\) è il costo minimo per avere un paio disponibile al mese \(t\) e il costo di assunzione lì non compare: il bound resta \(13500\), mentre l'ottimo sale da \(16660\) a \(19560\). È il caso limite di questo capitolo: un bound duale corretto può essere del tutto insensibile al dato che muove il problema, e qui è l'euristica a seguirlo.

valore che cos'è
\(\mathit{UB}\) \(24000\) soluzione euristica
\(\mathit{LB}\) \(13500\) certificato duale costruito a mano
\(z(\mathit{LP})\) \(17490\) rilassamento senza i bound
\(z(\mathit{LP}^+)\) \(17490\) rilassamento con i bound
\(z(\mathit{MILP})\) \(19560\) ottimo del MILP

Codice

Script completo — python/fam09_2_manodopera.py (riproducibile con python3 python/fam09_2_manodopera.py dalla cartella python/). Notebook — notebooks/fam09_2_manodopera.ipynb — che si apre in Colab dal badge in cima alla pagina.

Mostra lo script completo — python/fam09_2_manodopera.py (273 righe)
"""Problema 9.2 -- Produzione e manodopera: due formulazioni equivalenti.

La stessa decisione scritta due volte: con le *assunzioni* z_t (formulazione A)
oppure con l'*organico* y_t (formulazione B). Si dimostra che hanno lo stesso
insieme di piani ammissibili e lo stesso ottimo, e si confrontano i rilassamenti.
E' il tema del capitolo 4: due formulazioni si confrontano solo dopo aver
dimostrato che descrivono lo stesso insieme intero.
"""
import gurobipy as gp
import pandas as pd
from gurobipy import GRB

from mip import (ammissibile, dualita_forte, due_rilassamenti, frazione, nuovo_modello,
                 rilassamenti,
                 registra_bound, rilassamenti, rilassamento, risolvi, valuta)
from stile import ARANCIO, BLU, ROSSO, TEAL, intestazione, plt, salva_dati, salva_figura
from esteso import salva_modello

R = range

# ---------- 1. MODELLO E ISTANZA ----------
intestazione("9.2 Produzione e manodopera: assunzioni (A) oppure organico (B)")
d2 = [60, 100, 140]        # domanda dei tre mesi (paia)
p2 = [15, 15, 15]          # costo di produzione per paio
h2 = [3, 3]                # costo di magazzino a fine mese
w2, r2, g2, u2, m2, r0 = 1500, 160, 4, 100, 2, 0
n2 = len(d2)
salva_dati(pd.DataFrame({"mese": R(1, n2 + 1), "domanda": d2, "costo_paio": p2}), "fam09_2_dati")
print(f"  {m2} operai all'inizio, {r2} h al mese ciascuno, {g2} h per paio: la capacita'")
print(f"  iniziale e' {m2 * r2 // g2} paia al mese. Salario {w2}, assunzione {u2}.")


def modello_A(d, p, h, w, r, g, u, m0, r0):
    """Formulazione A: z_t = quanti operai si assumono all'inizio del mese t."""
    n = len(d)
    mm = nuovo_modello("manodopera_A")
    x = mm.addVars(n, vtype=GRB.INTEGER, name="x")
    s = mm.addVars(n - 1, vtype=GRB.INTEGER, name="s")
    z = mm.addVars(n, vtype=GRB.INTEGER, name="z")
    mm.setObjective(gp.quicksum(p[t] * x[t] for t in R(n))
                    + gp.quicksum(h[t] * s[t] for t in R(n - 1))
                    + gp.quicksum((u + w * (n - t)) * z[t] for t in R(n)), GRB.MINIMIZE)
    mm.addConstr(x[0] - s[0] == d[0] - r0, name="bilancio[0]")
    mm.addConstrs((x[t] + s[t - 1] - s[t] == d[t] for t in R(1, n - 1)), name="bilancio")
    mm.addConstr(x[n - 1] + s[n - 2] == d[n - 1], name=f"bilancio[{n - 1}]")
    mm.addConstrs((-g * x[t] + gp.quicksum(r * z[j] for j in R(t + 1)) >= -r * m0
                   for t in R(n)), name="ore")
    return mm, x, s, z


def modello_B(d, p, h, w, r, g, u, m0, r0):
    """Formulazione B: y_t = quanti operai lavorano nel mese t (organico)."""
    n = len(d)
    mm = nuovo_modello("manodopera_B")
    x = mm.addVars(n, vtype=GRB.INTEGER, name="x")
    s = mm.addVars(n - 1, vtype=GRB.INTEGER, name="s")
    y = mm.addVars(n, vtype=GRB.INTEGER, name="y")
    # l'organico paga il salario ogni mese; le assunzioni sono gli incrementi y_t - y_{t-1}
    mm.setObjective(gp.quicksum(p[t] * x[t] for t in R(n))
                    + gp.quicksum(h[t] * s[t] for t in R(n - 1))
                    + gp.quicksum(w * y[t] for t in R(n))
                    + u * (y[n - 1] - m0), GRB.MINIMIZE)   # assunzioni totali = y_n - m0
    mm.addConstr(x[0] - s[0] == d[0] - r0, name="bilancio[0]")
    mm.addConstrs((x[t] + s[t - 1] - s[t] == d[t] for t in R(1, n - 1)), name="bilancio")
    mm.addConstr(x[n - 1] + s[n - 2] == d[n - 1], name=f"bilancio[{n - 1}]")
    mm.addConstrs((-g * x[t] + r * y[t] >= 0 for t in R(n)), name="ore")
    mm.addConstr(y[0] >= m0, name="organico_iniziale")
    mm.addConstrs((-y[t - 1] + y[t] >= 0 for t in R(1, n)), name="mai_licenziamenti")
    return mm, x, s, y


def duale_A(d, p, h, w, r, g, u, m0, r0):
    """max sum_t b_t mu_t - r m0 sum_t nu_t;  mu_t - g nu_t <= p_t;
    -mu_t + mu_{t+1} <= h_t;  r sum_{t >= j} nu_t <= u + w (n - j);  mu libere, nu >= 0."""
    n = len(d)
    b = [d[0] - r0] + d[1:n - 1] + [d[n - 1]]
    dl = nuovo_modello("duale_manodopera")
    mu = dl.addVars(n, lb=-GRB.INFINITY, name="mu")
    nu = dl.addVars(n, name="nu")
    dl.setObjective(gp.quicksum(b[t] * mu[t] for t in R(n))
                    - r * m0 * gp.quicksum(nu[t] for t in R(n)), GRB.MAXIMIZE)
    dl.addConstrs((mu[t] - g * nu[t] <= p[t] for t in R(n)), name="rc_x")
    dl.addConstrs((-mu[t] + mu[t + 1] <= h[t] for t in R(n - 1)), name="rc_s")
    dl.addConstrs((r * gp.quicksum(nu[t] for t in R(j, n)) <= u + w * (n - j) for j in R(n)),
                  name="rc_z")
    return dl


mA, xA, sA, zA = modello_A(d2, p2, h2, w2, r2, g2, u2, m2, r0)
salva_modello(mA, "fam09_2_primale")
mB, xB, sB, yB = modello_B(d2, p2, h2, w2, r2, g2, u2, m2, r0)
salva_modello(mB, "fam09_2_primale_b")
costante_A = m2 * w2 * n2          # il salario degli operai iniziali, fuori dal modello A
zA_val = risolvi(mA) + costante_A
zB_val = risolvi(mB)
print(f"  Formulazione A (assunzioni): z = {frazione(zA_val)} "
      f"(di cui {costante_A} di salario degli operai iniziali, termine costante)")
print(f"  Formulazione B (organico):   z = {frazione(zB_val)}")
assert abs(zA_val - zB_val) < 1e-6, (zA_val, zB_val)
print("  I due ottimi coincidono: le due formulazioni descrivono lo stesso problema.")
print("  Piano A: produzione " + ", ".join(frazione(xA[t].X) for t in R(n2))
      + "; assunzioni " + ", ".join(frazione(zA[t].X) for t in R(n2)))
print("  Piano B: produzione " + ", ".join(frazione(xB[t].X) for t in R(n2))
      + "; organico " + ", ".join(frazione(yB[t].X) for t in R(n2)))

# ---------- 2. IL RILASSAMENTO LP ----------
zlp2, zlp2r, _ = rilassamenti(mA)

# ---------- 3. IL DUALE DEL RILASSAMENTO (LOWER BOUND) ----------
dl2 = duale_A(d2, p2, h2, w2, r2, g2, u2, m2, r0)
salva_modello(dl2, "fam09_2_duale")
# ricetta: nu = 0 (le ore non si pagano) e mu_t = costo minimo per avere un paio al mese t
mu = []
for t in R(n2):
    mu.append(p2[t] if t == 0 else min(mu[t - 1] + h2[t - 1], p2[t]))
mano = {f"mu[{t}]": mu[t] for t in R(n2)}
lb2_var, viol = valuta(dl2, mano)
assert viol <= 1e-9, viol
lb2 = lb2_var + costante_A
print("  Duale a mano: nu = 0 (le ore di lavoro non si pagano) e mu_t = min(mu_{t-1}+h, p_t)")
print(f"    mu = " + ", ".join(frazione(v) for v in mu)
      + f"  ->  lb = {frazione(lb2_var)} + {costante_A} = {frazione(lb2)}")
dualita_forte(dl2, zlp2)
zlp2, zlp2r = zlp2 + costante_A, zlp2r + costante_A

# ---------- 4. EURISTICA COSTRUTTIVA (UPPER BOUND) ----------
intestazione("9.2 Euristica, duale e bound")
# euristica costruttiva: si produce la domanda del mese, e si assume solo quando le ore non bastano
organico, assunzioni, prod = m2, [0] * n2, []
for t in R(n2):
    prod.append(d2[t])
    servono = -(-g2 * d2[t] // r2)               # ceil
    if organico < servono:
        assunzioni[t] = servono - organico
        organico = servono
    print(f"  Mese {t + 1}: si producono {d2[t]} paia, servono "
          f"ceil({g2} * {d2[t]} / {r2}) = {servono} operai; organico {organico - assunzioni[t]}"
          f" -> se ne assumono {assunzioni[t]}")
ub2 = sum(p2[t] * prod[t] for t in R(n2)) \
    + sum(assunzioni[t] * (u2 + w2 * (n2 - t)) for t in R(n2)) + costante_A
sol_eur = {f"x[{t}]": prod[t] for t in R(n2)} | {f"z[{t}]": assunzioni[t] for t in R(n2)} \
    | {f"s[{t}]": 0 for t in R(n2 - 1)}
assert ammissibile(mA, sol_eur)
print(f"  Costo dell'euristica: ub = {frazione(ub2)}")
riga = registra_bound("2 manodopera", ub2, lb2, zlp2, zlp2r, zA_val)
salva_dati(pd.DataFrame([riga]), "fam09_2_bound")
assert lb2 <= zlp2 <= zA_val <= ub2 + 1e-9

# ---------- 5. L'EQUIVALENZA, VERIFICATA ----------
intestazione("9.2 L'equivalenza fra le due formulazioni, verificata")
print("  La corrispondenza e' y_t = m0 + sum_{j <= t} z_j, cioe' z_t = y_t - y_{t-1}")
print("  (con y_0 = m0). Sui piani ottimi:")
yA = [m2 + sum(round(zA[j].X) for j in R(t + 1)) for t in R(n2)]
print("    da A: organico implicito = " + ", ".join(str(v) for v in yA))
print("    da B: organico           = " + ", ".join(str(round(yB[t].X)) for t in R(n2)))
zB_implicite = [round(yB[0].X) - m2] + [round(yB[t].X) - round(yB[t - 1].X) for t in R(1, n2)]
print("    da B: assunzioni implicite = " + ", ".join(str(v) for v in zB_implicite))
assert sum(v * (u2 + w2 * (n2 - t)) for t, v in enumerate(zB_implicite)) + costante_A \
    == sum(round(zA[t].X) * (u2 + w2 * (n2 - t)) for t in R(n2)) + costante_A
print("  Il costo del personale coincide: A paga ogni assunzione una volta per tutti i mesi")
print("  che restano, B paga l'organico mese per mese. Stessa somma, contata in due modi.")

# ---------- 6. CONFRONTO DEI RILASSAMENTI DELLE DUE FORMULAZIONI ----------
zlpA, _, _ = rilassamento(mA, rafforzato=True)
zlpB, _, _ = rilassamento(mB, rafforzato=True)
print(f"  Rilassamenti: A -> {frazione(zlpA + costante_A)}   B -> {frazione(zlpB)}   "
      f"z(MILP) = {frazione(zA_val)}")
salva_dati(pd.DataFrame([{"formulazione": "A (assunzioni)", "z_lp": zlpA + costante_A,
                          "z_milp": zA_val},
                         {"formulazione": "B (organico)", "z_lp": zlpB, "z_milp": zB_val}]),
           "fam09_2_formulazioni")

# ---------- 7. DOMANDE DI MODELLAZIONE AGGIUNTIVE ----------
varianti = {}


def variante(nome, m, costante=0.0):
    z = risolvi(m) + costante
    print(f"  {nome:70s} z = {frazione(z)}")
    return z


# 2a: assumere costa molto di piu' (3000 invece di 100)
m, x, s, y = modello_B(d2, p2, h2, w2, r2, g2, 3000, m2, r0)
varianti["2a"] = variante("2a. L'assunzione costa 3000 euro invece di 100", m)
print("     organico: " + ", ".join(str(round(y[t].X)) for t in R(n2))
      + ";  produzione: " + ", ".join(str(round(x[t].X)) for t in R(n2)))
# 2b: straordinari, fino a 40 ore in piu' per operaio al mese, a 25 euro l'ora
m, x, s, y = modello_B(d2, p2, h2, w2, r2, g2, u2, m2, r0)
o = m.addVars(n2, name="o")
m.update()
for t in R(n2):
    m.chgCoeff(m.getConstrByName(f"ore[{t}]"), o[t], 1.0)   # le ore disponibili aumentano
m.addConstrs((o[t] <= 40 * y[t] for t in R(n2)), name="max_straordinari")
m.setObjective(m.getObjective() + gp.quicksum(25 * o[t] for t in R(n2)), GRB.MINIMIZE)
varianti["2b"] = variante("2b. Straordinari: fino a 40 h per operaio, 25 euro l'ora", m)
print("     straordinari usati: " + ", ".join(frazione(o[t].X) for t in R(n2))
      + "  (nessuno: anticipare la produzione e tenerla a magazzino costa meno)")
salva_dati(pd.DataFrame({"variante": list(varianti), "z": list(varianti.values())}),
           "fam09_2_varianti")

# ---------- 8. IL SANDWICH SULLA VARIANTE 2a ----------
intestazione("9.2a Il sandwich sulla variante: l'assunzione costa 3000 euro")
U2A = 3000

# La variante non aggiunge vincoli: cambia un dato. Modello e duale sono gli
# stessi, con u = 3000, e anche la ricetta duale e' la stessa --- cambia il
# valore che restituisce, ed e' proprio questo il punto: il certificato segue i
# dati senza che il modello si tocchi.
m2a, x2a, s2a, z2a = modello_A(d2, p2, h2, w2, r2, g2, U2A, m2, r0)
salva_modello(m2a, "fam09_2a_primale")
d2a_ = duale_A(d2, p2, h2, w2, r2, g2, U2A, m2, r0)
salva_modello(d2a_, "fam09_2a_duale")

# -- euristica ammissibile: la stessa regola, con il costo di assunzione nuovo --
print("Euristica costruttiva: si produce la domanda del mese e si assume solo quando le ore")
print("degli operai in servizio non bastano. La regola non cambia; cambia quanto costa.")
operai = m2
assunti = [0] * n2
for tt in R(n2):
    servono = -(-g2 * d2[tt] // r2)                      # arrotondamento per eccesso
    if servono > operai:
        assunti[tt] = servono - operai
        operai = servono
        print(f"  mese {tt + 1}: servono {servono} operai, se ne assumono {assunti[tt]}")
    else:
        print(f"  mese {tt + 1}: i {operai} operai in servizio bastano")
ub2a = (sum(p2[tt] * d2[tt] for tt in R(n2))
        + sum((U2A + w2 * (n2 - tt)) * assunti[tt] for tt in R(n2)) + costante_A)
sol_2a = ({f"x[{tt}]": d2[tt] for tt in R(n2)} | {f"s[{tt}]": 0 for tt in R(n2 - 1)}
          | {f"z[{tt}]": assunti[tt] for tt in R(n2)})
assert ammissibile(m2a, sol_2a), "la soluzione euristica della variante deve essere ammissibile"
print(f"  ub = {frazione(ub2a)}")

# -- certificato duale: la ricetta del problema base, sui dati nuovi --
mu_2a = []
for tt in R(n2):
    tetto = p2[tt]
    mu_2a.append(tetto if tt == 0 else min(mu_2a[tt - 1] + h2[tt - 1], tetto))
mano_2a = {f"mu[{tt}]": mu_2a[tt] for tt in R(n2)}
lb2a_var, viol_2a = valuta(d2a_, mano_2a)
assert viol_2a <= 1e-9, viol_2a
lb2a = lb2a_var + costante_A                 # come nel problema base: il salario degli
                                             # operai iniziali sta fuori dalla formulazione A
print("Soluzione duale a mano: nu = 0 (le ore si regalano) e mu_t = il costo piu' basso per")
print("  avere un paio disponibile al mese t, cioe' min(mu_{t-1} + h_{t-1}, p_t), come nel")
print(f"  problema base: mu = {[frazione(v) for v in mu_2a]}")
print(f"  ->  lb = {frazione(lb2a_var)} + {costante_A} = {frazione(lb2a)}")
zlp2a, zlp2ar, _ = due_rilassamenti(m2a, d2a_)
zlp2a, zlp2ar = zlp2a + costante_A, zlp2ar + costante_A
z2a_val = risolvi(m2a) + costante_A
riga_2a = registra_bound("2a assunzione a 3000", ub2a, lb2a, zlp2a, zlp2ar, z2a_val)
salva_dati(pd.DataFrame([riga_2a]), "fam09_2a_bound")
assert lb2a <= zlp2a <= z2a_val <= ub2a + 1e-9

# ---------- 9. FIGURA ----------
fig, ax = plt.subplots(figsize=(7.0, 3.2))
mesi = list(R(1, n2 + 1))
ax.bar(mesi, [xB[t].X for t in R(n2)], color=TEAL, width=0.55, label="produzione $x_t$")
ax.plot(mesi, d2, "o--", color=ROSSO, label="domanda $d_t$")
ax2 = ax.twinx()
ax2.step(mesi, [yB[t].X for t in R(n2)], where="mid", color=BLU, lw=2, label="organico $y_t$")
ax2.set_ylabel("operai", color=BLU)
ax2.set_ylim(0, max(yB[t].X for t in R(n2)) + 1.5)
ax2.grid(False)
ax.set_xticks(mesi)
ax.set_xlabel("mese")
ax.set_ylabel("paia")
ax.set_title(f"9.2: piano ottimo (z = {frazione(zB_val)})")
ax.legend(fontsize=8, loc="upper left")
ax2.legend(fontsize=8, loc="lower right")
salva_figura(fig, "cap09_manodopera_ottimo")
print("Fine.")