Continuiamo la nostra presentazione di PuLP mostrando un esempio completo. A tal fine prendiamo in considerazione un problema reale: quello della schedulazione delle pompe idrauliche degli impianti di sollevamento dei serbatoi degli acquedotti.
Ipotizziamo di voler programmare l'accensione e spegnimento orario di 5 pompe nell'arco di 24 ore. Siccome il comportamento risultante di più pompe che funzionano contemporaneamente non può essere considerato lineare, le 5 pompe vengono modellate come 32 pompe virtuali distinte e mutuamente esclusive.
numPumps = 32 pumpSet = range(numPumps) numTimes = 24 timeSet = range(numTimes)
Ognuna delle 32 pompe virtuali ha un consumo orario di energia e una portata oraria:
energyForPump = [0, 595, 260, 855, 260, 855, 520, 1115, 445, 1040, 705, 1300, 705, 1300, 965, 1560, 595, 1190, 855, 1450, 855, 1450, 1115, 1710, 1040, 1635, 1300, 1895, 1300, 1895, 1560, 2155] dischargeForPump = [0, 1800, 828, 2600, 828, 2600, 1650, 3450, 1440, 3235, 2260, 4060, 2260, 4060, 3090, 4890, 1800, 3600, 2620, 4420, 2620, 4420, 3450, 5250, 3235, 5035, 4060, 5860, 4060, 5860, 4890, 6690]
Tra i vari obiettivi del modello c'è la minimizzazione del costo indotto dall'uso dell'energia elettrica consumata, ipotizzando un costo orario diversificato.
energyCost = [100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 200, 200, 200, 200, 200, 100, 100]
Il serbatoio nel quale stiamo pompando deve soddisfare una domanda oraria d'acqua giornaliera. Nel nostro esempio la domanda oraria è espressa come percentuale rispetto all'intera domanda giornaliera.
waterDemand = [1.92, 1.55, 1.55, 1.55, 1.92, 2.73, 3.91, 5.03, 5.84, 6.22, 6.22, 6.22, 5.47, 5.47, 5.84, 5.84, 5.84, 5.47, 5.03, 4.28, 3.91, 3.10, 2.73, 2.36] totalDemand = 10000 for i in range(numTimes): waterDemand[i] = totalDemand / 100 * waterDemand[i] pass
Ancora sul serbatoio ci sono degli ovvi limiti fisici, quali il volume minimo e massimo. In oltre c'è da specificare anche il volume d'acqua preesistente nel serbatoio al periodo considerato dall'ottimizzazione.
minVolume = 1 * 2600 maxVolume = 7 * 2600 initVolume = 3 * 2600
Finora abbiamo preso in considerazione i dati di input del modello che stiamo costruendo. Passiamo ora alle variabili.
Il primo insieme ovvio di variabili modella lo stato delle singole pompe. Si tratta di variabili binarie, dove il valore 0 denota una pompa spenta, mentre il valore 1 denota una pompa accesa. Queste variabili compongono una matrice, dove il primo indice rappresenta la pompa e il secondo indice l'ora relativa.
pumpStatus = LpVariable.matrix("pumpStatus", (pumpSet, timeSet), 0, 1, LpInteger)Il secondo insieme di variabili sono di tipo continuo e rappresentano il volume d'acqua all'interno del serbatoio nelle singole ore.
volumes = LpVariable.matrix("volumes", (timeSet), minVolume, maxVolume)Il terzo e ultimo insieme di variabili che consideriamo sono delle variabili di scarto che serviranno per modellare la variazione oraria dell'acqua all'interno del serbatoio:
volumesSPlus = LpVariable.matrix("volumesSPlus", (timeSet), 0, maxVolume)
volumesSMinus = LpVariable.matrix("volumesSMinus", (timeSet), 0, maxVolume)Come vediamo, PuLP ci ha offerto degli strumenti molto semplici e al tempo stesso efficaci per modellare delle situazioni tipiche dei modelli di programmazione lineare.
Passiamo ora ad analizzare i vincoli del nostro modello. Come già sappiamo questi si applicano arricchendo un “problema”:
prob = LpProblem("test1", LpMinimize)Per come abbiamo modellato il nostro sistema, una sola pompa alla volta può risultare accesa:
for t in timeSet: prob += lpSum( [pumpStatus[i][t] for i in pumpSet]) == 1 pass
La funzione lpSum fa parte del modulo pulp e serve per sommare tutti i termini della lista che ha come argomento.
Aggiungiamo poi i vincoli di conservazione del flusso d'acqua, che dicono che ad ogni ora il volume del serbatoio è pari al volume all'ora precedente, meno l'acqua usata per soddisfare la domanda, più l'acqua pompata dalla pompe.
for t in range(1,numTimes): prob += volumes[t] == volumes[t-1] - waterDemand[t] + \ lpSum( [pumpStatus[i][t] * dischargeForPump[i] for i in pumpSet]) pass prob += volumes[0] == initVolume - waterDemand[0] + \ lpSum( [pumpStatus[i][0] * dischargeForPump[i] for i in pumpSet])
Aggiungiamo infine i vincoli che modellano le variabili di scarto:
for t in range(1, numTimes): prob += volumes[t-1] - volumes[t] == volumesSPlus[t] - volumesSMinus[t] pass prob += initVolume - volumes[0] == volumesSPlus[0] – volumesSMinus[0]
Il nostro modello si completa con la funzione obiettivo. A titolo esemplificativo ne consideriamo 2. La prima serve a minimizzare il costo dell'energia elettrica consumato per il pompaggio; la seconda serve a minimizzare la variazione del livello d'acqua all'interno della pompa:
costEnergy = \ lpSum ( [ \ lpSum ( [ energyForPump[i]*energyCost[t]*pumpStatus[i][t] \ for i in pumpSet ] ) \ for t in timeSet ] ) levelVariation = lpSum ( [ volumesSPlus[t] + volumesSMinus[t] for t in timeSet ] )
Passiamo ora alla fase risolutiva del nostro modello, cioè invochiamo il nostro solver con in input il nostro modello. Nel caso specifico useremo CBC dell'organizzazione COIN.
myOptions = ["sec 90"] solver = COIN(path = ["C:\\usr\\bin\\cbc.exe", "C:\\usr\\bin\\cbc.exe"], keepFiles = 1, msg = 1, options = myOptions)
Di interessante c'è da notare che è possibile passare degli argomenti specifici del solver attraverso il parametro formale options; nel caso specifico abbiamo fissato un time-out di 90 secondi.
prob += levelVariation status = solver.solve(prob)
Il valore di ritorno dell'invocazione del metodo solve() ci indica lo stato della risoluzione, quindi ammissibile, non ammissibile, non limitato, etc. Notiamo pure come è semplice risolvere lo stesso modello con risolutori diversi.
Per analizzare il valore delle varie variabili possiamo scrivere cose del tipo:
print "Volume", for t in timeSet: print value(volumes[t]), pass print print "Discharge", for t in timeSet: for i in pumpSet: if value(pumpStatus[i][t])==1: print value(dischargeForPump[i]), break pass pass print print "Pumping", for t in timeSet: for i in pumpSet: if value(pumpStatus[i][t])==1: print i, break pass pass print
Mentre il valore della funzione obiettivo si ottiene con:
print "objective=", value(prob.objective)
Spero che l'esempio appena discusso vi abbia convinto della possibilità di usare PuLP all'interno di processi didattici per la programmazione lineare, ad esempio all'interno di laboratori e per organizzare delle prove pratiche.
Unico vero neo di PuLP è che non fornisce un modo per accedere al valore delle variabili duali. Ma, dopo aver analizzato il codice sorgente di PuLP posso dire che è solo un dettaglio tecnico.
I sorgenti completi dell'esempio presentato sono disponibili presso http://code.google.com/p/pump-scheduling/
Commenti