Ottimizzare con Pyomo: guida completa passo dopo passo

Configurare Pyomo

Prima di immergerci nella modellazione, dobbiamo predisporre l’ambiente installando Pyomo e scegliendo un solver appropriato.

1. Prerequisiti

Per usare pyomo devi avere Python 3.6 o superiore. Pyomo può essere installato tramite pip.

pip install pyomoOfferto da

Questo tutorial è stato creato utilizzando la versione di pyomo 6.8.0.

import pyomo
print(pyomo.__version__)Offerto da

Output:

>>> 6.8.0Offerto da

2. Scelta e installazione del solver giusto

Nell’ottimizzazione, i solver sono essenziali poiché sono gli algoritmi che trovano la soluzione ottimale al problema che hai definito. Solver diversi sono più adatti a seconda del tipo di problema (ad es. lineare, non lineare, intero). Pyomo è uno strumento di modellazione che si affida a solver esterni per eseguire il calcolo effettivo.

Rivediamo alcuni dei solver più comuni.

Solver open-source

1. GLPK (GNU Linear Programming Kit)

GLPK è uno strumento popolare per risolvere problemi di programmazione lineare (LP) e programmazione mista intera (MIP). 

È un’ottima scelta per compiti di ottimizzazione lineare di base ed è ampiamente usato in ambito accademico e industriale.

Installazione

• Windows: Scarica e installa GLPK.
• macOS: 

brew install glpkOfferto da

• Linux: 

sudo apt-get install glpk-utilsOfferto da

2. CBC (Coin-or Branch and Cut)

CBC è un solver open-source per problemi LP e MIP. 

Offre funzionalità avanzate e, in alcuni casi, prestazioni migliori rispetto a GLPK, rendendolo una valida opzione per compiti di ottimizzazione più complessi. 

CBC può essere installato tramite il gestore pacchetti conda.

conda install -c conda-forge coincbcOfferto da

3. IPOPT (Interior Point OPTimizer)

IPOPT è un solver potente progettato per la programmazione non lineare (NLP) su larga scala. 

È particolarmente adatto alla gestione di modelli non lineari complessi, rendendolo un’ottima scelta per problemi oltre l’ottimizzazione lineare. 

IPOPT può essere installato tramite il gestore pacchetti conda.

!conda install -c conda-forge ipoptOfferto da

Solver commerciali

1. CPLEX

CPLEX è un solver all’avanguardia che gestisce in modo efficiente problemi LP, MIP e QP (programmazione quadratica). 

Richiede una licenza IBM ma è disponibile gratuitamente per utenti accademici, risultando un’ottima scelta per ricerca e didattica.

2. Gurobi

Gurobi è un solver commerciale leader, noto per velocità ed efficienza nel risolvere problemi LP, MIP, QP e NLP (non lineari). 

Come CPLEX, richiede una licenza ma offre accesso gratuito agli utenti accademici. È quindi uno strumento standard di settore per l’ottimizzazione avanzata.

Open-source vs. solver commerciali

I solver open-source come GLPK e CBC sono gratuiti e sufficienti per molte esigenze di ottimizzazione di base. Sono ottime scelte per progetti di piccola scala e scopi educativi. 

Tuttavia, i solver commerciali come CPLEX e Gurobi offrono in genere prestazioni superiori, soprattutto per problemi più grandi e complessi. Questi solver hanno funzionalità avanzate, tra cui un miglior supporto per la programmazione quadratica e non lineare, e sono ottimizzati per applicazioni industriali su larga scala.

Sebbene i solver open-source possano gestire molti compiti di ottimizzazione di routine, i solver commerciali sono spesso una scelta migliore quando si affrontano esigenze più complesse e ad alte prestazioni. 

Ricorda che i solver commerciali richiedono una licenza, sebbene siano disponibili gratuitamente per utenti accademici.

Ora vediamo come configurare un solver in Pyomo. In questo caso userò GLPK. 

3. Configurare un solver in Pyomo

Per prima cosa, dopo l’installazione assicurati che l’eseguibile del solver sia nel PATH di sistema. 

Poi crea uno script Python e aggiungi quanto segue:

from pyomo.environ import SolverFactory
solver = SolverFactory('glpk')Offerto da

Per confermare che sia Pyomo sia il solver siano installati correttamente, risolviamo un semplice problema di test.

Problema di test: semplice programma lineare

Obiettivo: Minimizzare Z=x+y

Soggetto a:

• x + 2y ≥ 4
• x - y ≤ 1
• x ≥ 0
• y ≥ 0

Questo problema riguarda la ricerca del valore più piccolo di Z, che è la somma di due variabili, x e y. Tuttavia, x e y devono soddisfare alcune condizioni. 

Primo, quando sommi x e due volte y, il risultato deve essere almeno 4. Secondo, x meno y deve essere minore o uguale a 1. Infine, sia x sia y devono essere numeri maggiori o uguali a zero (non possono essere negativi). 

L’obiettivo è trovare valori di x e y che soddisfino queste condizioni rendendo Z il più piccolo possibile.

Implementazione con Pyomo:

import pyomo.environ as pyo
# Create a model
model = pyo.ConcreteModel()
# Define variables
model.x = pyo.Var(within=pyo.NonNegativeReals)
model.y = pyo.Var(within=pyo.NonNegativeReals)
# Define objective
model.obj = pyo.Objective(expr=model.x + model.y, sense=pyo.minimize)
# Define constraints
model.con1 = pyo.Constraint(expr=model.x + 2 * model.y >= 4)
model.con2 = pyo.Constraint(expr=model.x - model.y <= 1)
# Select solver
solver = pyo.SolverFactory('glpk')
# Solve the problem
result = solver.solve(model)
# Display results
print('Status:', result.solver.status)
print('Termination Condition:', result.solver.termination_condition)
print('Optimal x:', pyo.value(model.x))
print('Optimal y:', pyo.value(model.y))
print('Optimal Objective:', pyo.value(model.obj))Offerto da

Se tutto funziona correttamente, l’output atteso sarà:

Status: ok
Termination Condition: optimal
Optimal x: 0.0
Optimal y: 2.0
Optimal Objective: 2.0Offerto da

Esaminiamo il codice sopra: per prima cosa definisce due variabili, x e y, che possono assumere solo valori non negativi. L’obiettivo del modello è minimizzare la somma di x e y (x + y). Il codice definisce il solver come glpk per trovare i valori ottimali di x e y che soddisfano questi vincoli minimizzando l’obiettivo.

Dopo l’esecuzione, troviamo che i valori ottimali per le variabili sono x = 0.0 e y = 2.0, che minimizzano la funzione obiettivo Z = x + y. Pertanto, il valore minimo della funzione obiettivo è 2.0, che soddisfa i vincoli dati.

Basi della modellazione con Pyomo

Capire come definire i componenti di base di un modello di ottimizzazione in Pyomo è necessario per impostare e risolvere efficacemente i problemi.

1. Definire le variabili

Le variabili rappresentano le decisioni da prendere in un problema di ottimizzazione. In Pyomo, le variabili sono le quantità che il solver regolerà per ottimizzare la funzione obiettivo rispettando tutti i vincoli.

Variabili scalari

Una variabile scalare è una singola variabile non indicizzata su alcun insieme. Per definire una variabile scalare in Pyomo, si usa la classe Var dal modulo pyomo.environ.

from pyomo.environ import Var
model.x = Var()Offerto da

Per prima cosa importiamo Var e creiamo una variabile x usando Var(). Questa variabile non ha limiti specificati, quindi può assumere qualunque valore reale a meno che il modello non la vincoli diversamente.

Aggiungere limiti

Puoi restringere i valori che una variabile può assumere specificando dei limiti. I limiti sono definiti come una tupla (lower_bound, upper_bound):

from pyomo.environ import Var
model.x = Var(bounds=(0, None))Offerto da

Specificare i domini

Pyomo fornisce domini predefiniti per specificare il tipo di valori che una variabile può assumere, come NonNegativeReals, Integers o Binary: 

from pyomo.environ import Var, NonNegativeReals
model.x = Var(domain=NonNegativeReals)Offerto da

Variabili indicizzate

Quando si lavora con più variabili simili per natura, ad esempio variabili che rappresentano periodi di tempo o articoli diversi, è efficiente usare variabili indicizzate. Le variabili indicizzate sono definite su un insieme.

import pyomo.environ as pyo
model.I = pyo.Set(initialize=[1, 2, 3])
model.y = pyo.Var(model.I, domain=pyo.NonNegativeReals)Offerto da

Supponiamo che tu stia modellando le quantità prodotte per tre prodotti. Puoi definire:

model.Products = pyo.Set(initialize=['A', 'B', 'C'])
model.production = pyo.Var(model.Products, domain=pyo.NonNegativeReals)Offerto da

Ora, model.production['A'], model.production['B'] e model.production['C'] rappresentano le quantità prodotte per i prodotti A, B e C, rispettivamente.

2. Definire gli obiettivi

La funzione obiettivo è ciò che vogliamo ottimizzare (minimizzare o massimizzare). Definisce lo scopo del modello, come minimizzare i costi o massimizzare i profitti, ed è tipicamente espressa come un’equazione matematica che coinvolge le variabili decisionali. 

Queste sono definite tramite la classe Objective:

from pyomo.environ import ConcreteModel, Var, Objective, minimize, maximize, NonNegativeReals
# Create a model
model = ConcreteModel()
# Define variables
model.x = Var(within=NonNegativeReals)
model.y = Var(within=NonNegativeReals)
# Minimization (cost)
model.cost = Objective(expr=2 * model.x + 3 * model.y, sense=minimize)
# When Maximization profit - (can have one objective at a time)
# model.profit = Objective(expr=5 * model.x + 4 * model.y, sense=maximize)Offerto da

3. Aggiungere vincoli

I vincoli definiscono i limiti o i requisiti del problema:

from pyomo.environ import Constraint
model.con1 = Constraint(expr=model.x + model.y >= 10)Offerto da

L’esempio sopra definisce un vincolo nel modello Pyomo usando la classe Constraint. Il vincolo model.con1 specifica che la somma delle variabili x e y deve essere maggiore o uguale a 10.

4. Parametrizzare i modelli

I parametri sono valori fissi usati nel modello per rappresentare quantità note o costanti che non cambiano durante il processo di ottimizzazione. 

Aiutano a definire le relazioni tra variabili e vincoli, dando struttura al modello integrando dati reali o assunzioni:

from pyomo.environ import Param
model.p = Param(initialize=5)Offerto da

Il codice sopra definisce un parametro p nel modello Pyomo usando la classe Param e lo inizializza con il valore fisso 5. Il parametro p può ora essere usato nel modello per rappresentare una costante che non cambia durante l’ottimizzazione.

Ora lavoriamo su un problema di ottimizzazione end-to-end!

Esempio end-to-end con Pyomo

Vediamo un esempio end-to-end di risoluzione di un problema di ottimizzazione con Pyomo. Modelliamo uno scenario reale in cui una fabbrica produce due prodotti e l’obiettivo è massimizzare il profitto considerando i vincoli di tempo macchina. 

1. Descrizione del problema

Una fabbrica produce due prodotti, P1 e P2. Il profitto per unità è:

• P1: $40
• P2: $50

Tempo macchina disponibile:

• Macchina A: 100 ore
• Macchina B: 80 ore
• Macchina C: 90 ore

Tempo richiesto per unità:

Obiettivo: Massimizzare il profitto.

Variabili decisionali:

• x₁: unità di P1 da produrre.
• x₂: unità di P2 da produrre.

2. Formulazione matematica

Funzione obiettivo:

Massimizzare Z = 40x₁ + 50x₂

Vincoli:

1. Capacità A: 1x₁ + 2x₂ ≤ 100
2. Capacità B: 2x₁ + 1x₂ ≤ 80
3. Capacità C: 3x₂ ≤ 90
4. Non negatività: x₁, x₂ ≥ 0

3. Implementazione

In base all’obiettivo e ai vincoli del problema, ecco il codice Python per modellarlo, sempre usando GLPK.

# Step 1: Import Libraries
import pyomo.environ as pyo
# Step 2: Create a Concrete Model
model = pyo.ConcreteModel()
# Step 3: Define Decision Variables (Units of P1 and P2 to produce)
model.x1 = pyo.Var(within=pyo.NonNegativeReals)
model.x2 = pyo.Var(within=pyo.NonNegativeReals)
# Step 4: Define the Objective Function (Maximize profit)
model.profit = pyo.Objective(expr=40 * model.x1 + 50 * model.x2, sense=pyo.maximize)
# Step 5: Define Constraints
# Machine A capacity constraint: 1x1 + 2x2 <= 100
model.machine_a = pyo.Constraint(expr=1 * model.x1 + 2 * model.x2 <= 100)
# Machine B capacity constraint: 2x1 + 1x2 <= 80
model.machine_b = pyo.Constraint(expr=2 * model.x1 + 1 * model.x2 <= 80)
# Machine C capacity constraint: 3x2 <= 90
model.machine_c = pyo.Constraint(expr=3 * model.x2 <= 90)
# Step 6: Solve the Model using GLPK solver
solver = pyo.SolverFactory('glpk')
result = solver.solve(model)
# Step 7: Analyze Results
# Display Solver Status and Termination Condition
print('Solver Status:', result.solver.status)
print('Termination Condition:', result.solver.termination_condition)
# Get and display the optimal values for x1, x2, and the maximum profit
x1_opt = pyo.value(model.x1)
x2_opt = pyo.value(model.x2)
profit_opt = pyo.value(model.profit)
print(f'Optimal production of P1 (x1): {x1_opt}')
print(f'Optimal production of P2 (x2): {x2_opt}')
print(f'Maximum Profit: ${profit_opt}')Offerto da

Output:

>>> Solver Status: ok
>>> Termination Condition: optimal
>>> Optimal production of P1 (x1): 25.0
>>> Optimal production of P2 (x2): 30.0
>>> Maximum Profit: $2500.0Offerto da

Nel codice sopra, definiamo un modello di ottimizzazione lineare per massimizzare il profitto dalla produzione di due prodotti (P1 e P2). La funzione obiettivo massimizza il profitto, con ogni unità di P1 che contribuisce $40 e ogni unità di P2 che contribuisce $50.

Imponiamo tre vincoli che rappresentano i limiti di tempo macchina per le macchine A, B e C. 

Infine, usiamo il solver GLPK per risolvere il problema. 

La risposta finale è produrre 25 unità di P1 e 30 unità di P2, con un profitto massimo di $2.500.

Funzionalità avanzate in Pyomo

Nella sezione precedente abbiamo visto quanto sia facile implementare un problema end-to-end con Pyomo. Tuttavia, la maggior parte dei problemi reali non è così semplice da risolvere.

In questa sezione presento alcune funzionalità avanzate che puoi usare per scenari più complessi.

1. Ottimizzazione non lineare

L’ottimizzazione non lineare minimizza o massimizza una funzione obiettivo non lineare soggetta a vincoli non lineari. Vediamo un esempio in cui minimizziamo la somma delle differenze al quadrato con un vincolo circolare.

Descrizione del problema

Minimizza l’obiettivo: Z = (x - 1)² + (y - 2)²

Soggetto a:

• x² + y² ≤ 4
• x, y ≥ 0

In Pyomo possiamo definire le variabili decisionali x e y con limiti a partire da 0 per garantire la non negatività. La funzione obiettivo è scritta come somma di scarti quadratici da punti specifici e il vincolo assicura che la soluzione ricada entro un cerchio di raggio 2.

In questo caso, il solver IPOPT è adatto grazie alla sua capacità di risolvere problemi di ottimizzazione non lineare:

import pyomo.environ as pyo
model = pyo.ConcreteModel()
# Define variables with lower bounds
model.x = pyo.Var(bounds=(0, None))
model.y = pyo.Var(bounds=(0, None))
# Objective function: minimize (x - 1)² + (y - 2)²
model.obj = pyo.Objective(expr=(model.x - 1)**2 + (model.y - 2)**2, sense=pyo.minimize)
# Constraint: x² + y² ≤ 4 (circle of radius 2)
model.circle = pyo.Constraint(expr=model.x**2 + model.y**2 <= 4)
solver = pyo.SolverFactory('ipopt')
result = solver.solve(model)
print('Optimal x:', pyo.value(model.x))
print('Optimal y:', pyo.value(model.y))
print('Minimum Z:', pyo.value(model.obj))Offerto da

2. Programmazione mista intera (MIP)

La programmazione mista intera si usa quando alcune variabili decisionali sono intere (spesso binarie) mentre altre sono continue. È preziosa per problemi decisionali come localizzazione di impianti e pianificazione della produzione.

Descrizione del problema

Un’azienda deve decidere se aprire magazzini nelle località A, B e C. L’obiettivo è minimizzare il costo totale, che include i costi fissi di apertura e i costi di trasporto.

Iniziamo inizializzando i dati, inclusi i costi fissi di apertura dei magazzini, i costi di trasporto, i limiti di capacità e la domanda totale: 

locations = ['A', 'B', 'C']
FixedCost = {'A': 1000, 'B': 1200, 'C': 1500}
TransportCost = {'A': 5, 'B': 4, 'C': 6}
Capacity = {'A': 100, 'B': 80, 'C': 90}
Demand = 150
model = pyo.ConcreteModel()
# Binary variable: 1 if warehouse is open, 0 otherwise
model.y = pyo.Var(locations, domain=pyo.Binary)
# Continuous variable: amount of goods transported
model.x = pyo.Var(locations, domain=pyo.NonNegativeReals)
model.cost = pyo.Objective(
    expr=sum(FixedCost[i] * model.y[i] + TransportCost[i] * model.x[i] for i in locations),
    sense=pyo.minimize
)
# Demand constraint
model.demand = pyo.Constraint(expr=sum(model.x[i] for i in locations) >= Demand)
# Capacity constraints
def capacity_rule(model, i):
    return model.x[i] <= Capacity[i] * model.y[i]
model.capacity = pyo.Constraint(locations, rule=capacity_rule)
solver = pyo.SolverFactory('cbc')
result = solver.solve(model)
for i in locations:
    print(f"Warehouse {i}: Open={pyo.value(model.y[i])}, Transported={pyo.value(model.x[i])}")
print('Minimum Total Cost:', pyo.value(model.cost))Offerto da

Il modello include due tipi di variabili decisionali: una variabile binaria y che rappresenta se un magazzino è aperto (1 se aperto, 0 altrimenti) e una variabile continua x che rappresenta la quantità di merci trasportate da ciascun magazzino. 

La funzione obiettivo somma i costi fissi e i costi di trasporto di ciascun magazzino e minimizza il totale. I vincoli assicurano che le merci trasportate soddisfino la domanda e che la capacità di ciascun magazzino non venga superata se questo è aperto. 

3. Gestire obiettivi multipli

A volte i problemi di ottimizzazione coinvolgono obiettivi multipli potenzialmente in conflitto, come massimizzare il profitto minimizzando al contempo l'impatto ambientale. Un approccio comune è il metodo della somma pesata, in cui a ciascun obiettivo è assegnato un peso per bilanciarne l'importanza.

Descrizione del problema

Miriamo a massimizzare il profitto minimizzando l’impatto ambientale:

• Profitto: Z₁ = 3x + 5y
• Impatto ambientale: Z₂ = 2x + y

Possiamo combinare questi obiettivi usando pesi w1=0.6, w2=0.4, dove l’obiettivo totale diventa una somma pesata:

w1 = 0.6
w2 = 0.4
model.obj = pyo.Objective(
    expr=w1 * (3 * model.x + 5 * model.y) - w2 * (2 * model.x + model.y),
    sense=pyo.maximize
)Offerto da

In questo obiettivo combinato, massimizziamo il profitto minimizzando l’impatto ambientale tramite l’aggiustamento dei pesi.

4. Utilizzare fonti dati esterne

Quando si lavora con dataset ampi, è spesso utile importare dati da fonti esterne come file CSV. Pyomo lavora bene con Pandas per leggere e utilizzare dati esterni.

Possiamo leggere un file CSV con Pandas e usare i dati per inizializzare insiemi e parametri nel modello:

import pandas as pd
data = pd.read_csv('parameters.csv')
# Define set from CSV data
model.I = pyo.Set(initialize=data['index'].unique())
# Define parameter initialized from CSV data
param_dict = data.set_index('index')['value'].to_dict()
model.param = pyo.Param(model.I, initialize=param_dict)Offerto da

Suggerimenti e best practice per usare Pyomo

Quando lavori con Pyomo, è importante mantenere i modelli efficienti, ben documentati e facili da diagnosticare. 

1. Debug e troubleshooting

Durante la costruzione di modelli di ottimizzazione in Pyomo, è comune imbattersi in problemi come soluzioni non ammissibili, errori del solver o risultati errati. Ecco alcune best practice per il debug:

• Controlla i vincoli: se il tuo modello non produce una soluzione ammissibile, rivedi i vincoli. Vincoli troppo stretti possono rendere il problema inammissibile. Usa il metodo .display() di Pyomo per stampare i valori di variabili e vincoli e verificare che si comportino come previsto.
• Output del solver: abilita log dettagliati del solver passando tee=True quando chiami il metodo solve(). Questo può dare indicazioni su dove il solver potrebbe incontrare difficoltà, ad esempio variabili non limitate o inammissibilità.
• Prova prima modelli semplici: quando lavori con modelli complessi, prova una versione semplificata. Questo può aiutare a isolare i problemi senza il peso di un modello completamente specificato.

La risoluzione dei problemi è molto più semplice se la affronti in modo sistematico, analizzando vincoli, funzione obiettivo e feedback del solver.

2. Efficienza di modellazione

I problemi di ottimizzazione possono diventare onerosi computazionalmente con l’aumentare delle dimensioni del modello. Per garantire un modello efficiente, considera i seguenti suggerimenti:

• Sfrutta la sparsità: evita di iterare su indici non necessari quando definisci vincoli o obiettivi. Sfruttare la sparsità del problema riduce i tempi di calcolo.
• Variabili binarie vs. continue: quando possibile, riduci il numero di variabili binarie o intere. Le variabili continue sono più semplici da gestire per i solver e portano a soluzioni più rapide.
• Formulazione dei vincoli: mantieni i vincoli il più semplici possibile, sia nella forma matematica che nell’implementazione. Evita non linearità non necessarie e scomponi i vincoli complessi in parti più piccole e gestibili.

Modelli efficienti si risolvono più velocemente e sono più facili da debug e mantenere.

3. Documentazione e manutenzione

Mantenere modelli Pyomo ben documentati è una buona pratica per l’uso a lungo termine e la collaborazione. Una buona documentazione rende anche più facile rivedere e aggiornare i modelli nel tempo:

• Usa commenti inline: aggiungi sempre commenti per spiegare lo scopo di variabili, vincoli e funzione obiettivo. Questo è particolarmente importante nei modelli di ottimizzazione, dove la logica potrebbe non essere immediata.
• Modularizza il codice: suddividi il modello in sezioni logiche o persino in funzioni separate. Questo approccio modulare migliora la leggibilità e rende più facile il debug e la modifica di parti specifiche del modello.
• Tieni traccia delle modifiche: conserva una cronologia delle versioni del modello, specialmente se evolve. Usa strumenti di versionamento come Git per tracciare i cambiamenti e assicurarti che aggiornamenti o miglioramenti siano sempre rintracciabili.

Una corretta documentazione e codice strutturato renderanno i tuoi modelli Pyomo più accessibili ai futuri collaboratori e più facili da scalare o modificare con l’evolvere dei requisiti.

Conclusione

Pyomo è uno strumento potente e flessibile per costruire e risolvere modelli di ottimizzazione in Python. In tutto il tutorial, abbiamo visto come Pyomo consenta di modellare vari problemi di ottimizzazione, dalla programmazione lineare a quella non lineare e mista intera. 

Grazie alla sintassi intuitiva e all’integrazione con i solver, Pyomo rende la formulazione e la risoluzione di problemi reali accessibile sia ai principianti sia agli utenti avanzati.

Se vuoi approfondire la risoluzione di problemi reali con l’ottimizzazione, dai un’occhiata al corso gratuito di DataCamp Introduction to Optimization in Python!