Implementazione, Risoluzione e Visualizzazione del Problema del Commesso Viaggiatore con Python

Implementazione, Risoluzione e Visualizzazione del Problema del Commesso Viaggiatore utilizzando Python

Tradurre un modello di ottimizzazione dalla matematica a Python, ottimizzarlo e visualizzare la soluzione per ottenere un feedback rapido sugli errori di modellazione

Foto di John Matychuk su Unsplash

👁️ Questo è l’articolo n. 3 della serie che copre il progetto “Un sistema di supporto alle decisioni intelligenti per il turismo in Python”. Ti incoraggio a dare un’occhiata per avere una panoramica generale dell’intero progetto. Se sei interessato solo a come implementare un modello del TSP in Python, sei comunque nel posto giusto: questo articolo è autosufficiente e ti guiderò attraverso tutti i passaggi, dall’installazione delle dipendenze, all’analisi, al codice, all’interpretazione dei risultati e al debug del modello.

Il presente articolo continua il nostro viaggio esattamente dove si era interrotto sprint 2. Qui, prendiamo il modello matematico che abbiamo formulato nell’articolo precedente e lo implementiamo in Python, utilizzando Pyomo e seguendo buone pratiche. Successivamente, il modello viene ottimizzato e le sue soluzioni vengono visualizzate e analizzate. Per esigenze di didattica, scopriamo che la formulazione iniziale del modello è incompleta, quindi mostro come possiamo derivare, dai primi principi, i vincoli necessari per correggere la formulazione. Questi nuovi vincoli vengono aggiunti al modello Pyomo, e le nuove soluzioni vengono nuovamente analizzate e verificate.

Indice

1. Implementazione del modello in Python utilizzando Pyomo

  • 1.1. Installazione delle dipendenze
  • 1.2. La matematica diventa codice

2. Risoluzione e convalida del modello

  • 2.1. Risoluzione del modello
  • 2.2. Visualizzazione dei risultati
  • 2.3. Analisi dei risultati

3. Correzione della formulazione

  • 3.1. L’idea motivante
  • 3.2. Esprimere l’idea motivante come implicazioni logiche
  • 3.3. Formulare le implicazioni logiche come vincoli lineari
  • 3.4. Determinare un valore adeguato per “M grande”

4. Implementazione e verifica della nuova formulazione

  • 4.1. Aggiornamento del modello Pyomo
  • 4.2. Visualizzazione della soluzione del modello aggiornato

5. Conclusioni (per il prossimo sprint)

📖 Se non hai letto l’articolo precedente, non preoccuparti. La formulazione matematica è anche riportata (ma non derivata) qui, con ogni componente del modello accanto alla sua implementazione in codice. Se non comprendi da dove provengono o cosa significano le cose, ti preghiamo di leggere l’articolo di “sprint 2” e, se desideri ulteriori informazioni sul contesto, leggere l’articolo per “sprint 1”.

1. Implementazione del modello in Python utilizzando Pyomo

Python è utilizzato in quanto è il linguaggio principale per la scienza dei dati, e Pyomo in quanto è una delle migliori librerie (open-source) che si occupano in modo efficace di modelli su larga scala.

In questa sezione, andrò attraverso ogni componente del modello definito nella formulazione, e spiegherò come viene tradotto in codice Pyomo. Ho cercato di non lasciare vuoti, ma se pensi diversamente, lascia una domanda nei commenti.

Avviso: Il lettore di destinazione si presume che sia nuovo a Pyomo, e anche alla modellazione, quindi per ridurre il loro carico cognitivo, viene data priorità a una implementazione concisa e semplice piuttosto che alle migliori pratiche di programmazione. L’obiettivo ora è insegnare la modellazione dell’ottimizzazione, non l’ingegneria del software. Il codice viene migliorato incrementalmente nelle iterazioni future, man mano che questa prova di concetto si sviluppa in un MVP più sofisticato.

1.1. Installazione delle dipendenze

Per le persone di fretta

Installa (o assicurati di aver già installato) le librerie pyomo, networkx e pandas, e il pacchetto glpk.

📝 Il pacchetto glpk contiene il risolutore GLPK, che è un (open source) risolutore esterno che useremo per ottimizzare i modelli che creiamo. Pyomo viene utilizzato per creare modelli di problemi e passarli a GLPK, che eseguirà gli algoritmi che effettuano il processo di ottimizzazione stesso. GLPK quindi restituisce le soluzioni all’oggetto del modello Pyomo, che vengono memorizzate come attributi del modello, in modo da poterle utilizzare 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 tutte in una volta, esegui:

conda install -y -c conda-forge pyomo=6.5.0 pandas networkx glpk

Per le persone organizzate

Raccomando di creare un ambiente virtuale separato in cui installare tutte le librerie necessarie per seguire gli articoli di questa serie. Copia questo testo

name: ttp  # traveling tourist problemchannels:  - conda-forgedependencies:  - python=3.9  - pyomo=6.5.0  - pandas  - networkx  - glpk  # risolutore esterno utilizzato per ottimizzare i modelli  - jupyterlab  # commenta questa riga se non userai Jupyter Lab come IDE

e salvalo in un file YAML chiamato environment.yml. Apri un prompt di Anaconda nella stessa posizione e esegui il comando

conda env create --file environment.yml

Dopo alcuni minuti, l’ambiente viene creato con tutte le dipendenze installate al suo interno. Esegui conda activate ttp per “entrare” nell’ambiente, avvia Jupyter Lab (eseguendo jupyter lab nel terminale) e inizia a scrivere il codice!

1.2. La matematica diventa codice

Innanzitutto, assicuriamoci che il risolutore GLPK sia rilevabile da Pyomo

### =====  Blocco di codice 3.1  ===== ###import pandas as pdimport pyomo.environ as pyoimport pyomo.versionimport networkx as nxsolver = pyo.SolverFactory("glpk")solver_available = solver.available(exception_flag=False)print(f"Solver '{solver.name}' disponibile: {solver_available}")if solver_available:    print(f"Solver version: {solver.version()}")print("Versione di pyomo:", pyomo.version.__version__)print("Versione di networkx:", nx.__version__)

Risolutore 'glpk' disponibile: TrueVersione del risolutore: (5, 0, 0, 0)Versione di pyomo: 6.5Versione di networkx: 2.8.4

⛔ Se hai ricevuto il messaggio 'glpk' disponibile: False, il risolutore non si è installato correttamente. Prova una di queste soluzioni per risolvere il problema:

– ripeti i passaggi di installazione con attenzione

– esegui conda install -y -c conda-forge glpk nell’ambiente base (predefinito)

– prova a installare un diverso risolutore che funzioni per te

Poi leggi il file CSV dei dati della distanza

### =====  Blocco di codice 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 [workflow di Agile Operations Research], contrassegnata come il blocco verde qui sotto:

Figura 1. Workflow minimalista per la risoluzione dei problemi in OR. 4a fase: modello informatico (Immagine dell'autore)

Il compito è prendere il modello matematico creato in precedenza e implementarlo nel codice esattamente come è stato definito matematicamente.

👁️ Ci è permesso creare quanti oggetti Python vogliamo se ciò semplifica l’implementazione del modello, ma non possiamo modificare in alcun modo il modello sottostante mentre lo implementiamo. Causerebbe una discrepanza tra il modello matematico e il modello informatico, rendendo piuttosto difficile il debug del modello in seguito.

Istanziamo un modello Pyomo vuoto, in cui memorizzeremo i componenti del modello come attributi:

model = pyo.ConcreteModel("TSP")

1.2.1. Insiemi

Per creare l’insieme dei siti 𝕊 = {Louvre, Tour Eiffel, …, hotel}, estraiamo i loro nomi dall’indice del dataframe e li utilizziamo per creare un Set Pyomo chiamato sites:

### =====  Blocco di codice 3.3  ===== ###list_of_sites = df_distances.index.tolist()model.sites = pyo.Set(initialize=list_of_sites,                       domain=pyo.Any,                       doc="insieme di tutti i siti da visitare (𝕊)")

Per creare l’insieme derivato

Espressione 3.1. Insieme derivato dei possibili archi del tour (traiettorie da sito a sito).

Memorizziamo il filtro 𝑖 ≠ 𝑗 all’interno di una regola di costruzione (la funzione Python _rule_domain_arcs), e passiamo questa regola al parametro filter durante l’inizializzazione del Set. Nota che questo filtro verrà applicato al prodotto cartesiano dei siti (𝕊 × 𝕊), e filtrerà i membri del prodotto cartesiano che non soddisfano la regola.

### =====  Blocco di codice 3.4  ===== ###def _rule_domain_arcs(model, i, j):    """ Tutti gli archi possibili che collegano i siti (𝔸) """    # crea solo la coppia (i, j) se il sito i e il sito j sono diversi    return (i, j) if i != j else None rule = _rule_domain_arcsmodel.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.Param, che prende come il primo (posizionale) argomento il dominio 𝔸 (model.valid_arcs) che lo indicizza, e come argomento chiave initialize un’altra regola di costruzione, _rule_distance_between_sites, che viene valutata per ogni coppia (𝑖, 𝑗) ∈ 𝔸. In ogni valutazione, il valore numerico della distanza viene recuperato dal dataframe df_distances, e collegato internamente alla coppia (𝑖, 𝑗):

### =====  Blocco di codice 3.5  ===== ###def _regola_distanza_tra_siti(modello, i, j):    """ Distanza tra il sito i e il sito j (𝐷𝑖𝑗) """    return df_distanze.at[i, j]  # recupera la distanza dal dataframe
regola = _regola_distanza_tra_siti
modello.distanza_ij = pyo.Param(modello.archi_validi,                               initialize=regola,                               doc=regola.__doc__)

1.2.3. Variabili di decisione

Poiché 𝛿ᵢⱼ ha lo stesso “dominio degli indici” di 𝐷ᵢⱼ, il modo di costruire questo componente è molto simile, con l’eccezione che qui non è necessaria una regola di costruzione.

Espressione 3.2. Variabili di decisione binarie
modello.delta_ij = pyo.Var(modello.archi_validi, within=pyo.Binary,                          doc="Se passare dal sito i al sito j (𝛿𝑖𝑗)")

Il primo argomento posizionale di pyo.Var è riservato per il suo insieme di indici 𝔸, e il “tipo” di variabile è specificato con l’argomento chiave within. Con “tipo di variabile” intendo l’intervallo di valori che la variabile può assumere. Qui, l’intervallo di 𝛿ᵢⱼ è semplicemente 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 al momento della creazione della variabile.

1.2.4. Funzione obiettivo

Espressione 3.3. La funzione obiettivo da minimizzare: distanza totale del tour

Per costruire una funzione obiettivo possiamo “memorizzare” l’espressione in una funzione che verrà utilizzata come regola di costruzione per l’obiettivo. Questa funzione prende solo un argomento, il modello, che viene utilizzato per recuperare eventuali componenti del modello necessarie per costruire l’espressione.

### =====  Blocco di codice 3.6  ===== ###def _regola_distanza_totale_percorsa(modello):    """ distanza totale percorsa """    return pyo.summation(modello.distanza_ij, modello.delta_ij)regola = _regola_distanza_totale_percorsamodello.obiettivo_distanza_totale = pyo.Objective(regola=regola,                                          sense=pyo.minimize,                                          doc=regola.__doc__)

Osserva il parallelismo tra l’espressione matematica e l’istruzione di ritorno della funzione. Specifichiamo che si tratta di un problema di minimizzazione con la parola chiave sense.

1.2.5. Vincoli

Se ti ricordi dal precedente articolo, un modo conveniente per garantire che ogni sito venga visitato solo una volta è garantire che ogni sito venga “entrato” solo una volta e “uscito” solo una volta, contemporaneamente.

Ogni sito viene visitato solo una volta

Espressione 3.4. Insieme di vincoli che garantisce che ogni sito venga
def _regola_sito_entrato_una_volta(modello, j):    """ ogni sito j deve essere visitato da esattamente un altro sito """    return sum(modello.delta_ij[i, j] for i in modello.siti if i != j) == 1regola = _regola_sito_entrato_unamodello.vincolo_ogni_sito_entrato_una_volta = pyo.Constraint(                                          modello.siti,                                           rule=regola,                                           doc=regola.__doc__)

Ogni sito viene lasciato una sola volta

Espressione 3.5. Insieme di vincoli che garantiscono che ogni sito venga lasciato solo una volta.
def _rule_site_is_exited_once(model, i):    """ ogni sito i deve partire esattamente verso un altro sito """    return sum(model.delta_ij[i, j] for j in model.sites if j != i) == 1rule = _rule_site_is_exited_oncemodel.constr_each_site_is_exited_once = pyo.Constraint(                                          model.sites,                                           rule=rule,                                           doc=rule.__doc__)

1.2.6. Ispezione finale del modello

L’implementazione del modello è completa. Per vedere com’è nel complesso, dovremmo eseguire model.pprint() e navigare un po’ per vedere se ci sono state delle dichiarazioni omesse o degli errori commessi.

Per avere un’idea delle dimensioni del modello, stampiamo il numero di variabili e vincoli che lo compongono:

def print_model_info(model):    print(f"Nome: {model.name}",           f"Numero di variabili: {model.nvariables()}",          f"Numero di vincoli: {model.nconstraints()}", sep="\n- ")print_model_info(model)#[Out]:# Nome: TSP#  - Numero di variabili: 72#  - Numero di vincoli: 18

Con meno di 100 vincoli o variabili, si tratta di un problema di piccole dimensioni e verrà ottimizzato relativamente velocemente dal risolutore.

2. Risoluzione e convalida del modello

2.1. Risoluzione del modello

Il passo successivo nel [flusso di lavoro del RAO] è ottimizzare il modello e ispezionare le soluzioni:

res = solver.solve(model)  # ottimizzazione del modelloprint(f"Soluzione ottimale trovata: {pyo.check_optimal_termination(res)}")# [Out]: Soluzione ottimale trovata: True

Ottima notizia! Il risolutore ha trovato una soluzione ottimale a questo problema! Controlliamola per sapere quale percorso seguire quando arriveremo a Parigi!

Per un controllo molto veloce, possiamo eseguire model.delta_ij.pprint(), che stampa tutti i valori (ottimali) delle variabili 𝛿ᵢⱼ, che possono essere 0 o 1:

Figura 3.1. Estratto dei valori delle variabili decisionali come stampato dal modello (Immagine dell'Autore)

È difficile visualizzare un percorso semplicemente guardando un elenco di archi scelti, senza parlare di analizzarlo per convalidare che abbiamo formulato correttamente il modello.

Per capire veramente la soluzione, dobbiamo visualizzarla.

2.2. Visualizzazione dei risultati

Un’immagine vale più di mille registri

Dal momento che stiamo trattando nodi e archi, il modo più semplice per visualizzare la soluzione è rappresentarla come un grafico. Ricorda che questa è una prova di concetto, quindi un feedback rapido ed efficace prevale sulla bellezza. Le visualizzazioni più approfondite possono essere fatte dopo aver ottenuto un MVP funzionante. Per ora, scriviamo alcune funzioni di aiuto per rappresentare le soluzioni in modo efficiente. La

La funzione extract_solution_as_arcs riceve in input un modello Pyomo risolto ed estrae gli “archi scelti” dalla soluzione. Successivamente, la funzione plot_arcs_as_graph memorizza la lista di archi attivi in un oggetto Graph per un’analisi più semplice e rappresenta quel grafico in modo che l’hotel sia l’unico nodo rosso, per riferimento. Infine, la funzione plot_solution_as_graph chiama le due funzioni precedenti per mostrare la soluzione del modello dato come un grafico.

### =====  Blocco di codice 3.7  ===== ###def estrai_soluzione_come_archi(modello):    """ Estrai la lista degli archi attivi (selezionati) dal modello risolto,    mantenendo solo gli archi i cui variabili binarie sono 1 """    archi_attivi = [(i, j)                   for i, j in modello.archi_validi                   if modello.delta_ij[i, j].valore == 1]    return archi_attividef grafico_archi_come_grafo(tour_come_archi):    """ Prendi una lista di tuple che rappresentano gli archi, convertila     in un grafo networkx e disegnala    """    G = nx.DiGraph()    G.add_edges_from(tour_come_archi)  # memorizza la soluzione come grafo    colori_nodi = ['red' if nodo == 'hotel' else 'skyblue'                    for nodo in G.nodi()]    nx.draw(G, colori_nodo=colori_nodi, con_etichette=True, dimensionefonte=6,             forma_nodo='o', dimensione_frecce=5, stile='solid')def grafico_soluzione_come_grafo(modello):    """ Disegna la soluzione del modello dato come un grafo """    print(f"Distanza totale: {modello.obj_distanza_totale()}")        archi_attivi = estrai_soluzione_come_archi(modello)    grafico_archi_come_grafo(archi_attivi)

Ora possiamo vedere come appare realmente la soluzione:

grafico_soluzione_come_grafo(modello)
Figura 3.2. La soluzione della prima formulazione, mostrando sottotour indesiderati invece di un singolo tour. (Immagine dell'autore)

Bene, questo è ovviamente non ciò che ci aspettavamo! Non c’è un singolo tour che visita tutti i siti e torna all’hotel. È vero che tutti i siti sono visitati, ma solo come parte di piccoli cluster di siti sconnessi. Tecnicamente, sì, i vincoli che abbiamo specificato vengono rispettati attentamente: ogni sito viene visitato solo una volta ed esce solo una volta, ma il risultato complessivo non è un singolo tour, come avevamo intenzione, ma un gruppo di sottotour. Ciò significa che l’assunzione che abbiamo fatto nell’articolo precedente) è sbagliata, quindi qualcos’altro manca nel modello che codifichi l’esigenza che “i sottotour non sono ammessi nella soluzione”.

2.3. Analisi dei risultati

Cosa è andato storto?

Quando la soluzione di un modello non ha senso, c’è una sola spiegazione possibile: il modello è sbagliato¹.

Il solutore ci ha dato ciò che è veramente la soluzione ottimale del modello, ma gli abbiamo dato un modello che non mappa il problema che vogliamo risolvere. Il compito ora è scoprire perché e dove abbiamo commesso un errore. Riflettendo, il candidato ovvio è l’assunzione dubbia che abbiamo fatto negli ultimi due paragrafi della sezione “4.4. Creazione dei vincoli” nell’articolo precedente, dove abbiamo progettato il modello matematico. Noi (sbagliando, come sappiamo ora) abbiamo supposto che la formazione di un singolo tour derivasse naturalmente dai due vincoli che garantiscono che ogni sito sia visitato solo una volta. Ma come abbiamo appena visualizzato, non succede. Perché?

La causa principale dell’errore risiede in ciò che chiamo “senso comune non detto“: conoscenza che abbiamo del mondo che è così ovvia da dimenticare di specificarla in un modello

Sapevamo, implicitamente, che la teletrasportazione non è possibile durante la visita dei siti, ma non l’abbiamo esplicitamente detto al modello. Ecco perché osserviamo quei piccoli sottotour, che collegano alcuni dei siti, ma non tutti. Il modello “pensa” che sia normale teletrasportarsi da un sito all’altro, purché, una volta sul sito, si esca una volta e si entri una volta (vedi nuovamente la Figura 3.2). Se vediamo dei sottotour è solo perché abbiamo detto al modello di minimizzare la distanza del tour, e per caso la teletrasportazione è utile per risparmiare distanza.

Quindi, dobbiamo prevenire la formazione di questi sottotour per ottenere una soluzione realistica. Dobbiamo progettare nuovi vincoli che “dicano al modello” che i sottotour sono vietati, o, equivalente, che la soluzione deve essere un unico percorso. Optiamo per quest’ultimo e deduciamo, dai principi fondamentali, un insieme di vincoli intuitivi che svolgano egregiamente il compito.

3. Correzione della formulazione

Facendo riferimento al [workflow di Agile Operations Research], siamo ora nella fase di riformulazione del modello. Una riformulazione di un modello può riguardare un miglioramento o una correzione. La nostra riformulazione si occuperà di correggere il modello.

Sappiamo cosa vogliamo: forzare la soluzione ad essere un unico percorso, che inizia e termina nel nostro sito iniziale, l’hotel. La sfida è come codificare tale requisito in un insieme di vincoli lineari. Di seguito è riportata un’idea, derivante dalle proprietà di un percorso.

3.1. L’idea motivante

Abbiamo 𝑁 siti da “attraversare”, compreso l’hotel. Poiché partiamo dall’hotel, ci rimangono 𝑁 − 1 siti da visitare. Se teniamo traccia dell'”ordine cronologico” in cui visitiamo questi siti, in modo che la prima destinazione (dopo l’hotel) venga assegnata il numero 1, la seconda destinazione venga assegnata il numero 2 e così via, allora l’ultima destinazione prima di tornare all’hotel verrà assegnata il numero 𝑁 − 1. Se chiamiamo questi numeri utilizzati per tenere traccia dell’ordine delle visite “rango”, allora il pattern che si verifica nel percorso è che il rango di ogni sito (diverso dall’hotel) è sempre di 1 unità superiore al rango del sito che lo precede nel percorso. Se potessimo formulare un insieme di vincoli che imponga tale pattern su qualsiasi soluzione ammissibile, introdurremmo, in parole povere, un requisito di “ordine cronologico” nel modello. E si scopre che possiamo farlo effettivamente.

3.2. Esprimere l’idea motivante come implicazioni logiche

💡 Questo è il “pattern” che vogliamo che qualsiasi soluzione ammissibile soddisfi:

Il rango di ogni sito (diverso dall’hotel) deve essere sempre di 1 unità superiore al rango del sito che lo precede nel percorso

Possiamo esprimere nuovamente questo pattern come un’implicazione logica, così: “il rango del sito 𝑗 deve essere di 1 unità superiore al rango del sito 𝑖 se e solo se 𝑗 viene visitato subito dopo 𝑖, per tutti gli archi (𝑖, 𝑗) che non includono l’hotel 𝐻”. Questa formulazione viene espressa matematicamente come:

Espressione 3.6. Implicazione logica per le variabili di rango: aumentano di 1 con ogni nuovo sito visitato.

dove 𝑟ᵢ sono le nuove variabili che dobbiamo seguire per tenere traccia dell’ordine (ancora sconosciuto) delle visite. Per distinguerle dalle variabili di decisione, chiamiamole “variabili di rango“. Il lato destro si legge come “per tutti 𝑖 e 𝑗 che appartengono all’insieme di tutti i siti escludendo l’hotel”. Per comodità notazionale, definiamo il nuovo insieme 𝕊* per memorizzare tutti i siti tranne l’hotel (denotato da 𝐻):

Espressione 3.7. L'insieme di tutti i siti di interesse: tutti i siti tranne l'hotel.

che ci permette di definire le variabili di rango in modo conciso come:

Espressione 3.8. Definizione delle variabili di rango, definite solo per i siti di interesse.

👁️ È fondamentale che l’hotel non abbia una variabile di rango associata perché sarà contemporaneamente l’origine e la destinazione finale di ogni tour, una condizione che violerebbe lo schema di “variabili di rango sempre crescenti nel tour”. In questo modo, il rango di ogni sito è sempre costretto ad aumentare con ogni nuovo arco preso, il che garantisce che i loop chiusi siano vietati a meno che i loop si chiudano all’unico sito che non ha una variabile di rango: l’hotel

I limiti di 𝑟ᵢ sono derivati dalla sua descrizione: il rango parte da 1 e aumenta in modo monotono fino a quando tutti i siti in 𝕊* sono visitati, terminando quindi a | 𝕊* | (la dimensione dell’insieme dei siti non hotel). Inoltre, permettiamo loro di assumere qualsiasi valore reale positivo:

Espressione 3.9. Limiti e intervallo delle variabili di rango

3.3. Formulazione delle implicazioni logiche come vincoli lineari

La sfida ora è tradurre questa implicazione logica in un insieme di vincoli lineari. Fortunatamente, le disuguaglianze lineari servono anche a scopi di applicazione di implicazioni logiche, non solo a limitazioni di risorse finite.

Un modo per farlo è tramite il cosiddetto metodo del “Big-M”, che consiste nel dichiarare un vincolo in modo tale che, quando la condizione che ti interessa è soddisfatta, il vincolo si applica (diventa attivo), e quando la condizione che ti interessa non è soddisfatta, il vincolo diventa ridondante (diventa inattivo). La tecnica viene chiamata “Big-M” perché fa uso di un valore costante 𝑀 che è sufficientemente grande da rendere il vincolo ridondante per ogni caso. Quando 𝑀 non compare nel vincolo, il vincolo è “attivo” nel senso che sta applicando l’implicazione desiderata.

Ma cosa determina se un vincolo è “attivo” o no? La risposta breve sono i valori delle variabili decisionali stesse a cui il vincolo si applica. Vediamo come funziona.

L’implicazione desiderata è avere 𝑟ⱼ = 𝑟ᵢ + 1 solo quando 𝛿ᵢⱼ = 1. Possiamo sostituire l’1 nell’espressione con 𝛿ᵢⱼ, ottenendo così

Questa è la relazione che vogliamo che valga quando 𝛿ᵢⱼ = 1, ma non quando 𝛿ᵢⱼ = 0. Per “correggere” il caso in cui 𝛿ᵢⱼ = 0, aggiungiamo un termine di ridondanza, 𝑅𝑇, il cui compito è “disattivare il vincolo” solo quando 𝛿ᵢⱼ = 0. Pertanto, questo termine di ridondanza deve includere la variabile 𝛿ᵢⱼ, poiché dipende da essa.

In questo contesto, “disattivare il vincolo” significa “renderlo ridondante”, poiché un vincolo ridondante ha lo stesso effetto di un vincolo non esistente in un modello.

Vediamo come possiamo dedurre l’espressione per RT. L’espressione per 𝑅𝑇(𝛿ᵢⱼ) deve soddisfare queste proprietà:

Espressione 3.10. Proprietà che il termine di ridondanza deve soddisfare per garantire la ridondanza valida.

Per soddisfare il punto (1) abbiamo bisogno che 𝑅𝑇(𝛿ᵢⱼ = 1) = 0, e quindi l’espressione per 𝑅𝑇 deve contenere il moltiplicatore (1 – 𝛿ᵢⱼ), poiché diventa 0 quando 𝛿ᵢⱼ = 1. Questa forma fa sì che 𝑅𝑇 “svanisca” quando 𝛿ᵢⱼ = 1, o “si riduca” a una costante (chiamiamola M) quando 𝛿ᵢⱼ = 0. Pertanto, un candidato per il termine di ridondanza è

Espressione 3.11. Definizione del

dove 𝑀 dovrebbe essere determinato dai dati del problema (ne parleremo più avanti).

Per soddisfare il punto (2) per tutti i possibili i e j, dobbiamo rendere l’uguaglianza nell’espressione (3.11) una disuguaglianza (il simbolo = diventa ≥), e trovare una costante 𝑀 abbastanza grande (in valore assoluto) in modo che, indipendentemente dai valori 𝑟ⱼ e 𝑟ᵢ assumano, la restrizione sia sempre soddisfatta. È qui che deriva l'”M grande” (big M).

Una volta trovata una costante sufficientemente grande 𝑀, la nostra restrizione di “implicazione logica” assumerà la forma

Espressione 3.12. Restrizioni per l'implicazione che

Introdurre queste restrizioni nel modello forzerà presumibilmente la soluzione ad essere un unico percorso. Ma le restrizioni non avranno l’effetto desiderato a meno che non specifichiamo prima un buon valore per M.

3.4. Dedurre un valore appropriato per il “big M”

Dato che l’obbiettivo è avere 𝑅𝑇(𝛿ᵢⱼ = 0) = 𝑀, possiamo dedurre un valore appropriato per 𝑀 impostando 𝛿ᵢⱼ = 0 nella restrizione generale definita nell’Espressione (3.12):

Espressione 3.13. Deduzione del valore minimo per M.

Per soddisfare 𝑟ⱼ − 𝑟ᵢ ≥ 𝑀 per tutti i siti non hotel 𝑖, 𝑗, abbiamo bisogno che il limite inferiore di 𝑟ⱼ − 𝑟ᵢ sia anch’esso maggiore di 𝑀. Il limite inferiore (LB) di 𝑟ⱼ − 𝑟ᵢ è il valore minimo che 𝑟ⱼ − 𝑟ᵢ può assumere, e può essere ottenuto da

dove 𝑟ᵐⁱⁿ è il rango minimo possibile e 𝑟ᵐᵃˣ il rango massimo possibile. Pertanto, affinché l’ineguaglianza (1) nell’Espressione (3.13) sia vera per tutti i ranghi di tutti i siti, deve valere anche la seguente disuguaglianza:

Grazie a questa disuguaglianza sappiamo il valore minimo che 𝑀 deve assumere affinché il metodo big-M funzioni, ovvero

Espressione 3.14. Limite inferiore per il big-M.

E quali sono i valori di 𝑟ᵐⁱⁿ e 𝑟ᵐᵃˣ? Nella nostra convenzione, abbiamo assegnato al primo sito visitato il rango 1, che è ovviamente il rango minimo (cioè 𝑟ᵐⁱⁿ = 1). Poiché il rango cresce di 1 unità in ogni sito visitato, l’ultimo sito non alberghiero nel tour avrà il rango massimo, pari al numero di tutti i siti non alberghieri. Poiché il numero di siti non alberghieri può variare, abbiamo bisogno di un’espressione generale. Se ricordi, abbiamo definito l’insieme 𝕊* che contiene tutti i siti non alberghieri, quindi il numero che stiamo cercando è la dimensione (cioè il numero di elementi) dell’insieme 𝕊*, che in notazione matematica è espresso come | 𝕊* |. Pertanto, abbiamo dedotto che 𝑟ᵐⁱⁿ = 1 e 𝑟ᵐᵃˣ = | 𝕊* |. Sostituendo nell’Espressione (3.14), arriviamo infine a un valore adeguato per M:

Espressione 3.15. Il valore del big-M, dedotto dai dati del problema.

Poiché 𝕊* avrà sempre più di 2 siti da visitare (altrimenti, non ci sarebbe nessun problema di decisione in primo luogo), il valore “big M”, in questo modello, è sempre un valore “big” negativo.

📝 Valori teorici rispetto ai valori computazionali

Teoricamente, siamo autorizzati a scegliere arbitrariamente valori “più negativi” per 𝑀 rispetto a quello dedotto qui – persino inventare numeri enormi per essere sicuri che siano abbastanza grandi, per evitare questo calcolo – ma questa non è una buona pratica. Se 𝑀 diventa troppo grande (in valore assoluto), può creare problemi di prestazioni negli algoritmi del risolutore o, nel caso peggiore, persino far sì che il risolutore consideri soluzioni infattibili come fattibili. Ecco perché la pratica consigliata è derivare, dai dati del problema, un valore preciso, ma sufficientemente grande per 𝑀.

Ora che abbiamo dedotto un valore adeguato per “big M”, lo memorizzeremo in un nuovo parametro del modello, per un riutilizzo più semplice. Con questo, la serie di vincoli di eliminazione dei sotto-tour è pronta, e la sua “forma completa” è

Espressione 3.16. Vincoli di eliminazione dei sotto-tour.

Per tenere le cose in prospettiva, nota che questo è effettivamente l’equivalente di “vincolo” dell’implicazione logica originale che abbiamo formulato in precedenza nell’Espressione (3.6):

Congratulazioni! Finalmente abbiamo un insieme di vincoli che possono essere aggiunti al nostro modello. Trovarli è stata la parte difficile. Ora, verifichiamo che la loro aggiunta al modello comporti realmente la scomparsa dei sotto-tour.

4. Implementazione e verifica della nuova formulazione

4.1. Aggiunta del modello Pyomo con la nuova formulazione

La revisione e la correzione del modello hanno richiesto l’aggiunta di alcuni altri insiemi, parametri, variabili e vincoli. Aggiungiamo questi nuovi componenti del modello al modello Pyomo, seguendo lo stesso ordine della fase di formulazione iniziale.

4.1.1. Insiemi e parametri

  • Siti di interesse, 𝕊*: insieme di tutti i siti escluso l’hotel:
Espressione 3.17. Definizione dell'insieme di siti di interesse (noti anche come siti non alberghieri)

Gli oggetti Set di Pyomo hanno operazioni compatibili con gli insiemi Python, quindi possiamo fare direttamente la differenza tra un insieme Pyomo e un insieme Python:

model.sites_except_hotel = pyo.Set(    initialize=model.sites - {'hotel'},     domain=model.sites,    doc="Siti di interesse, ovvero tutti i siti tranne l'hotel (𝕊*)")
  • big M, il nuovo parametro per il vincolo di eliminazione dei sotto-tour:

model.M = pyo.Param(initialize=1 - len(model.sites_except_hotel),                    doc="big M per rendere alcuni vincoli ridondanti")

4.1.2. Variabili ausiliarie

  • Variabili di rango, rᵢ: per tenere traccia dell’ordine in cui vengono visitati i siti:

model.rank_i = pyo.Var(    model.sites_except_hotel,  # i ∈ 𝕊* (indice)    within=pyo.NonNegativeReals,  # rᵢ ∈ ℝ₊ (dominio)    bounds=(1, len(model.sites_except_hotel)),  # 1 ≤ rᵢ ≤ |𝕊*|    doc="Rango di ogni sito per tenere traccia dell'ordine di visita")

Nei commenti, puoi vedere come gli elementi della definizione matematica completa delle variabili corrispondono agli argomenti della funzione di dichiarazione delle variabili di Pyomo pyo.Var. Spero che questo ti aiuti a comprendere il valore di avere un modello matematico ben definito prima di iniziare a costruire un modello Pyomo. L’implementazione si svolgerà in modo naturale e gli errori saranno meno probabili.

4.1.3. Vincoli

  • La soluzione deve essere un tour singolo che inizia e termina all’hotel:

def _rule_path_is_single_tour(model, i, j):    """ Per ogni coppia di siti non alberghieri (i, j),     se il sito j è visitato dal sito i, il rango di j deve essere     strettamente maggiore del rango di i. """    if i == j:  # se i siti coincidono, salta la creazione di un vincolo        return pyo.Constraint.Skip        r_i = model.rank_i[i]    r_j = model.rank_i[j]    delta_ij = model.delta_ij[i, j]    return r_j >= r_i + delta_ij + (1 - delta_ij) * model.M# prodotto cartesiano dei siti non alberghieri, per indicizzare il vinconlo non_hotel_site_pairs = model.sites_except_hotel * model.sites_except_hotelrule = _rule_path_is_single_tourmodel.constr_path_is_single_tour = pyo.Constraint(    non_hotel_site_pairs,    rule=rule,     doc=rule.__doc__)

Il modello Pyomo è stato aggiornato. Di quanto è cresciuto dopo l’aggiunta dei vincoli di eliminazione dei sottotour?

print_model_info(model)# [Out]:# Nome: TSP#  - Numero di variabili: 80#  - Numero di vincoli: 74

Siamo passati da 72 a 80 variabili e da 18 a 74 vincoli. Chiaramente, questa formulazione è più pesante in termini di vincoli che di variabili, poiché ha quadruplicato il numero di vincoli che avevamo prima. Questo è il “prezzo che paghiamo”, in generale, per rendere i modelli più “realistici”, poiché il realismo di solito implica, se i dati non cambiano, di limitare il numero di soluzioni ammissibili.

Come sempre, possiamo analizzare la struttura del modello con model.pprint(). Anche se questa “stampa” perde rapidamente valore man mano che il numero di componenti del modello aumenta, può comunque darci una rapida panoramica di ciò di cui è composto il modello e di quanto è grande.

4.2. Plotting della soluzione aggiornata del modello

Risolviamo il modello aggiornato e plottiamo la nuova soluzione. Incrociamo le dita.

res = solver.solve(model)  # ottimizza i modellif solution_found = pyo.check_optimal_termination(res)print(f"Soluzione ottimale trovata: {solution_found}")# [Out]: Soluzione ottimale trovata: Trueif solution_found:    plot_solution_as_graph(model)

Adesso stiamo parlando! Questo è esattamente ciò che ci aspettavamo dopo l’aggiunta dei nuovi vincoli: nessun sottotour si è formato, rendendo il percorso della soluzione un singolo tour adesso.

Nota come, ovviamente, il valore dell’obiettivo di questo modello risolto sia ora di 14.9 km, invece dei 5.7 km inaffidabili che abbiamo ottenuto con il modello incompleto.

👁️ Un disegno di un grafo non è un percorso su una mappa

Notare che questa immagine è solo una possibile rappresentazione grafica di un grafo, non una traiettoria geografica. I cerchi e i collegamenti che vedete non corrispondono a un percorso reale nello spazio geografico (come potrebbero, se non abbiamo utilizzato alcuna informazione geografica nella sua creazione?). È possibile verificarlo eseguendo plot_solution_as_graph(model) più volte: ogni volta che lo si esegue, i nodi assumeranno diverse posizioni. I grafici sono strutture matematiche astratte che collegano “punti” tramite “collegamenti”, rappresentando relazioni di qualsiasi tipo tra entità di qualsiasi tipo. Abbiamo utilizzato un grafo qui per studiare la validità della soluzione, non per visualizzare il vero tour a Parigi. Lo facciamo in [questo articolo].

5. Conclusioni (o pianificazione per lo sprint successivo)

Con questa verifica finale della soluzione, concludiamo che questa versione aggiornata del modello può risolvere qualsiasi istanza del Problema del Commesso Viaggiatore, quindi possiamo considerarlo un proof-of-concept (POC) di successo.

💡 Evolvere le soluzioni, uno sprint alla volta

Questo POC non risolve (ancora) il nostro problema turistico complesso originale, ma risolve il problema minimo a valore proposto come primo passo verso una soluzione più sofisticata. Quindi, ci avvicina provabilmente (ossia basato su evidenze) a una soluzione di valore del problema più complesso. Con un esempio di lavoro minimo a portata di mano, possiamo valutare meglio cosa deve essere fatto per avanzare nella direzione desiderata e cosa può essere semplificato provvisoriamente fino ad arrivare a una versione più matura della soluzione. Mentre abbiamo qualcosa di utile in qualsiasi momento, svilupperemo un sistema sempre più prezioso, fino a essere soddisfatti. Questa è l’essenza della “strada agile” verso soluzioni efficaci.

Con la validità di questo approccio dimostrata, dobbiamo espanderlo e perfezionarlo in modo che possa gradualmente includere più elementi del nostro problema originale, con ogni iterazione che fornisce soluzioni incrementali preziose. In questa POC ci siamo concentrati sulla progettazione e formulazione di un modello di base, quindi abbiamo dovuto assumere un insieme fisso di siti e la loro matrice di distanza come dato. Naturalmente, questo è limitante e il prossimo passo dovrebbe essere avere un modello che accetti qualsiasi numero di siti. Per questo, abbiamo bisogno di un modo per generare automaticamente la matrice delle distanze dato un elenco di siti e le loro coordinate geografiche. Questo è l’obiettivo del nostro prossimo sprint.

5.1. Quali sono i prossimi passi

Nel prossimo articolo (sprint 4) lavoreremo su una classe che genera automaticamente una matrice delle distanze da qualsiasi elenco di siti. Questa funzionalità, combinata con il modello appena creato qui, ci permetterà, tra le altre cose, di risolvere molti modelli per diversi input in modo rapido e confrontarli. Inoltre, generalizzare la soluzione in questo modo ci faciliterà la vita in futuro quando faremo un po’ di analisi di sensibilità e scenario nei prossimi sprint. Inoltre, mentre aggiorniamo questa prova concettuale a uno stato di “MVP”, inizieremo a utilizzare codice orientato agli oggetti per mantenere tutto ben organizzato e pronto per l’estendibilità. Non perdere il ritmo e inizia subito!

Note a piè di pagina

  1. Effettivamente, c’è un’altra causa di un risultato errato: i dati con cui il modello viene alimentato possono anche essere errati, non solo la formulazione del modello. Ma, in un certo senso, se pensi a un “modello” come a un “istanza del modello”, cioè un modello con dati concreti, allora il modello sarà ovviamente errato se i dati sono sbagliati, che è quello che intendevo con l’affermazione.

Grazie per aver letto, e a presto! 📈😊

Non esitate a seguirmi, a farmi domande, a darmi feedback nei commenti o a contattarmi su LinkedIn.