Primi modelli con Python-MIP

Il pacchetto Python-MIP (modulo mip, documentazione su python-mip.com) permette di creare, manipolare e risolvere modelli con vincoli lineari e variabili continue, intere o binarie. Tutti i problemi che incontreremo nei laboratori hanno la stessa struttura: una funzione obiettivo, un insieme di vincoli e la natura delle variabili. Questo capitolo mostra come ciascuno di questi tre pezzi si traduce in una riga di codice, prima su esempi minimi e poi su un problema di produzione modellato da zero.

1. Il primo modello: due variabili, un vincolo#

Consideriamo il problema

$$\begin{array}{ll} \max & x_1 + x_2\\ \text{s.t.} & 2 x_1 + x_2 \le 10\\ & x_1, x_2 \ge 0 \end{array}$$

Prima di scrivere codice conviene capire la soluzione a occhio. Le due variabili hanno lo stesso peso nell’obiettivo, ma $x_1$ «costa» il doppio nel vincolo: ogni unità di $x_1$ consuma 2 unità di capacità, ogni unità di $x_2$ ne consuma una. Conviene quindi mettere tutto su $x_2$: $x_1 = 0$, $x_2 = 10$, valore dell’obiettivo 10. Il problema è banale di proposito, per concentrarsi sulla sintassi.

1.1 Installazione e importazione#

Ogni notebook comincia installando e importando il modulo:

!pip install mip
import mip

1.2 Il modello passo per passo#

Si crea un oggetto modello con il costruttore mip.Model(). Il nome (m nei notebook) è arbitrario; a questo oggetto si aggiungono poi variabili, vincoli e obiettivo.

m = mip.Model()

x1 = m.add_var()
x2 = m.add_var()

add_var() senza argomenti crea una variabile continua e non negativa: il limite inferiore predefinito è 0, quello superiore $+\infty$. Osserviamo che dopo queste righe non è stato risolto nulla: print(x1) mostra un oggetto variabile, print(m) un modello vuoto.

Il vincolo si aggiunge con add_constr, scrivendolo come espressione algebrica delle variabili; l’obiettivo si assegna all’attributo objective, avvolto in mip.maximize o mip.minimize:

m.add_constr(2*x1 + x2 <= 10)
m.objective = mip.maximize(x1 + x2)

I vincoli possono essere molti, l’obiettivo è uno solo: per questo i vincoli si aggiungono con un metodo mentre l’obiettivo si assegna a un attributo. Infine si chiama il solver e si leggono i risultati:

m.optimize()
print('solution:', x1.x, ',', x2.x)      # solution: 0.0 , 10.0
print('objective:', m.objective_value)   # 10.0

Il programma completo, nell’ordine che useremo sempre:

import mip

m = mip.Model()

x1 = m.add_var()
x2 = m.add_var()

m.objective = mip.maximize(x1 + x2)
m.add_constr(2*x1 + x2 <= 10)

m.optimize()

print('solution:', x1.x, ',', x2.x)

2. Un knapsack: variabili binarie e forma vettoriale#

Il secondo esempio è un problema di knapsack (zaino): ci sono oggetti con un valore e un peso, uno zaino con una capacità massima, e si vuole massimizzare il valore trasportato senza superare la capacità. La decisione per ogni oggetto è «lo prendo o no», quindi le variabili sono binarie.

$$\begin{array}{lll} \max & 3 x_1 + 4 x_2 + 7 x_3 + 5 x_4\\ \text{s.t.} & 4 x_1 + 5 x_2 + 6 x_3 + 4 x_4 \le 13\\ & x_1, x_2, x_3, x_4 \in \{0,1\} \end{array}$$

2.1 La natura delle variabili#

La differenza rispetto al primo esempio è tutta nel parametro var_type di add_var:

In tutti i casi la variabile è non negativa per default. Si può anche dare un nome alla variabile con name='x1', utile quando si stampa il modello.

2.2 Dati, variabili in lista e xsum#

Con quattro variabili scrivere x1 = m.add_var(...), x2 = ... e così via è già noioso; con centinaia diventa impraticabile. La soluzione è vettorizzare: si definiscono i dati come liste e si creano le variabili con una list comprehension, nell’ordine dati, modello, variabili, vincoli, obiettivo.

import mip

# dati: due vettori (valore e peso) e la capacità
value = [3, 4, 7, 5]
weight = [4, 5, 6, 4]
max_weight = 13

m = mip.Model()

x = [m.add_var(var_type=mip.BINARY) for i in range(4)]

Ora x è una lista di quattro variabili, x[0], …, x[3]: gli indici partono da zero, come sempre in Python, quindi la variabile matematica $x_1$ è x[0]. Per scrivere la somma $\sum_i w_i x_i$ si usa mip.xsum, con dentro un ciclo for in una riga:

m.add_constr(mip.xsum(weight[i] * x[i] for i in range(4)) <= max_weight)

m.objective = mip.maximize(mip.xsum(value[i] * x[i] for i in range(4)))

L’espressione weight[i] * x[i] for i in range(4) genera i quattro prodotti $w_i x_i$; xsum li somma in un’unica espressione lineare, che viene poi confrontata con max_weight. Per esempio con weight = [4, 5, 6, 4] la somma è $4x_0 + 5x_1 + 6x_2 + 4x_3$. Il modulo accetta anche l’operatore += al posto di add_constr: m += mip.xsum(...) <= max_weight aggiunge il vincolo al modello, e nei notebook le due forme si alternano.

m.optimize()
print([x[i].x for i in range(4)])      # [1.0, 1.0, 0.0, 1.0]
print(m.objective_value)               # 12.0

La soluzione prende gli oggetti 1, 2 e 4 (peso $4 + 5 + 4 = 13$, valore $3 + 4 + 5 = 12$). L’oggetto 3, il più prezioso, resta fuori: con i suoi 6 di peso lascerebbe spazio per un solo altro oggetto e il valore massimo sarebbe $7 + 5 = 12$ con peso 10, oppure $7 + 4 = 11$. A parità di valore il solver riporta una delle soluzioni ottime. Si può stampare la soluzione in due modi equivalenti: ciclando sugli indici, [x[i].x for i in range(4)], oppure ciclando direttamente sulle variabili, [v.x for v in x].

import mip

m = mip.Model(sense=mip.MAXIMIZE)

x1 = m.add_var(var_type=mip.BINARY)
x2 = m.add_var(var_type=mip.BINARY)
x3 = m.add_var(var_type=mip.BINARY)
x4 = m.add_var(var_type=mip.BINARY)

m.add_constr(4*x1 + 5*x2 + 6*x3 + 4*x4 <= 13)
m.objective = mip.maximize(3*x1 + 4*x2 + 7*x3 + 5*x4)

m.optimize()
print("Solution:", x1.x, "", x2.x, "", x3.x, "", x4.x)
# print("Solution: {0}, {1}, {2:.8f}, {3}".format(x1.x, x2.x, x3.x, x4.x))

3. La fabbrica di cellulari: dal testo al modello#

Da qui in avanti i laboratori seguono sempre lo stesso schema: una descrizione del problema, la sua modellazione matematica, la traduzione in Python-MIP e la soluzione.

Un modo per leggere la tabella: guardando solo i display si possono produrre al più 10 M1 e nessun M2, oppure 5 M2 e nessun M1. Ogni riga impone un vincolo analogo, e la soluzione deve rispettarli tutti insieme.

3.1 I quattro passi della modellazione#

Conviene seguire sempre la stessa scaletta: insiemi, parametri, variabili decisionali, modello.

Insiemi. Due: l’insieme dei modelli di cellulare $I = \{1, 2\}$, con indice $i$, e l’insieme dei componenti $J = \{1, \ldots, 6\}$, con indice $j$.

Parametri. Sono tutte le informazioni numeriche del testo, e solo quelle:

• la disponibilità $b_j$ del componente $j$, il vettore $(10, 18, 12, 21, 9, 10)$;
• il fabbisogno $a_{ji}$ del componente $j$ per produrre un’unità del modello $i$, la matrice $6 \times 2$ della tabella;
• il prezzo di vendita $c_i$ del modello $i$, il vettore $(110, 130)$.

Variabili decisionali. $x_i$ = numero di cellulari del modello $i$ da produrre. Non si producono 2,3 cellulari: $x_i$ è intera e non negativa. Per convenzione le lettere iniziali dell’alfabeto ($a, b, c$) si riservano ai parametri e quelle finali ($x, y, z$) alle variabili.

Modello. Si massimizza il ricavo, somma dei prezzi per le quantità, rispettando la disponibilità di ogni componente:

$$\begin{array}{lll} \max & \displaystyle\sum_{i \in I} c_i x_i\\ \text{s.t.} & \displaystyle\sum_{i \in I} a_{ji} x_i \le b_j & \forall j \in J\\ & x_i \in \mathbb{Z}_+ & \forall i \in I \end{array}$$

Il «per ogni $j \in J$» genera un vincolo per componente. Scritto per esteso, con i numeri della tabella, il modello è

$$\begin{array}{llll} \max & 110 x_1 + 130 x_2\\ \text{s.t.} & x_1 + 2 x_2 & \le 10&\quad\text{(display)}\\ & 2 x_1 + 2 x_2 & \le 18&\quad\text{(memoria)}\\ & x_1 + 3 x_2 & \le 12&\quad\text{(fotocamera)}\\ & 2 x_1 + 3 x_2 & \le 21&\quad\text{(CPU)}\\ & x_1 & \le 9 &\quad\text{(fotocamera termica)}\\ & x_1 & \le 10&\quad\text{(rx/tx satellitare)}\\ & x_1, x_2 \in \mathbb Z_+. \end{array}$$

Le due scritture sono equivalenti: la forma con le sommatorie è quella parametrica, che si legge allo stesso modo con due prodotti o con duecento; la forma esplicita è comoda solo quando il modello è piccolo.

3.2 Traduzione esplicita#

La traduzione riga per riga del modello esplicito. Le variabili sono intere, quindi var_type=mip.INTEGER; l’esercizio del notebook lascia da completare proprio questo argomento e i vincoli.

import mip

m = mip.Model()

x1 = m.add_var(var_type=mip.INTEGER)
x2 = m.add_var(var_type=mip.INTEGER)

m.objective = mip.maximize(110*x1 + 130*x2)

m.add_constr(  x1 + 2*x2 <= 10)
m.add_constr(2*x1 + 2*x2 <= 18)
m.add_constr(  x1 + 3*x2 <= 12)
m.add_constr(2*x1 + 3*x2 <= 21)
m.add_constr(  x1        <=  9)
m.add_constr(  x1        <= 10)

m.optimize()

print(x1.x, x2.x)        # 8.0 1.0

La soluzione ottima produce 8 unità di M1 e 1 di M2, con ricavo $110 \cdot 8 + 130 \cdot 1 = 1010$. Vale la pena verificarla sui vincoli: i display usati sono $8 + 2 = 10$ su 10 e la memoria $16 + 2 = 18$ su 18, entrambi saturi; fotocamere $8 + 3 = 11$ su 12, CPU $16 + 3 = 19$ su 21, termiche 8 su 9. Sono display e memoria a limitare la produzione.

Osserviamo che l’ottimo intero coincide qui con il vertice del poligono: il rilassamento continuo del problema ha già soluzione intera $(8, 1)$. Non è la regola, ed è il motivo per cui la natura delle variabili va sempre dichiarata.

4. La forma parametrica#

Lo stesso problema in forma matriciale: definiti

$$x = \begin{pmatrix}x_1\\x_2\end{pmatrix},\qquad c = \begin{pmatrix}110\\130\end{pmatrix},\qquad A = \begin{pmatrix}1&2\\2&2\\1&3\\2&3\\1&0\\1&0\end{pmatrix},\qquad b = \begin{pmatrix}10\\18\\12\\21\\9\\10\end{pmatrix},$$

il modello si scrive

$$\begin{array}{lll} \max & c^\top x\\ \text{s.t.} & A x \le b\\ & x \in \mathbb Z_+^2 \end{array}$$

che è la forma standard con cui si presenta ogni problema di programmazione lineare. In codice, i dati stanno in liste separate dal modello, così potrebbero anche essere letti da un file. Una matrice è semplicemente una lista di liste: A[j] è la riga $j$, A[j][i] l’elemento di riga $j$ e colonna $i$.

import mip

n_model = 2
n_component = 6

c = [110, 130]

A = [[1,2],
     [2,2],
     [1,3],
     [2,3],
     [1,0],
     [1,0]]

b = [10,18,12,21,9,10]

Le variabili si creano in una riga, anche se sono solo due; l’obiettivo è una xsum sui modelli; i vincoli si generano con un ciclo for sui componenti, e dentro ognuno la somma corre sui modelli:

m = mip.Model()
x = [m.add_var(var_type=mip.INTEGER) for i in range(n_model)]

m.objective = mip.maximize(mip.xsum(c[i]*x[i] for i in range(n_model)))

for j in range(n_component):
    m += mip.xsum(A[j][i]*x[i] for i in range(n_model)) <= b[j]

m.optimize()
print([x[i].x for i in range(n_model)])     # [8.0, 1.0]

Il ciclo for j in range(n_component) è la traduzione letterale del «$\forall j \in J$» del modello, e mip.xsum(A[j][i]*x[i] for i in range(n_model)) è la traduzione di $\sum_{i \in I} a_{ji} x_i$. Un modo per non sbagliare gli indici è partire dal termine di destra: si scrive prima <= b[j], poi la somma che deve rispettarlo, e ci si chiede su quale indice corre la somma (qui $i$, i modelli) e quale invece è fissato dal vincolo (qui $j$, il componente). La soluzione è la stessa, $(8, 1)$, perché è lo stesso modello scritto in un altro modo.

import mip

n_model = 2
n_component = 6

c = [110, 130]

A = [[1,2],
     [2,2],
     [1,3],
     [2,3],
     [1,0],
     [1,0]]

b = [10,18,12,21,9,10]

m = mip.Model()
x = [m.add_var(var_type=mip.INTEGER) for i in range(n_model)]

m.objective = mip.maximize(mip.xsum(c[i]*x[i] for i in range(n_model)))

for j in range(n_component):
    m.add_constr(mip.xsum(A[j][i]*x[i] for i in range(n_model)) <= b[j])

m.optimize()

print([x[i].x for i in range(n_model)])

5. Errori tipici della prima lezione#

Durante il laboratorio sono emersi alcuni inciampi ricorrenti, tutti facili da evitare una volta nominati.

• Natura delle variabili dimenticata. Senza var_type la variabile è continua: il solver può restituire 8,5 cellulari. Nel knapsack, senza mip.BINARY, la soluzione riempirebbe lo zaino con frazioni di oggetti.
• Matrice trascritta male. Nella versione parametrica l’ultima riga di A è stata discussa a lungo in aula, perché la formula del notebook e la tabella non concordano (si veda la nota sul refuso più sopra). Conviene sempre rileggere A riga per riga contro la tabella del testo, che è la fonte dei dati, e non contro una formula intermedia.
• Celle non rilanciate. Dopo aver corretto un dato bisogna rieseguire la cella dei dati e poi tutte le celle successive, altrimenti il modello continua a usare i valori vecchi (si veda il troubleshooting del capitolo 0).
• Completamento automatico. Colab suggerisce codice mentre si scrive; è comodo, ma il suggerimento non è sempre quello giusto. Ogni riga proposta va letta prima di accettarla.
• Uguaglianza nei vincoli. Un vincolo di uguaglianza si scrive con ==, non con =.

6. Riepilogo#

Il messaggio da portare via è duplice. Da un lato la scaletta di modellazione: insiemi, parametri, variabili decisionali, modello (obiettivo, vincoli, natura delle variabili). Dall’altro la traduzione meccanica in Python-MIP:

Nel prossimo capitolo la stessa scaletta si applica a problemi su grafi, dove insiemi e variabili sono indicizzati da coppie di nodi e le liste lasciano il posto ai dizionari.

Elaborato dai tre notebook del Lab 1 (firstMIP, cellphone, cellphone-parametric) e dalla registrazione della prima lezione di laboratorio.