📖 Se non hai letto il articolo precedentenon preoccuparti. Anche la formulazione matematica lo è dichiarato (ma non derivato) qui, con ciascun componente del modello accanto alla relativa implementazione del codice. Se non capisci da dove vengono o significano le cose, leggi l’articolo di “sprint 2“, e se desideri più contesto, leggi l’articolo per “sprint 1“.
Python viene utilizzato poiché è il linguaggio migliore nella scienza dei dati e Piomo poiché è una delle migliori librerie (open source) che gestiscono efficacemente modelli su larga scala.
In questa sezione, li esaminerò ciascuno componente del modello definito nella formulazione, e spiegare come viene tradotto nel codice Pyomo. Ho cercato di non lasciare alcuna lacuna, ma se la pensi diversamente, lascia una domanda nei commenti.
Disclaimer: Si prevede che il lettore target sia nuovo a Pyomo e persino alla modellazione, in modo da ridurre il carico cognitivo, un’implementazione concisa e semplice ha la priorità rispetto alle migliori pratiche di programmazione. L’obiettivo ora è insegnare la modellazione di ottimizzazione, non l’ingegneria del software. Il codice viene migliorato in modo incrementale nelle iterazioni future man mano che questa prova di concetto si evolve in un MVP più sofisticato.
1.1. Installazione delle dipendenze
Per chi ha fretta
Installa (o assicurati di aver già installato) le librerie pyomo, networkx E pandase il pacchetto glpk.
📝 Il pacchetto
glpkcontiene il Risolutore GLPKche è un (open source) risolutore esterno utilizzeremo per ottimizzare i modelli che creiamo. Pyomo è abituato creare modelli di problemi e passarli a GLPK, che eseguirà gli algoritmi che eseguono il processo di ottimizzazione stesso. GLPK restituisce quindi le soluzioni all’oggetto del modello Pyomo, che vengono archiviate come attributi del modello, in modo che possiamo usarle comodamente senza uscire da Python.
Il modo consigliato per installare GLPK è tramite Conda in modo che Pyomo possa trovare facilmente il risolutore GLPK. Per installare tutte le dipendenze insieme in una volta sola, esegui:
conda install -y -c conda-forge pyomo=6.5.0 pandas networkx glpk
Per persone organizzate
Consiglio di crearne uno separato ambiente virtuale in cui installare tutte le librerie necessarie per seguire gli articoli di questa serie. Copia questo testo
name: ttp # traveling tourist problem
channels:
- conda-forge
dependencies:
- python=3.9
- pyomo=6.5.0
- pandas
- networkx
- glpk # external solver used to optimize models
- jupyterlab # comment this line if you won't use Jupyter Lab as IDE
e salvalo in un file YAML denominato environment.yml. Apri un prompt di Anaconda nella stessa posizione ed esegui il comando
conda env create --file environment.yml
Dopo alcuni minuti, viene creato l’ambiente con tutte le dipendenze installate al suo interno. Correre conda activate ttp per “entrare” nell’ambiente, avvia Jupyter Lab (eseguendo jupyter lab nel terminale) e inizia a programmare!
1.2. La matematica diventa codice
Innanzitutto, assicuriamoci che Pyomo possa trovare il risolutore GLPK
### ===== Code block 3.1 ===== ###
import pandas as pd
import pyomo.environ as pyo
import pyomo.version
import networkx as nxsolver = pyo.SolverFactory("glpk")
solver_available = solver.available(exception_flag=False)
print(f"Solver '{solver.name}' available: {solver_available}")
if solver_available:
print(f"Solver version: {solver.version()}")
print("pyomo version:", pyomo.version.__version__)
print("networkx version:", nx.__version__)
Solver 'glpk' available: True
Solver version: (5, 0, 0, 0)
pyomo version: 6.5
networkx version: 2.8.4
⛔ Se hai ricevuto il messaggio
'glpk' available: Falseil risolutore non è stato installato correttamente. Prova uno di questi per risolvere il problema:– ripetere attentamente i passaggi di installazione
– correre
conda install -y -c conda-forge glpknell’ambiente di base (predefinito).– prova a installare un risolutore diverso che funzioni per te
Quindi leggere il file CSV dei dati sulla distanza
### ===== Code block 3.2 ===== ###
path_data = (
"https://raw.githubusercontent.com/carlosjuribe/"
"traveling-tourist-problem/main/"
"Paris_sites_spherical_distance_matrix.csv"
)
df_distances = pd.read_csv(path_data, index_col="site")df_distances
Ora entriamo nella “fase 4” del (flusso di lavoro Agile Operations Research), contrassegnato dal blocco verde di seguito:
Il compito è prendere il modello matematico creato in precedenza e implementarlo nel codice esattamente come è stato definito matematicamente.
👁️ Possiamo creare tutti gli oggetti Python che vogliamo se ciò rende più semplice l’implementazione del modello, ma non ci è consentito modificare il modello sottostante in alcun modo, mentre lo stiamo codificando. Ciò causerebbe la non sincronizzazione del modello matematico e del modello computerizzatorendendo così piuttosto difficile il debug successivo del modello.
Istanziamo un modello Pyomo vuoto, in cui memorizzeremo i componenti del modello come attributi:
model = pyo.ConcreteModel("TSP")
1.2.1. Imposta
Per creare l’insieme dei siti 𝕊 = {Louvre, Tour Eiffel, …, hotel}, estraiamo i loro nomi dall’indice del dataframe e lo utilizziamo per creare un Pyomo Set di nome sites:
### ===== Code block 3.3 ===== ###
list_of_sites = df_distances.index.tolist()model.sites = pyo.Set(initialize=list_of_sites,
domain=pyo.Any,
doc="set of all sites to be visited (𝕊)")
Per creare l’insieme derivato
Conserviamo il filtro 𝑖 ≠ 𝑗 all’interno di a regola di costruzione (la funzione Python _rule_domain_arcs) e passare questa regola al file filter parola chiave durante l’inizializzazione del file Set. Tieni presente che questo filtro verrà applicato al prodotto incrociato dei siti (𝕊 × 𝕊) e filtrerà i membri del prodotto incrociato che non soddisfano la regola.
### ===== Code block 3.4 ===== ###
def _rule_domain_arcs(model, i, j):
""" All possible arcs connecting the sites (𝔸) """
# only create pair (i, j) if site i and site j are different
return (i, j) if i != j else None rule = _rule_domain_arcs
model.valid_arcs = pyo.Set(
initialize=model.sites * model.sites, # 𝕊 × 𝕊
filter=rule, doc=rule.__doc__)
1.2.2. Parametri
Il parametro
𝐷ᵢⱼ ≔ Distanza tra il sito 𝑖 e il sito 𝑗
viene creato con il costruttore pyo.Paramche prende come il primo argomento (posizionale) il dominio 𝔸 (model.valid_arcs) che lo indicizza e come il parola chiave discussione initialize un’altra regola di costruzione, _rule_distance_between_sitesche viene valutato per ogni coppia (𝑖, 𝑗) ∈ 𝔸. In ogni valutazione, il valore numerico della distanza viene recuperato dal dataframe df_distancese collegato internamente alla coppia (𝑖, 𝑗):
### ===== Code block 3.5 ===== ###
def _rule_distance_between_sites(model, i, j):
""" Distance between site i and site j (𝐷𝑖𝑗) """
return df_distances.at(i, j) # fetch the distance from dataframerule = _rule_distance_between_sites
model.distance_ij = pyo.Param(model.valid_arcs,
initialize=rule,
doc=rule.__doc__)
1.2.3. Variabili decisionali
Poiché 𝛿ᵢⱼ ha lo stesso “dominio dell’indice” di 𝐷ᵢⱼ, il modo di costruire questo componente è molto simile, con l’eccezione che qui non è necessaria alcuna regola di costruzione.
model.delta_ij = pyo.Var(model.valid_arcs, within=pyo.Binary,
doc="Whether to go from site i to site j (𝛿𝑖𝑗)")
Il primo argomento posizionale di pyo.Var è riservato per il suo insieme di indici 𝔸 e il “tipo” di variabile è specificato con l’argomento parola chiave within. Con “tipo di variabile” intendo il gamma di valori che la variabile può assumere. Qui, l’intervallo di 𝛿ᵢⱼ è solo 0 e 1, quindi è di tipo binario. Matematicamente scriveremmo 𝛿ᵢⱼ ∈ {0, 1}, ma invece di creare vincoli separati per indicarlo, possiamo indicarlo direttamente in Pyomo impostando within=pyo.Binary nel momento in cui creiamo la variabile.
1.2.4. Funzione obiettivo
Per costruire una funzione obiettivo possiamo “memorizzare” l’espressione in una funzione che verrà utilizzata come regola di costruzione per l’obiettivo. Questa funzione accetta solo un argomento, il modello, che viene utilizzato per recuperare qualsiasi componente del modello necessario per creare l’espressione.
### ===== Code block 3.6 ===== ###
def _rule_total_distance_traveled(model):
""" total distance traveled """
return pyo.summation(model.distance_ij, model.delta_ij)rule = _rule_total_distance_traveled
model.obj_total_distance = pyo.Objective(rule=rule,
sense=pyo.minimize,
doc=rule.__doc__)
Osservare il parallelismo tra l’espressione matematica e l’istruzione return della funzione. Precisiamo che si tratta di un problema di minimizzazione con il sense parola chiave.
1.2.5. Vincoli
Se ricordi l’articolo precedente, un modo conveniente per imporre che ciascun sito venga visitato una sola volta è imporre che ogni sito venga “entrato” una volta e “uscito” una volta, contemporaneamente.
Ogni sito viene inserito una sola volta
def _rule_site_is_entered_once(model, j):
""" each site j must be visited from exactly one other site """
return sum(model.delta_ij(i, j) for i in model.sites if i != j) == 1rule = _rule_site_is_entered_once
model.constr_each_site_is_entered_once = pyo.Constraint(
model.sites,
rule=rule,
doc=rule.__doc__)
Ogni sito viene chiuso una sola volta
def _rule_site_is_exited_once(model, i):
""" each site i must departure to exactly one other site """
return sum(model.delta_ij(i, j) for j in model.sites if j != i) == 1rule = _rule_site_is_exited_once
model.constr_each_site_is_exited_once = pyo.Constraint(
model.sites,
rule=rule,
doc=rule.__doc__)
1.2.6. Controllo finale del modello
L’implementazione del modello è terminata. Per vedere come appare nel suo complesso, dovremmo eseguire model.pprint()e navigare un po’ per vedere se abbiamo tralasciato qualche dichiarazione o commesso qualche errore.
Per avere un’idea delle dimensioni del modello, stampiamo il numero di variabili e vincoli di cui è composto:
def print_model_info(model):
print(f"Name: {model.name}",
f"Num variables: {model.nvariables()}",
f"Num constraints: {model.nconstraints()}", sep="\n- ")print_model_info(model)
#(Out):
# Name: TSP
# - Num variables: 72
# - Num constraints: 18
Avendo meno di 100 vincoli o variabili, questo è un problema di piccole dimensioni e verrà ottimizzato in tempi relativamente brevi dal risolutore.
Fonte: towardsdatascience.com
