Il problema della dispersione delle strutture modelli di programmazione mista-intera

Il problema della dispersione delle strutture modelli di programmazione mista-intera' - 'The problem of dispersion in mixed-integer programming model structures

Un tutorial completo di Python sui modelli di dispersione p e maxisum

Foto di Z su Unsplash

In alcuni problemi di posizionamento delle strutture, le strutture devono essere posizionate in modo che l’influenza di una non offuschi o influenzi negativamente le altre. Che sia guidato da motivi di mitigazione del rischio, evitamento della concorrenza o altre considerazioni, padroneggiare le soluzioni a queste sfide è una competenza cruciale nel toolkit della ricerca operativa.

Kuby (1987) ha proposto due formulazioni distinte di programmazione lineare mista (MIP) per questi problemi: il problema di dispersione p e il problema di maxisum. In questo articolo, entrambi saranno implementati utilizzando la libreria Python Pyomo e il solver HiGHS.

Oltre ai due modelli, verranno presentate alcune utili risorse di modellazione. In primo luogo, una strategia per linearizzare il prodotto di variabili binarie nel contesto di MIP, anche considerando gli obiettivi di massimizzazione, non è necessaria in modo esplicito in questo problema. In secondo luogo, una formulazione MIP max-min, che intende massimizzare qualcosa di più piccolo di qualsiasi parametro di un insieme di elementi se l’elemento viene selezionato. Infine, una strategia per risolvere obiettivi multipli con una gerarchia ben definita di priorità che combina elementi dei due modelli.

Coloro interessati che ancora non sono familiari con l’ottimizzazione numerica potrebbero trovare utile dare un’occhiata ai miei precedenti articoli su Programmazione Lineare e il metodo del Branch & Bound prima di procedere con questo.

Come al solito, è possibile trovare il codice completo in questo repository git.

Quale di queste posizioni sceglieresti per posizionare 5 strutture?

Possibili posizioni nel problema di dispersione delle strutture. (Immagine dell'autore).

Il prodotto di variabili binarie

Quando si definiscono gli elementi di base di questo problema, è possibile utilizzare variabili binarie che assumono il valore 1 se una posizione viene selezionata e 0 altrimenti. Indichiamo queste variabili come xᵢ. Supponiamo che le distanze (o qualche altra metrica di dispersione) tra due posizioni siano calcolate in anticipo e le indichiamo come dᵢⱼ. Come potremmo calcolare la dispersione di qualsiasi coppia di strutture selezionate?

In una situazione del genere, è in qualche modo intuitivo formulare il calcolo utilizzando il prodotto delle variabili binarie xᵢ e xⱼ per calcolare la dissimilarità ottenuta quando entrambe sono incluse nella soluzione. Tuttavia, ciò porterebbe a una formulazione quadratica, quindi non sarebbe possibile applicare solver lineari. Fortunatamente, esiste un modo per formulare il prodotto di variabili binarie in MIP utilizzando alcune vincoli.

Consideriamo un grafo diretto G(V, A) e zᵢⱼ una nuova variabile binaria che indica che entrambi i nodi i e j sono selezionati. Se sia i che j non sono selezionati, zᵢⱼ deve essere 0. Questo porta al primo gruppo di vincoli:

Primo gruppo di vincoli nella forma linearizzata del prodotto di variabili binarie. (Immagine dell'autore).

Fino ad ora, anche se sia i che j sono selezionati, zᵢⱼ può assumere il valore di 0. Pertanto, dobbiamo includere un vincolo aggiuntivo tale che quando i e j sono selezionati, zᵢⱼ diventa 1.

Secondo gruppo di vincoli nella forma linearizzata del prodotto di variabili binarie. (Immagine dell'autore).

Quando si massimizza qualcosa proporzionale a zᵢⱼ, come nel problema di maxisum, il secondo gruppo di vincoli presentato non è necessario, poiché non avrebbe senso non calcolare un guadagno proporzionale a zᵢⱼ se è fattibile. Tuttavia, questo può essere utile nella formulazione di altri modelli MIP.

Nella seguente sezione, applichiamo questi concetti al problema del maxisum.

Il modello maxisum

Il problema del maxisum discreto deve definire la posizione di p strutture tra N nodi discreti per massimizzare la somma delle distanze (o la media delle distanze) calcolate tra tutte le coppie di strutture selezionate. Pertanto, consideriamo le strutture come nodi disposti in un grafo diretto G(V, A). Il peso di ogni arco da i a j è la metrica di distanza (dispersione) dᵢⱼ nota in anticipo. Inoltre, consideriamo le variabili binarie xᵢ per indicare se una posizione i è selezionata e zᵢⱼ per indicare se sia i che j sono selezionati.

L’obiettivo può essere formulato come segue:

Oltre ai vincoli per linearizzare il prodotto di variabili binarie presentato nella sezione precedente, è necessario includere un vincolo per garantire che vengano selezionate p posizioni.

Pertanto, possiamo enunciare i vincoli del problema come segue:

Mettiamolo in codice usando Python. Per farlo, dobbiamo iniziare importando le librerie che verranno utilizzate. La libreria numpy verrà utilizzata per i calcoli di algebra lineare e per creare punti di coordinate casuali; le funzioni squareform e pdist da scipy verranno utilizzate per calcolare le distanze date una matrice di coordinate; matplotlib sarà il nostro strumento di visualizzazione; e pyomo sarà il linguaggio di modellazione algebrica (AML) (con il solver HiGHS).

import numpy as npfrom scipy.spatial.distance import squareform, pdistimport matplotlib.pyplot as pltimport pyomo.environ as pyofrom pyomo.contrib.appsi.solvers.highs import Highs

Dobbiamo creare N punti e definire quanti di quelli devono essere selezionati come posizioni delle strutture. Fissando il seme casuale (12) in ogni esecuzione del codice, si dovrebbero ottenere i punti rappresentati nell’introduzione.

# Fissare i semi casualinp.random.seed(12)# Creare punti casualiN = 25p = 5coordinates = np.random.random((N, 2))# Calcolare le distanze tra coppie di puntidistances = squareform(pdist(coordinates))

Ora abbiamo gli elementi necessari per iniziare il nostro modello pyomo.

Ci sono due approcci per modellare un problema in pyomo: modelli astratti e concreti. Nel primo approccio, le espressioni algebriche del problema vengono definite prima che vengano forniti alcuni valori dei dati, mentre nel secondo, l’istanza del modello viene creata immediatamente man mano che i suoi elementi vengono definiti. È possibile trovare ulteriori informazioni su questi approcci nella documentazione della libreria o nel libro di Bynum et al. (2021). In tutto questo articolo, adotteremo la formulazione del modello concreto.

model = pyo.ConcreteModel()

In questo modello, abbiamo due insiemi: nodi e archi. Considerando un grafo diretto completo, esistono archi tra ogni coppia di due nodi eccetto il nodo stesso.

# Insiemi di nodi e archimodel.V = pyo.Set(initialize=range(N))model.A = pyo.Set(initialize=[(i, j) for i in model.V for j in model.V if i != j])

I parametri forniti in anticipo sono il numero di nodi che devono essere selezionati e le distanze tra le coppie di nodi.

# Parametrimodel.d = pyo.Param(model.A, initialize={(i, j): distances[i, j] for (i, j) in model.A})model.p = pyo.Param(initialize=p)

Quindi, includiamo le variabili decisionali.

# Variabili decisionalimodel.x = pyo.Var(model.V, within=pyo.Binary)model.z = pyo.Var(model.A, within=pyo.Binary)

E i vincoli…

# Sono selezionati p nodidef p_selection(model):    return sum(model.x[:]) == model.p# Se il nodo di partenza non è selezionato, l'arco è 0def dispersion_c1(model, i, j):    return model.z[i, j] <= model.x[i]# Se il nodo di arrivo non è selezionato, l'arco è 0def dispersion_c2(model, i, j):    return model.z[i, j] <= model.x[j]# Includi i vincoli nel modellomodel.p_selection = pyo.Constraint(rule=p_selection)model.dispersion_c1 = pyo.Constraint(model.A, rule=dispersion_c1)model.dispersion_c2 = pyo.Constraint(model.A, rule=dispersion_c2)

Infine, dobbiamo creare la funzione obiettivo.

def disp_obj(model):    return sum(model.z[i, j] * model.d[i, j] for (i, j) in model.A)model.obj = pyo.Objective(rule=disp_obj, sense=pyo.maximize)

Ora il modello maxisum è pronto per essere risolto. Dobbiamo quindi istanziare un risolutore compatibile con pyomo per gestirlo. Il risolutore Highs è disponibile in Pyomo (controlla gli import) dalla versione 6.4.3 se il pacchetto highspy è installato. Quindi assicurati di eseguire un pip install highspy .

solver = Highs()solver.solve(model)

E, dopo circa 120 secondi di tempo di calcolo, abbiamo i seguenti risultati:

Risultati del modello maxisum. (Immagine dell'autore).

Nota che, anche se la dispersione totale è massimizzata, due strutture nell’angolo in alto a sinistra sono ancora abbastanza vicine tra loro, il che può essere indesiderabile. Pertanto, in alternativa alla formulazione maxisum abbiamo il modello di p-dispersione in cui massimizziamo la distanza minima tra qualsiasi coppia di nodi selezionati.

Il modello di p-dispersione

Nel modello di p-dispersione, abbiamo bisogno di una variabile decisionale aggiuntiva per calcolare la distanza minima tra tutte le coppie di nodi selezionati, che è il nostro obiettivo da massimizzare. Indichiamo questa variabile con D. Dobbiamo creare un vincolo big M che garantisca che se entrambi i e j sono selezionati, D sia minore o uguale a dᵢⱼ; altrimenti, dobbiamo garantire che D sia illimitato. Se manteniamo la formulazione con zᵢⱼ come il prodotto di variabili binarie, possiamo formulare questo vincolo aggiuntivo come segue.

Vincolo max-min con prodotto di variabili binarie. (Immagine dell'autore).

In questa formulazione, M è un parametro fisso arbitrariamente grande utilizzato per formulare una regola disgiuntiva. Una buona scelta di M dovrebbe essere sufficientemente grande da garantire la fattibilità del vincolo (i, j) per qualsiasi valore di D quando zᵢⱼ è zero, ma il più piccolo possibile in modo che il rilassamento lineare del modello sia simile alla versione intera, il che facilita la convergenza. In questo problema, il valore più grande di dᵢⱼ può essere una scelta interessante.

Anche se questa formulazione funzionerebbe bene per questo modello, è possibile utilizzare un approccio più conciso, in cui le variabili zᵢⱼ non sono nemmeno incluse nel modello.

Vincolo max-min con variabili di nodo. (Immagine dell'autore).

In questa formulazione, sia xᵢ che xⱼ uguali a zero sono una condizione sufficiente per garantire che l’ineguaglianza sia valida per qualsiasi valore di D. E l’obiettivo diventa semplicemente massimizzare D.

Nel codice Python successivo, consideriamo di avere un nuovo modello con gli stessi insiemi e parametri del precedente, così come il gruppo di variabili decisionali x.

# Vincolo max-mindef maxmin_rule(model, i, j):    return model.D <= model.d[i, j] + model.M * (1 - model.x[i]) + model.M * (1 - model.x[j])# Nuovo parametro big Mmodel.M = max(model.d[:, :])# Nuova variabilemodel.D = pyo.Var(within=pyo.NonNegativeReals)# Nuovo vincolomodel.maxmin_rule = pyo.Constraint(model.A, rule=maxmin_rule)# Obiettivomodel.obj = pyo.Objective(expr=model.D, sense=pyo.maximize)

E dopo aver chiamato il risolutore, ci sono voluti meno di 1.2 secondi per ottenere i seguenti risultati.

Risultati del modello di p-dispersione. (Immagine dell'autore).

Che sembra ottimo poiché le posizioni sono ben distribuite nello spazio.

C’è un modo per migliorare questa distribuzione?

Un approccio multicriterio

Ricorda che la funzione obiettivo del modello di dispersione p dipende solo dalla distanza minima tra ogni coppia di nodi selezionati. Pertanto, è possibile ottenere diverse soluzioni utilizzando combinazioni dei due punti che definiscono questa distanza e altri con distanze tra di loro maggiori o uguali all’obiettivo stesso. Possiamo migliorare la nostra soluzione selezionando la migliore di queste alternative? Questo porta all’approccio “multicriterio a due fasi” proposto da Kuby (1987).

Nella prima fase, il modello di dispersione p viene risolto per ottimalità e l’obiettivo corrente viene memorizzato in un parametro d_opt. Quindi, viene risolto un modello di maxisum con un vincolo aggiuntivo che garantisce che D ≥ d_opt, per ottenere, tra le soluzioni ottimali del modello di dispersione p, quella con la massima dispersione totale.

Per farlo in Python, considera di avere un modello di dispersione p istanziato con anche variabili decisionali e vincoli del modello di maxisum.

# D deve essere ottimaledef composed_constr(model):    return model.D >= model.d_opt# Risolvi dispersioni p solver.solve(model)# Nuovo parametromodel.d_opt = pyo.Param(initialize=model.obj())# Disattiva vecchio obiettivomodel.obj.deactivate()# Più variabilimodel.z = pyo.Var(model.A, within=pyo.Binary)# La soluzione non peggiora il valore corrente di Dmodel.composed_constr = pyo.Constraint(rule=composed_constr)# Nuovo obiettivomodel.obj_disp = pyo.Objective(rule=disp_obj, sense=pyo.maximize)solver.solve(model)

E questo è sorprendentemente veloce perché quando il risolutore passa al secondo obiettivo, lo spazio delle alternative ammissibili è sostanzialmente piccolo. In meno di un (ulteriore) secondo, si arriva ai seguenti risultati.

Problema multicriterio: modello di dispersione p seguito da obiettivo di maxisum. (Immagine dell'autore).

Ulteriori letture

Quando i clienti non sono uniformemente distribuiti, le strutture hanno una capacità limitata o forse il numero adeguato di strutture non è noto a priori, probabilmente ti trovi di fronte a un diverso Problema di Localizzazione delle Strutture. Puoi trovare la formulazione di Balinski (1965) implementata in Python utilizzando PuLP in questa incredibile storia di Nicolo Cosimo Albanese.

Ottimizzazione: Problema di Localizzazione delle Strutture con Capacità in Python

Trova il numero e la posizione ottimali dei magazzini per ridurre i costi e soddisfare la domanda

towardsdatascience.com

Conclusioni

In questo articolo, sono stati implementati in Python utilizzando Pyomo due modelli di programmazione intera mista per il problema di dispersione delle strutture. Come verificato in precedenza nella letteratura, il modello di maxisum può portare a una distribuzione irregolare degli elementi nello spazio, mentre il modello di dispersione p ha prodotto soluzioni con posizioni ben distribuite nello spazio e uniformemente distribuite. Un obiettivo di maxisum può essere utilizzato per affinare le soluzioni del modello di dispersione p selezionando quella con la maggiore dispersione totale dall’insieme delle soluzioni ottimali. Il codice completo utilizzato in questi esempi è disponibile per un ulteriore utilizzo.

Riferimenti

Balinski, M. L. 1965. Integer programming: methods, uses, computations. Management science, 12(3), 253–313.

Bynum, M. L., Hackebeil, G. A., Hart, W. E., Laird, C. D., Nicholson, B. L., Siirola, J. D., … & Woodruff, D. L. 2021. Pyomo-optimization modeling in python (Vol. 67). Berlin/Heidelberg, Germany: Springer.

Kuby, M. J. 1987. Programming models for facility dispersion: The p‐dispersion and maxisum dispersion problems. Geographical Analysis, 19(4), 315–329.