Analisi di sensitività e pricing

Risolvere un problema di programmazione lineare produce molto più della soluzione ottima. Il solver restituisce anche, gratuitamente, le variabili duali dei vincoli e i costi ridotti delle variabili, e queste due quantità rispondono alle domande che un decisore si pone subito dopo aver visto la soluzione: quale risorsa conviene aumentare? Quanto dovrebbe salire il prezzo di un prodotto non conveniente perché lo diventi? Quale nuovo prodotto vale la pena aggiungere al catalogo? Il laboratorio risponde prima nel modo ingenuo, risolvendo molti problemi modificati, e poi mostra che la teoria della dualità dà le stesse risposte con una sola risoluzione.

1. Il problema del product mix#

1.1 Formulazione#

Una sola classe di variabili, $x_p \ge 0$ litri del profumo $p$, e un vincolo di disponibilità per ingrediente:

$$\begin{array}{lrrrrrr} \max & 300 x_1 &+ 255 x_2 &+ 260 x_3 &+ 390 x_4\\ \text{s.t.} &0.01 x_1 &+ 0.05 x_2 &+ 0.07 x_3 &+ 0.04 x_4 &\le & 30\\ &0.34 x_1 &+ 0.45 x_2 &+ 0.36 x_3 &+ 0.51 x_4 &\le &400\\ &0.08 x_1 &+ 0.06 x_2 &+ 0.12 x_3 &+ 0.12 x_4 &\le & 90\\ &0.55 x_1 &+ 0.35 x_2 &+ 0.29 x_3 &+ 0.32 x_4 &\le &450\\ &0.02 x_1 &+ 0.09 x_2 &+ 0.16 x_3 &+ 0.01 x_4 &\le & 70\\ &x_1, &x_2, &x_3, &x_4 &\ge& 0 \end{array}$$

ovvero, in forma compatta, $\max\{c^\top x : Ax \le b,\ x \ge 0\}$ con $c$ il vettore dei prezzi, $A$ la matrice delle composizioni e $b$ il vettore delle disponibilità. (Il notebook scrive per un refuso $\min c^\top x$ nella forma compatta; il problema è di massimo, come dice la forma estesa e come fa il codice.) Le variabili sono continue: i litri si possono frazionare.

1.2 Una funzione che risolve il problema#

Poiché il problema verrà risolto molte volte con dati diversi, conviene incapsularlo in una funzione che riceve $A$, $b$, $c$ e restituisce il modello, la soluzione e il valore ottimo. Il notebook lascia da completare le quattro righe del modello:

import mip

c = [300, 255, 260, 390]
b = [30, 400, 90, 450, 70]

A = \
[[0.01, 0.05, 0.07, 0.04],
 [0.34, 0.45, 0.36, 0.51],
 [0.08, 0.06, 0.12, 0.12],
 [0.55, 0.35, 0.29, 0.32],
 [0.02, 0.09, 0.16, 0.01]]

def solve_productmix(A, b, c):
    """Risolve max{cx : Ax <= b, x >= 0} con A, b, c passati come argomenti"""
    n = len(c)     # numero di variabili
    k = len(b)     # numero di vincoli

    m = mip.Model()
    x = [m.add_var() for j in range(n)]

    for i in range(k):                          # un vincolo per ingrediente
        m.add_constr(mip.xsum(A[i][j] * x[j] for j in range(n)) <= b[i])

    m.objective = mip.maximize(mip.xsum(c[j] * x[j] for j in range(n)))
    m.optimize()

    return (m, [x[j].x for j in range(n)], m.objective_value)

_, solution, objective = solve_productmix(A, b, c)
print("Solution:", solution)
print("Objective: {0:10.2f}".format(objective))

L’underscore _ è il nome convenzionale per un valore restituito che non interessa (qui il modello). Output:

Solution: [573.115003808073, 89.74358974359002, 0.0, 323.0515359228229]
Objective:  320809.22

Si producono 573.1 litri del profumo 1, 89.7 del 2, 323.1 del 4 e niente del profumo 3, con ricavo 320809.22. Osserviamo che, con A[i][j], la riga i è l’ingrediente e la colonna j il profumo, e che gli indici Python partono da zero: il profumo 3 della tabella è solution[2].

2. Aumentare le risorse: la variazione di $b$#

Task 2. Se prezzi e composizioni non si toccano, l’unica leva è la disponibilità $b$. Aumentandola del 10% in tutte le componenti il ricavo cresce esattamente del 10%:

b2 = [b[i]*1.1 for i in range(len(b))]
_, solution2, objective2 = solve_productmix(A, b2, c)

print("Solution:", solution2)
print(f"Objective: {objective2}, increase: {((objective2 - objective) / objective) * 100: 6.2f}%")

Solution: [630.4265041888804, 98.71794871794896, 0.0, 355.3566895151053]
Objective: 352890.1370906322, increase:  10.00%

Il risultato è esatto, non approssimato: se tutti i termini di destra vengono moltiplicati per 1.1, la soluzione moltiplicata per 1.1 resta ammissibile e ottima, e con essa l’obiettivo. Ogni componente della soluzione è infatti cresciuta del 10% (573.1 diventa 630.4).

3. Quale vincolo è più restrittivo: le variabili duali#

Task 3. Per capire quale risorsa «stringe» di più conviene modificare le disponibilità una alla volta e osservare la variazione dell’obiettivo. Il notebook aumenta ogni $b_i$ di 0.1 e riottimizza:

for i in range(len(b)):
    b3 = b[:]      # [:] crea una copia: modificare b3 non deve toccare b
    b3[i] += .1
    _, sol, obj = solve_productmix(A, b3, c)
    print(f"Objective: {obj:10.2f}, Increase: {obj - objective:8.2f}, {(obj - objective)/objective * 100:5.2f}%")

Objective:  320809.22, Increase:     0.00,  0.00%
Objective:  320828.41, Increase:    19.19,  0.01%
Objective:  321020.96, Increase:   211.75,  0.07%
Objective:  320821.10, Increase:    11.88,  0.00%
Objective:  320809.22, Increase:     0.00,  0.00%

Due letture. Gli ingredienti 1 e 5 non producono alcuna variazione: i loro vincoli non sono attivi all’ottimo (la soluzione ne usa rispettivamente 23.14 litri su 30 e 22.77 su 70), quindi averne di più non serve. Gli ingredienti 2, 3 e 4 sono invece saturi, e fra questi il terzo è di gran lunga il più restrittivo: un decilitro in più vale 211.75 di ricavo, contro 19.19 e 11.88 degli altri due.

Perché aumentare e non diminuire? Perché il problema è di massimo e vogliamo sapere quanto guadagneremmo con più risorse. In un problema di minimo si farebbe il contrario: si ridurrebbe la risorsa per vedere di quanto scende il costo.

3.1 Non serve risolvere cinque problemi#

La stessa informazione è già contenuta nella soluzione del problema originale. A ogni vincolo di un LP corrisponde una variabile duale, il cui valore all’ottimo misura la variazione dell’obiettivo per un aumento unitario del termine di destra. In Python-MIP si legge dall’attributo .pi di ogni vincolo, e la lista dei vincoli del modello è m.constrs:

m, sol, obj = solve_productmix(A, b, c)
print("Dual values:", [con.pi for con in m.constrs])

Dual values: [-0.0, 191.92688499619217, 2117.479055597866, 118.81188118811873, -0.0]

Sono esattamente le variazioni misurate prima, moltiplicate per 10 (là l’incremento era 0.1, qui è unitario): $19.19 \cdot 10 = 191.9$, $211.75 \cdot 10 = 2117.5$, $11.88 \cdot 10 = 118.8$. Le variabili duali dei vincoli non attivi sono zero.

Il prezzo ombra vale localmente: aumentando molto una risorsa la soluzione cambia base e il tasso cambia. Per l’ingrediente 3 il tasso 2117.48 vale fra $b_3 \approx 73$ e $b_3 \approx 94.1$; oltre 94.1 litri l’ingrediente 3 smette di essere scarso, il suo vincolo diventa inattivo e il ricavo si ferma a 329528.25. Sotto i 73 litri il tasso è ancora più alto, perché la risorsa è ancora più preziosa.

4. Un prodotto non conveniente: i costi ridotti#

Task 4. Il profumo 3 non viene prodotto: $x_3 = 0$. A che prezzo comincerebbe a convenire? L’approccio ingenuo aumenta $c_3$ di 5 alla volta finché $x_3$ diventa positivo:

c_new = c[:]
sol = solution[:]
while sol[2] <= 1e-5:         # sol[2] è x_3 (indici da zero)
    c_new[2] += 5
    _, sol, obj = solve_productmix(A, b, c_new)
    print(f"With c_3 = {c_new[2]}: produce {sol[2]} liters of 3")

print(f"We need to set the price c_3 to {c_new[2]}")

With c_3 = 265: produce 0.0 liters of 3
With c_3 = 270: produce 0.0 liters of 3
...
With c_3 = 355: produce 0.0 liters of 3
With c_3 = 360: produce 116.60624892111176 liters of 3
We need to set the price c_3 to 360

Il confronto sol[2] <= 1e-5 al posto di == 0 protegge dagli arrotondamenti del solver. Il metodo funziona ma non è esatto: sappiamo solo che la soglia sta fra 355 e 360, e ha richiesto venti risoluzioni.

4.1 Il costo ridotto#

Anche questa informazione è già nel modello risolto. Il costo ridotto di una variabile si legge dall’attributo .rc; le variabili del modello sono nella lista m.vars, e la terza è m.vars[2]:

print(m.vars[2].rc)     # -97.64661081492756

Aggiungendo 97.6466 al prezzo originale si ottiene $260 + 97.6466 = 357.6466$, che sta proprio fra 355 e 360. La verifica fine del notebook prova due prezzi vicinissimi:

c_test = c[:]
c_test[2] = c[2] + 97.64
_, sol, obj = solve_productmix(A, b, c_test)
print(f"before threshold: x_3 = {sol[2]}")      # 0.0

c_test[2] = c[2] + 97.65
_, sol, obj = solve_productmix(A, b, c_test)
print(f"after threshold: x_3 = {sol[2]}")       # 116.60624892111176

Sotto la soglia non si produce, sopra si producono d’un colpo 116.6 litri: la soluzione salta da un vertice all’altro del poliedro.

Il calcolo per il profumo 3, con i duali $y^* = (0,\ 191.93,\ 2117.48,\ 118.81,\ 0)$ e la colonna $a_{\cdot 3} = (0.07, 0.36, 0.12, 0.29, 0.16)$: le risorse consumate valgono $0.36 \cdot 191.93 + 0.12 \cdot 2117.48 + 0.29 \cdot 118.81 = 69.09 + 254.10 + 34.46 = 357.65$, e $rc_3 = 260 - 357.65 = -97.65$. Il profumo 3 vale 260 ma «brucia» ingredienti che, usati altrove, rendono 357.65: finché il suo prezzo non copre quel valore non conviene. È lo stesso numero restituito da .rc, che il pacchetto ottiene dalla formula, non per tentativi.

c = [300, 255, 260, 390]
A = [[0.01, 0.05, 0.07, 0.04],
     [0.34, 0.45, 0.36, 0.51],
     [0.08, 0.06, 0.12, 0.12],
     [0.55, 0.35, 0.29, 0.32],
     [0.02, 0.09, 0.16, 0.01]]
y = [0.0, 191.92688499619217, 2117.479055597866, 118.81188118811873, 0.0]   # duali ottimi

for j in range(4):
    rc = c[j] - sum(A[i][j] * y[i] for i in range(5))
    print(f"profumo {j+1}: prezzo {c[j]}, valore risorse {c[j]-rc:8.2f}, costo ridotto {rc:8.2f}")

I profumi 1, 2 e 4, prodotti all’ottimo, hanno costo ridotto zero (a meno di arrotondamenti): il loro prezzo eguaglia esattamente il valore delle risorse che consumano.

5. Pricing: risolvere con un sottoinsieme di colonne#

La seconda parte del laboratorio allarga il catalogo: la fabbrica ha ricette per $n = 40$ profumi, $P_0 = \{0, 1, \ldots, 39\}$, con gli stessi 5 ingredienti e le stesse disponibilità. La matrice $A$ è ora $5 \times 40$ e $c$ ha 40 componenti; i primi quattro profumi coincidono con quelli della prima parte. (Da qui in avanti gli indici seguono la convenzione Python, da 0.)

Il vettore dei prezzi è

c = np.array([300, 255, 260, 390, 243, 93, 310, 117, 286, 267,
              250, 276, 324, 262, 216, 310, 285, 266, 280, 190,
              241, 301, 271, 252, 326, 279, 310, 278, 250, 248,
              319, 281, 287, 279, 299, 274, 305, 209, 250, 296])
b = np.array([30, 400, 90, 450, 70])

e la matrice $A$ delle composizioni è riportata nel notebook (5 righe di 40 valori). L’istruzione print(f"A has size {A.shape}, ...") verifica le dimensioni: (5, 40), (5,), (40,).

5.1 Una funzione con l’insieme delle colonne come argomento#

Task 1. Poiché risolveremo il problema con sottoinsiemi diversi di profumi, la funzione riceve anche la lista indices delle colonne da includere. Le variabili diventano un dizionario indicizzato dai profumi scelti, così x[p] ha sempre il numero «vero» del profumo:

import mip

I = range(5)   # ingredienti

def solve_productmix_indices(A, b, c, indices):

    m = mip.Model()

    P = [p for p in range(len(c)) if p in indices]   # solo i profumi richiesti

    x = {p: m.add_var() for p in P}

    for i in I:
        m.add_constr(mip.xsum(A[i][p]*x[p] for p in P) <= b[i])

    m.objective = mip.maximize(mip.xsum(c[p]*x[p] for p in P))

    m.optimize()

    return (m, [x[p].x for p in P], m.objective_value)

_, solution, objective = solve_productmix_indices(A, b, c, [0,1,2,3])

print("Solution:", solution)
print("Objective: {0:10.2f}".format(objective))

Con [0,1,2,3] si riottiene $x = (573.1,\ 89.7,\ 0,\ 323.1)$ e ricavo 320809.22: nulla è cambiato, sono gli stessi dati.

5.2 Quale colonna aggiungere: prima per tentativi#

Task 2. Il modo ingenuo: per ogni $p$ da 4 a 39 si risolve il problema con le colonne [0,1,2,3,p] e si guarda la variazione dell’obiettivo. Trentasei LP.

for p in range(4, 40):
    m, sol2, obj2 = solve_productmix_indices(A, b, c, [0,1,2,3,p])
    rndsol2 = [round(s, ndigits=1) for s in sol2]
    print(f"Added perfume {p:2d} --> {obj2:10.2f} (diff: {obj2 - objective:8.2f}). Sol: {rndsol2}")

Solo quattro profumi migliorano il ricavo; per tutti gli altri l’obiettivo resta 320809.22 e la nuova variabile vale zero.

Tre osservazioni del notebook. Ogni volta che entra un nuovo profumo, il profumo 1 esce dalla produzione (da 89.7 a 0): il nuovo prodotto usa meglio gli stessi ingredienti. Il profumo 3 aumenta sempre; il profumo 0 a volte aumenta e a volte diminuisce. Il ricavo non può mai diminuire aggiungendo una colonna: la vecchia soluzione, con la nuova variabile a zero, resta ammissibile.

5.3 Quale colonna aggiungere: con i costi ridotti#

Anche qui trentasei risoluzioni erano superflue. Il costo ridotto di una colonna si calcola dalla formula usando i duali della soluzione corrente, anche se la colonna non fa parte del modello: basta conoscere $c_p$ e la colonna $a_{\cdot p}$. Il vincolo duale associato a $x_p$ è $\sum_{i \in I} a_{ip}\, \eta_i \ge c_p$; se i duali correnti $\eta^*$ lo violano, cioè se

$$rc_p = c_p - \left(a_{0p} \eta^*_0 + a_{1p} \eta^*_1 + a_{2p} \eta^*_2 + a_{3p} \eta^*_3 + a_{4p} \eta^*_4\right) > 0,$$

allora aggiungere la colonna $p$ migliora l’obiettivo; se $rc_p \le 0$ la colonna non serve.

m, solution, objective = solve_productmix_indices(A, b, c, [0, 1, 2, 3])

dual = [c.pi for c in m.constrs]

for p in range(4, 40):
    rc = c[p] - sum(dual[i]*A[i][p] for i in I)
    print(f"Reduced cost for {p:2d} = {rc:9.2f}", end='')
    if rc > 0:
        print(" <======")
    else:
        print("")

I duali sono gli stessi della prima parte, $(0,\ 191.93,\ 2117.48,\ 118.81,\ 0)$, e i costi ridotti positivi sono esattamente quattro: profumo 14 (38.48), 17 (17.98), 23 (96.75) e 35 (40.27). Sono le stesse quattro colonne trovate per tentativi, con una sola risoluzione al posto di trentasei.

Il costo ridotto dice se una colonna migliora l’obiettivo, non di quanto: il profumo 35 ha costo ridotto 40.27, minore del 96.75 del profumo 23, ma la sua introduzione fa crescere il ricavo di più (7657 contro 6841). Il costo ridotto è il miglioramento per unità della nuova variabile; quanto la variabile potrà crescere dipende dalle risorse, e lo si scopre solo riottimizzando.

import numpy as np
A = np.array([[0.01, 0.05, 0.07, 0.04, 0.07, 0.17, 0.08, 0.14, 0.03, 0.01, 0.03, 0.03, 0.06, 0.06, 0.02, 0.05, 0.04, 0.01, 0.03, 0.04, 0.01, 0.02, 0.01, 0.02, 0.04, 0.02, 0.02, 0.11, 0.03, 0.02, 0.02, 0.03, 0.04, 0.06, 0.05, 0.04, 0.01, 0.05, 0.06, 0.04],
              [0.34, 0.45, 0.36, 0.51, 0.41, 0.15, 0.68, 0.39, 0.65, 0.27, 0.25, 0.11, 0.03, 0.69, 0.5 , 0.04, 0.33, 0.69, 0.32, 0.46, 0.53, 0.03, 0.39, 0.29, 0.68, 0.26, 0.53, 0.14, 0.3 , 0.76, 0.49, 0.49, 0.46, 0.15, 0.43, 0.27, 0.58, 0.53, 0.23, 0.63],
              [0.08, 0.06, 0.12, 0.12, 0.14, 0.04, 0.12, 0.05, 0.14, 0.14, 0.11, 0.15, 0.11, 0.19, 0.02, 0.17, 0.1 , 0.04, 0.1 , 0.04, 0.13, 0.13, 0.07, 0.01, 0.14, 0.14, 0.11, 0.2 , 0.12, 0.07, 0.14, 0.11, 0.09, 0.19, 0.15, 0.05, 0.24, 0.06, 0.16, 0.15],
              [0.55, 0.35, 0.29, 0.32, 0.36, 0.49, 0.12, 0.09, 0.12, 0.51, 0.56, 0.59, 0.79, 0.06, 0.33, 0.68, 0.48, 0.26, 0.55, 0.4 , 0.31, 0.75, 0.53, 0.66, 0.1 , 0.56, 0.31, 0.34, 0.55, 0.05, 0.33, 0.33, 0.41, 0.58, 0.31, 0.64, 0.01, 0.3 , 0.51, 0.08],
              [0.02, 0.09, 0.16, 0.01, 0.02, 0.15, 0.  , 0.33, 0.06, 0.07, 0.05, 0.12, 0.01, 0.  , 0.13, 0.06, 0.05, 0.  , 0.  , 0.06, 0.02, 0.07, 0.  , 0.02, 0.04, 0.02, 0.03, 0.21, 0.  , 0.1 , 0.02, 0.04, 0.  , 0.02, 0.06, 0.  , 0.16, 0.06, 0.04, 0.1 ]])
c = np.array([300, 255, 260, 390, 243, 93, 310, 117, 286, 267, 250, 276, 324, 262, 216, 310, 285, 266, 280, 190,
              241, 301, 271, 252, 326, 279, 310, 278, 250, 248, 319, 281, 287, 279, 299, 274, 305, 209, 250, 296])
dual = np.array([0.0, 191.92688499619217, 2117.479055597866, 118.81188118811873, 0.0])  # con le colonne 0..3

rc = c - dual @ A                      # tutti i 40 costi ridotti in una riga
for p in range(4, 40):
    if rc[p] > 0:
        print(f"profumo {p:2d}: costo ridotto {rc[p]:7.2f}  <== conviene aggiungerlo")

5.4 La sesta colonna#

Task 3. Aggiunto il profumo 14 (la scelta del notebook: non necessariamente quella con il costo ridotto massimo), si riottimizza e si ricalcolano i costi ridotti con i nuovi duali:

m, solution, objective = solve_productmix_indices(A, b, c, [0, 1, 2, 3, 14])
print(f"New objective: {objective}; solution = {solution}")

dual = [c.pi for c in m.constrs]
for p in range(4, 40):
    rc = c[p] - sum(dual[i]*A[i][p] for i in I)
    ...

Il nuovo ricavo è 322431.83 con $x = (589.1,\ 0,\ 0,\ 350.2,\ 42.2)$; i duali diventano $(0,\ 284.65,\ 1723.42,\ 118.81,\ 0)$. Il costo ridotto del profumo 14, ora in soluzione con valore positivo, è zero, come dev’essere per una variabile in base. Restano tre colonne con costo ridotto positivo: 12 (32.02), 23 (73.80) e 35 (34.93). Osserviamo che il profumo 12, prima non conveniente ($-8.54$), lo è diventato perché i prezzi ombra sono cambiati, e che il 17 non lo è più. Aggiungere il 14 non ha esaurito il potenziale di miglioramento.

5.5 Tutte le colonne#

Task 4. Con tutte e 40 le colonne il ricavo non può che essere maggiore o uguale: ogni soluzione ammissibile del problema ridotto lo è anche di quello completo, che quindi ha un ottimo almeno altrettanto buono. È l’argomento visto al task 2, applicato a molte colonne insieme.

m, solution, objective = solve_productmix_indices(A, b, c, range(40))

print(f"With all perfumes, objective: {objective}")
for p in range(40):
    if solution[p] > 0:
        print(f"{p:3d}: {solution[p]:10.2f}")

With all perfumes, objective: 337301.13572587486
  0:     121.50
  3:     573.16
 21:      71.43
 23:     221.50

Il ricavo sale a circa 337k producendo solo i profumi 0, 3, 21 e 23. Il profumo 14, aggiunto al passo precedente, è stato soppiantato da 21 e 23; il profumo 1 non si produce più. La sequenza di aggiunte non è quindi monotona: una colonna che entra può uscire quando ne arrivano altre, e le combinazioni finali restano combinatorie.

5.6 Verso la generazione di colonne#

La domanda naturale è se il procedimento «calcola i costi ridotti, aggiungi una colonna conveniente, riottimizza» possa essere reso automatico e portare all’ottimo del problema completo senza mai caricarlo tutto. La risposta è sì, ed è l’idea della generazione di colonne (column generation): si parte da poche colonne, a ogni iterazione si cerca la colonna con costo ridotto massimo (il problema di pricing) e la si aggiunge; quando nessuna colonna ha costo ridotto positivo, la soluzione corrente è ottima anche per il problema completo, perché tutte le variabili assenti soddisfano il criterio di ottimalità.

Sull’istanza dei 40 profumi, scegliendo a ogni passo la colonna di costo ridotto massimo:

Dopo tre aggiunte il ricavo è 337301.14, lo stesso del problema con tutte le 40 colonne, e il criterio di arresto lo certifica senza aver mai costruito le altre 33. Il metodo diventa decisivo quando le colonne sono troppe per essere elencate, come nei problemi di taglio o di scheduling dove ogni colonna è un pattern.

6. Riepilogo e glossario#

Elaborato dai due notebook del Lab 3 (sensitivity analysis, pricing) e dalla registrazione della terza lezione di laboratorio. Tutti i valori numerici sono stati ricalcolati.