Appunti del corso di Analisi delle Decisioni
Ottimizzazione con Incertezza
Prof. Giampaolo Liuzzi
a.a. 2004/2005
1
Premessa
Queste note costituiscono il materiale didattico per il corso di Analisi delle
Decisioni erogato per gli studenti della Laurea magistrale in Ingegneria Elettrica presso la Facoltà di Ingegneria dell’Università degli Studi di Cassino. Gli
appunti sono in larga parte stati estratti dal testo [1] di Birge e Louveaux e
dalle dispense del corso tenuto dal Prof. Rockafellar [2] presso l’Università di
Washington. Per eventuali approfondimenti si rimanda dunque ai sopracitati
riferimenti [1, 2] ai quali va tributata anche la gran parte degli argomenti coperti da questa deispensa. Tutti gli errori, qualora ve ne fossero, sono invece da
imputarsi esclusivamente a disattenzione dello scrivente che per questo si scusa
anticipatamente.
2
Nozioni preliminari
È noto che un generico problema di ottimizzazione può essere posto nella forma
minx
c.v.
f◦ (x)
x ∈ D,
ove x ∈ ℜn è detto vettore delle variabili di decisione, f : ℜn → ℜ è detta
funzione obiettivo e D ⊆ ℜn è detto insieme ammissibile. A seconda delle proprietà della funzione f◦ (x) e dell’insieme D parleremo di ottimizzazione lineare
o non lineare, ottimizzazione vincolata o non vincolata, ottimizzazione continua
o combinatoria e cosı̀ via.
Per quanto concerne queste note è comodo introdurre la seguente funzione
obiettivo estesa f : ℜn → ℜ̄ tale che
f◦ (x)
se x ∈ D,
f (x) =
+∞
se x 6∈ D,
per cui è possibile riscrivere il generico problema di ottimizzazione come
min f (x).
x
1
Dato uno spazio di probabilità (Ω, F , P ), sia ω : (Ω, F ) → (ℜn , B n ) una v.a.
e supponiamo che la funzione obiettivo f◦ (x, ω) e l’insieme ammissibile D(ω)
dipendano dalle realizzazioni di ω. Anche in questo caso possiamo pensare alla
funzione obiettivo estesa definita come
f◦ (x, ω)
se x ∈ D(ω),
f (x, ω) =
+∞
se x 6∈ D(ω).
Nel caso attuale, in cui i dati del problema dipendono esplicitamente dalla v.a.
ω, non è più immediato riconoscere un problema di ottimizzazione nel senso che
non è ben chiaro cosa dobbiamo ottimizzare e con quali vincoli dato che sia la
funzione obiettivo che la regione ammissibile dipendono da un elemento incerto
ovvero dalle realizzazioni della v.a. ω. In altri termini, non ha certamente senso
considerare il problema
min f (x, ω)
(1)
x
2.1
Formulazione deterministica
Ci sono diversi modi di trattare l’incertezza in un problema di ottimizzazione.
Un primo modo consiste nello stimare un valore della v.a. (come per esempio il
suo valore atteso IEω [ω]), sia esso ω̄, e quindi risolvere il problema
min f (x, ω̄)
x
(2)
In questo caso, la funzione obiettivo f (x, ω̄) non dipende più dalla v.a. ω
che è fissata al valore ω̄ e quindi il problema è nuovamente un problema, per
cosı̀ dire, deterministico che può essere risolto come usualmente si farebbe.
Questo approccio, benchè perfettamente giustificabile in alcuni ambiti applicativi, potrebbe non avere molto senso e, a volte, portare a soluzioni certamente
non ottime come vedremo più avanti con alcuni esempi.
2.2
Osservazione e ottimizzazione
In alcuni casi è possibile posticipare la scelta delle variabili di decisione x ad un
momento successivo alla osservazione della v.a. ω il che corrisponde a quella che
viene comunemente definita strategia wait-and-see. In questo caso, ci troviamo
a dover risolvere un problema di ottimizzazione in x che dipende da ω come
parametro. Potremmo cosı̀ pensare di risolvere il problema (2) solo per il valore
ω̂ corrispondente alla realizzazione della v.a. oppure potremmo concentraci sulla
determinazione del valore ottimo
ψ(ω) = inf f (x, ω)
x
e dell’insieme di soluzioni ottime
Ψ(ω) = arg min f (x, ω),
x
entrambi funzione della variabile aleatoria ω.
2
2.3
Ottimizzazione e osservazione
Contrariamente al caso precedente, potrebbe non essere possibile posticipare la
scelta di x fin dopo aver osservato la realizzazione della v.a. ovvero è necessario adottare quella che comunemente viene definita strategia here-and-now.
In questo assetto le decisione vanno prese subito avendo a disposizione solo la
conoscenza sulla v.a. ω che ci viene dalla sua f.d. In questo caso è necessario
riflettere sul fatto che per ogni scelta delle variabili x non abbiamo un valore
certo f (x) ma, piuttosto, un valore incerto f (x, ω). In altri termini, per ogni
realizzazione della v.a. ω abbiamo una funzione ω 7→ f (x, ω) ovvero, potremmo
pensare ad f (x, ω) come ad una v.a. composta.
Il caso peggiore
In questo assetto possiamo pensare di risolvere il problema minx f˜(x) dove
f˜(x) = sup f (x, ω).
ω
Ovviamente, questo significa concentrarsi esclusivamente sul peggiore risultato
ottenibile avendo fissato le variabili x senza alcun tentativo di distinguere le
realizzazioni della v.a. in base alla sua distribuzione di probabilità.
Programmazione stocastica
Questo approccio consiste nell’interpretare, per ogni x, f (x, ω) a sua volta come
una v.a. che eredita la sua f.d. da quella della v.a. ω. È dunque possibile,
per ogni x, calcolare il valore atteso fˆ(x) = IEω [f (x, ω)] e quindi risolvere il
problema
min fˆ(x).
(3)
x
2.4
Esempio: il problema del venditore di giornali
Questo problema, meglio noto con il nome di newsvendor problem, è probabilmente il più semplice esempio di ottimizzazione in presenza di incertezza
ed è, pertanto, particolarmente adatto per familiarizzare con alucni fondamentali concetti della programmazione in presenza di incertezza. Il problema è il
seguente.
Il ragazzo che vende i giornali al semaforo deve decidere quanti giornali x acquistare alle prime luci dell’alba dall’editore al costo di c euro per giornale.
Durante la mattinata potrà rivendere i giornali al prezzo di s euro l’uno e, a
fine mattinata, può rivendere all’editore i giornali non venduti ricavando r euro
l’uno. Ovviamente vale la seguente relazione
0 ≤ r < c < s.
Il grosso problema del venditore è che, al momento in cui deve scegliere quanti giornali (x) acquistare dall’editore, non conosce con esattezza il numero di
clienti che acquisteranno il giornale da lui al semaforo. In altri termini, la sua
3
domanda giornaliera di quotidiani può essere considerata come una variabile
aleatoria la cui realizzazione accadrà in un momento successivo a quello in cui
deve acquistare i giornali dall’editore.
Indichiamo con D la v.a. che descrive il numero di giornali effettivamente venduti al semaforo durante la mattinata. La tabella che segue riporta i costi, i
ricavi e il profitto del venditore nella mattinata in funzione della variabile di
decisione x e della v.a. D.
Costi
acquisto giornali: cx
Ricavi
sx se x < D
vendite al semaforo:
sD se x ≥ D
0
se x < D
vendite all’editore:
r(x − D) se x ≥ D
profitto: −f (x, D) = −cx + s min{x, D} + r max{0, x − D}
Il problema da risolvere sarebbe dunque il seguente
min f (x, D),
x≥0
e cioè un problema della forma di (1) in cui la v.a. ω = D.
Supponiamo ora che il venditore possa decidere quanti giornali acquistare dall’editore dopo aver avuto notizia di quanti giornali potrà vendere nella mattinata. Siamo nel caso wait-and-see (cfr. sottosezione 2.2) in cui la v.a. D è trattata
come un parametro. In questo caso la f (x, D) è funzione della sola variabile x
e dipende parametricamente da D.
f (x, D)
D
x
Ovviamente, la soluzione ottima risulta essere x = D e il valore ottimo (c − s)D
ovvero, coerentemente con la notazione fin qui adottata,
ψ(D) = (c − s)D,
e
Ψ(D) = {x ∈ ℜ : x = D}.
Ovviamente, questa tecnica risolutiva ha senso solo quando è lecito supporre
che ci sia una sorta di oracolo che predice al venditore il numero di giornali che
potrà vendere al semaforo. In ogni altro caso, questa tecnica risolutiva non ha
molto senso.
Supponiamo ora che D sia una v.a. discreta cioè abbia la seguente f.d. di
probabilità
4
evento
D = 30
D = 40
D = 50
D = 60
D = 100
FD (·)
1/7
2/7
2/7
1/7
1/7
Un altro modo di approcciare il problema consiste, come accennato nella sottosezione 2.1, nel risolvere un problema in cui si fissa il valore della v.a., per
esempio, al suo valore atteso IE[D] = D̄ = 370/7 e quindi risolvere il problema per D = D̄. In questo caso la soluzione è, banalmente, Ψ(D̄) cioè x = D̄,
y1 = D̄ e y2 = 0 con valore ottimo ψ(D̄) = (c − s)D̄.
Come visto nella sottosezione 2.3, un ulteriore modo di procedere è quello di
risolvere il problema
min max f (x, D).
x≥0 30≤D≤100
Considerato il fatto che la funzione f (x, D) è in questo caso non crescente
rispetto alla variabile D,
f (x, D)
x
D
risulta max30≤D≤100 f (x, D) = f (x, D̂), con D̂ = 30. D’altra parte, siccome,
come si verifica facilmente, la funzione f (x, 30)
f (x, 30)
30
x
è decrescente per 0 ≤ x ≤ 30 e crescente per x ≥ 30, il problema minx≥0 f (x, 30)
ha soluzione ottima x̄ = D̂ = 30 con valore ottimo (c − s)x̄ = 30(c − s).
Questa soluzione corrisponde perfettamente all’atteggiamento del venditore che,
cautamente, sceglie di acquistare il maggior numero di giornali nell’ipotesi di
venderne il minor numero possibile. I limiti di questo approccio sono evidenti.
Nel caso, per esempio, in cui sia possibile anche l’evento D = 0 la strategia
che dovremmo adottare sarebbe x = 0 ovvero il venditore non acquista alcun
giornale per limitare i danni nel caso peggiore in cui gli dovesse capitare di non
venderne nessuno.
5
Per concludere la rassegna dei vari approcci al problema rimane da analizzare
quello detto di stochastic programming (programmazione stocastica). Secondo
quanto visto nella sottosezione 2.3, questo tipo di approccio consiste nel considerare, per ogni valore x ammissibile, la f (x, D) come una nuova v.a. per la quale
è possibile calcolare il valore atteso IED [f (x, D)] e quindi risolvere il problema
min IED [f (x, D)].
(4)
x≥0
Nel caso in esame, in cui D è una v.a. discreta, la funzione che otteniamo è
f (x) = IED [f (x, D)] =
1
2
2
1
1
f (x, 30)+ f (x, 40)+ f (x, 50)+ f (x, 60)+ f (x, 100)
7
7
7
7
7
che è una funzione convessa, lineare a tratti e continua ma non differenziabile.
Sia ∂f (x) il subdifferenziale di f in x cosı̀ definito:
∂f (x) = {t ∈ ℜ : f (y) ≥ f (x) + t(y − x) per ogni y ∈ ℜ}.
Per la funzione f (x) risulta

(c − s)


 1

(c − r) +


 73
7 (c − r) +
∂f (x) ∋
5

7 (c − r) +


6

(c − r) +


 7
(c − r)
6
7 (c
4
7 (c
2
7 (c
1
7 (c
− s)
− s)
− s)
− s)
se
se
se
se
se
se
x < 30
30 ≤ x < 40
40 ≤ x < 50
50 ≤ x < 60
60 ≤ x < 100
x ≥ 100
ovvero, più sinteticamente,
∂f (x)
∋ (c − r)P [D ≤ x] + (c − s)P [D > x]
= (c − r)FD (x) + (c − s)(1 − FD (x))
= (c − s) + (s − r)FD (x).
Come è noto, considerata la convessità di f (x), x⋆ è un punto di minimo di f (x)
se e solo se
0 ∈ ∂f (x⋆ ).
Nel caso in esame in cui f (x) è funzione di una sola variabile e ∂f (x) ∋ c−s < 0,
per ogni x < 30 , è possibile determinare il punto di minimo di f (x) cercando il
più piccolo x⋆ ∈ ℜ tale per cui
∂f (x⋆ ) ∋ (c − s) + (s − r)FD (x⋆ ) > 0,
da cui otteniamo
−1
x⋆ = FD
s−c
s−r
−1
avendo indicato con FD
(α) = min{x : FD (x) ≥ α} ovvero l’α-quantile della
distribuzione di probabilità FD (ω). Notiamo che x⋆ corrisponde al punto in
6
cui la pendenza di f (x) passa da negativa (f (x) decrescente) a positiva (f (x)
crescente).
La figura che segue riporta l’andamento della funzione f (x), quando r = 0.8,
c = 1.6 e s = 1.8, da cui è facile riconoscere che il valore ottimo del problema
(4) è
−1
x⋆ = FD
(0.2) = 40.
IED [f (x, D)]
30
40
50
60
100
x
Supponiamo ora che D sia una v.a. continua, cioè una v.a. con f.d. FD (ω)
continua, per cui risulta
Z ∞
f (x) = IE[f (x, D)] =
f (x, ω)dFD (ω).
0
Notiamo che, anche in questo caso, essendo la funzione f (x, D) convessa (e
lineare a tratti rispetto ad x), anche la funzione f (x) è convessa rispetto ad x.
Allora
Z x
Z ∞
f (x) = cx +
[−rx − (s − r)ω]dFD (ω) − sx
dFD (ω),
0
x
da cui, mediante integrazione per parti, otteniamo
Z x
f (x) = (c − s)x + (s − r)
FD (ω)dω.
0
Derivando rispetto a x, otteniamo, come nel caso precedente,
f ′ (x) = (c − s) + (s − r)FD (x).
Percui, ricordando che f (x) è convessa rispetto a x, x⋆ è punto di minimo se e
solo se
f ′ (x⋆ ) = (c − s) + (s − r)FD (x⋆ ) = 0,
da cui, anche in questo caso, otteniamo
FD (x⋆ ) =
s−c
s−r
s−c
s−r
e quindi, nuovamente,
−1
x⋆ = FD
7
.
3
Fondamenti di Programmazione Stocastica
Come gia in parte anticipato nel corso della precedente sezione, le principali
caratteristiche della programmazione stocastica sono le seguenti:
1. here-and-now, vale a dire che, non ostante l’ingente coinvolgimento di
accadimenti futuri, ogni cosa è mirata al miglioramento di decisioni che
devono essere prese nel presente ovvero prima che qualunque v.a. possa
realizzarsi;
2. ricorsione, ovvero l’opportunità che è data al modellista di tenere conto
del fatto che le decisioni prese nel presente possano essere, almeno in parte,
corrette in tempi successivi;
3. informazione e osservazione, le decisioni in tempi successivi possono rispondere a informazioni che si sono rese disponibili da quando sono state prese
le decisioni iniziali. Questa informazioni è modellata mediante l’osservazione di v.a.;
4. convessità, teoria e metodi risolutivi di programmazione stocastica sono
sviluppati sotto l’ipotesi restrittiva di convessità;
5. indipendenza delle misure di probabilità dalle decisioni, si assume cioè che
le v.a. siano indipendenti dalle decisioni.
Uno degli aspetti più caratteristici della programmazione stocastica ed anche
quello che maggiormente la differenzia dalla programmazione matematica è certamente quello dinamico. Più precisamente, in un contesto di programmazione
stocastica è opportuno familiarizzare con il cosı̀ detto processo ricorsivo nel quale
le decisioni si alternano con le osservazioni. In particolare possiamo pensare al
seguente processo ricorsivo
u0 ∈ ℜn0
ω 1 ∈ Ω1
decisione iniziale
osservazione
u1 ∈ ℜn1
ω 2 ∈ Ω2
decisione di ricorsione
osservazione
u2 ∈ ℜn2
..
.
ω N ∈ ΩN
decisione di ricorsione
uN ∈ ℜnN
decisione di ricorsione
osservazione
In questo processo ricorsivo, abbiamo uno stadio iniziale nel quale viene presa la
decisione u0 ∈ ℜn0 e successivamente N stadi di ricorsione ciascuno composto
da una osservazione, nella quale una v.a. si realizza, e da una decisione, che
8
dovrebbe essere la risposta all’osservazione appena avvenuta. Alla fine di questo
processo ricorsivo, avremo un vettore di decisioni
u = (u0 , u1 , . . . , uN ) ∈ ℜn con n = n0 + n1 + . . . + nN
e un vettore di osservazioni che potremmo definire la storia delle osservazioni e
rappresentare come
ω = (ω1 , . . . , ωN ) ∈ Ω con Ω = Ω1 × . . . × ΩN .
Un altro ingrediente essenziale è il costo del processo ricorsivo che può essere
pensato come una funzione f : ℜn × Ω → ℜ̄ con
f (u, ω) = f (u0 , u1 , . . . , uN , ω1 , . . . , ωN ).
Naturalmente la funzione f (u, ω) ha valori nello spazio esteso dei reali ovvero
può assumere anche il valore +∞. Più precisamente, f (u, ω) = +∞ quando
u 6∈ C(ω).
A questo punto, potremmo pensare che lo scopo dell’ottimizzatore sia quello di
risolvere il seguente problema
min IEω [f (u, ω)],
u
ovvero un problema come il problema (3). Ovviamente, questo non è il caso. Infatti se facessimo cosı̀, ci troveremmo, nuovamente, a dover scegliere
u = (u0 , u1 , . . . , un ) prima di conoscere ω = (ω1 , . . . , ωN ), il che distruggerebbe
immediatamente tutti i vantaggi di avere decisioni e osservazioni mescolate le
une con le altre. In particolare, nel nostro caso, se ci poniamo nello stadio
k-esimo possiamo partizionare il vettore ω = (ω1 , . . . , ωN ) nei due sottovettori
(ω1 , . . . , ωk )
(ωk+1 , . . . , ωN )
informazione disponibile
incertezza residua.
Nel fare questo, bisogna sempre tenere presente che il fine ultimo dell’ottimizzazione è quello di produrre una decisione iniziale u0 e che le v.a. ωi , i = 1, . . . , N ,
non saranno mai veramente osservate ma, piuttosto, si potrà al più supporre che
lo siano state. Per esempio, al generico stadio k ci troveremo a dover determinare
uk supponendo che siano state osservate ω1 , . . . , ωk . Questo è il motivo per cui,
allo stadio k, siamo interessati a determinare, non tanto un vettore uk ∈ ℜnk ,
quanto piuttosto una funzione di ricorsione definita sullo spazio Ω1 × . . . × Ωk
uk (·) : (ω1 , . . . , ωk ) 7→ uk (ω1 , . . . , ωk ) ∈ ℜnk .
Nello scegliere una tale funzione uk (·) stiamo, in sostanza, specificando nel presente (here-and-now) come risponderemo ad ogni possibile realizzazione delle
prime k osservazioni.
Definiamo strategia una funzione u : Ω → ℜn ovvero
u(ω) = (u0 , u1 (ω1 ), u2 (ω1 , ω2 ), . . . , uN (ω1 , . . . , ωN )).
9
(5)
Indichiamo con U l’insieme di tutte le possibili funzioni u(·) di questo tipo.
Notiamo che la componente uk (ω) dipende da ω1 , . . . , ωk e non da ωk+1 , . . . , ωN .
Questa proprietà, detta non anticipo delle strategie, serve a garantire che le
decisioni non siano basate sulla conoscenza di eventi futuri prima che questi
accadano. Il problema al quale arriviamo è il seguente
min IEω [f (u(ω), ω)],
u(ω)∈U
dove u(ω) è dato da (5). È importante notare che il precedente è un problema
di ottimizzazione in cui le variabili di decisione sono particolari funzioni distinte
dalla proprietà sopra citata di “non anticipo”.
3.1
Esempio: financial planning and control
Una famiglia americana sta progettando di mandare il proprio figlio al college.
I genitori sanno che tra Y = N v anni, e cioè quando il figlio avrà raggiunto l’età
per potersi iscrivere, la retta per l’intero periodo degli studi sarà di G = 80000
euro. Al momento i genitori dispongono di un budget di b = 55000 euro (b < G)
e devono decidere come investire questi risparmi. Al termine degli Y anni i
genitori potranno
1. prendere in prestito (con un interesse r = 4%) i soldi che servono per
arrivare a coprire la retta di G euro per l’iscrizione al college; oppure
2. depositare in un libretto di risparmio (con un interesse q = 1%) i soldi
avanzati dopo il pagamento della retta di iscrizione.
La famiglia può investire su I tipi diversi di investimento e può cambiare investimento ogni v anni. I genitori hanno, quindi, N = Y /v differenti periodi di
investimento.
Supponiamo, per semplicità di trattazione, che I = 2 (due soli tipi di investimento: stock o bond), N = 3 (tre periodi) e v = 5 (5 anni di durata minima
di ogni investimento). A seconda dell’andamento dei mercati finanziari, si può
avere un interesse di 1.25 per gli stock e 1.14 per i bond oppure 1.06 per gli stock
e 1.12 per i bond con uguali probabilità (in condizione di massima incertezza).
Il processo decisionale è quella gia visto in precedenza in cui abbiamo N = 3
stadi di ricorsione ovvero
u0 = (u10 , u20 )T ∈ ℜ2
ω1 ∈ Ω
decisione iniziale
osservazione
u1 (ω1 ) = (u11 (ω1 ), u21 (ω1 ))T ∈ ℜ2
ω2 ∈ Ω
decisione di ricorsione
osservazione
u2 (ω1 , ω2 ) = (u12 , u22 )T ∈ ℜn2
ω3 ∈ Ω
decisione di ricorsione
osservazione
u3 (ω1 , ω2 , ω3 ) = (u13 (ω1 , ω2 , ω3 ), u23 (ω1 , ω2 , ω3 ))T ∈ ℜ2
decisione di ricorsione
10
dove Ω = {up, down} contiene i due soli andamenti possibili per i mercati
finanziari a cui corrispondono i tassi di interesse degli investimenti disponibili.
Stadio iniziale: Le variabili di decisione sono indipendenti dalle v.a. (non
anticipo) e rappresentano l’ammontare degli investimenti in stock e bond,
rispettivamente, all’inizio degli Y = 15 anni. Esse devono soddisfare il
seguente vincolo
u10 + u20 = b,
ovvero, inizialmente, gli investimenti devono essere esattamente uguali al
budget di b euro disponibili all’inizio del periodo di investimento. Potremmo scrivere questo vincolo, sinteticamente, come
X
ui0 = b.
i∈I
Primo stadio: Le variabili di decisione di primo stadio dipendono dalla
v.a. ω1 ma non da ω2 e ω3 (deve valere la proprietà di non anticipo).
Come visto nel corso della sezione e considerato il fatto che ω1 è una v.a.
discreta, è possibile definire le funzioni u11 (ω1 ) e u21 (ω1 ) in forma tabellare
ovvero definendo le variabili u11 (up), u11 (down), u21 (up) e u21 (down) che
rappresentano le azioni intraprese nei due scenari di primo stadio possibili
(ed equiprobabili), ovvero ω1 = up e ω1 = down. Queste variabili di primo
stadio devono soddisfare i vincoli
1.25u10 + 1.14u20
=
u11 (up) + u21 (up)
1.06u10 + 1.12u20
=
u11 (down) + u21 (down),
cioè, qualunque sia la realizzazione della v.a. ω1 , quello che si investe nel
quinto anno deve essere uguale al capitale disponibile nel quinto anno.
Sintetizzando, possiamo scrivere questi vincoli come
X
X
tω1 ,i ui0 =
ui1 (ω1 ), ∀ ω1 ∈ Ω.
i∈I
i∈I
Secondo stadio: Le variabili di decisione di secondo stadio dipendono dalle
v.a. ω1 e ω2 ma non (sempre per la proprietà di non anticipo) dalla
v.a. ω3 . Anche in questo caso è possibile definire le funzioni u12 (ω1 , ω2 )
e u22 (ω1 , ω2 ) mediante l’introduzione di tante variabili di decisione quanti
sono i possibili scenari di secondo stadio. Dato che i possibili scenari sono
dati da tutte le possibili combinazioni di una realizzazione di ω1 con una
realizzazione di ω2 , dobbiamo introdurre complessivamente 8 variabili e
precisamente:
u12 (up, up), u12 (up, down), u12 (down, up), u12 (down, down),
u22 (up, up), u22 (up, down), u22 (down, up), u22 (down, down).
11
Queste variabili devono soddisfare i vincoli
1.25u11(up) + 1.14u21 (up)
1.06u11(up) + 1.12u21 (up)
1
1.25u1 (down) + 1.14u21 (down)
1.06u11 (down) + 1.12u21 (down)
= u12 (up, up) + u22 (up, up)
= u12 (up, down) + u22 (up, down)
= u12 (down, up) + u22 (down, up)
= u12 (down, down) + u22 (down, down)
cioè, qualunque siano le realizzazione delle v.a. ω1 e ω2 , quello che si investe nel decimo anno deve essere uguale al capitale disponibile nel decimo
anno. Anche qui, possiamo riscrivere i vincoli di secondo stadio nel modo
seguente.
X
X
tω2 ,i ui1 (ω1 ) =
ui2 (ω1 , ω2 ), ∀ ω1 , ω2 ∈ Ω.
i∈I
i∈I
Ultimo stadio: Anche nell’ultimo stadio vale un discorso analogo a quello
fatto nei due stadi precedenti. Infatti, le variabili di decisione dell’ultimo
stadio sono funzione di ω1 , ω2 e ω3 . La variabile u13 (ω1 , ω2 , ω3 ) rappresenta
quanto denaro deve essere preso in prestito con interesse di r% per poter
avere G euro per pagare la rata di iscrizione al college. Al contrario, la
variabile u23 (ω1 , ω2 , ω3 ) rappresenta l’ammontare che può essere versato su
un libretto di risparmio con interesse di q%, dopo aver pagato la retta di G
euro per il college. Anche in questo ultimo caso, le funzioni u13 (ω1 , ω2 , ω3 )
e u23 (ω1 , ω2 , ω3 ) possono essere definite mediante l’introduzione di tante
variabili quante sono le possibili combinazioni delle realizzazioni delle v.a.
ωi , i = 1, 2, 3. Tali variabili devono soddisfare i seguenti vincoli
1
2
1.25u2 (up, up) + 1.14u2 (up, up)
1
2
1.25u2 (up, down) + 1.14u2 (up, down)
1
2
1.25u2 (down, up) + 1.14u2 (down, up)
1
2
1.25u2 (down, down) + 1.14u2 (down, down)
1
2
1.06u2 (up, up) + 1.12u2 (up, up)
1
2
1.06u2 (up, down) + 1.12u2 (up, down)
1
2
1.06u2 (down, up) + 1.12u2 (down, up)
1
2
1.06u2 (down, down) + 1.12u2 (down, down)
=
=
=
=
=
=
=
=
1
2
G − u3 (up, up, up) + u3 (up, up, up)
1
2
G − u3 (up, down, up) + u3 (up, down, up)
1
2
G − u3 (down, up, up) + u3 (down, up, up)
1
2
G − u3 (down, down, up) + u3 (down, down, up)
1
2
G − u3 (up, up, down) + u3 (up, up, down)
1
2
G − u3 (up, down, down) + u3 (up, down, down)
1
2
G − u3 (down, up, down) + u3 (down, up, down)
1
2
G − u3 (down, down, down) + u3 (down, down, down)
cioè, qualunque siano le realizzazione delle v.a. ω1 , ω2 e ω3 , quello di cui
si dispone nel quindicesimo anno deve essere pari all’importo della retta
G, eventualmente prendendo a prestito la quantità u13 oppure depositando
i contanti in avanzo u23 in un libretto di risparmio. Anche questo ultimo
gruppo di vincoli possono essere riscritti sinteticamente come
X
tω3 ,i ui2 (ω1 , ω2 ) + u13 (ω1 , ω2 , ω3 ) − u23 (ω1 , ω2 , ω3 ) = G, ∀ ω1 , ω2 , ω3 ∈ Ω.
i∈I
Supponiamo che le tre v.a. siano indipendenti l’una dalle altre, cosicché avremo
p(ω1 , ω2 , ω3 ) = p(ω1 )p(ω2 )p(ω3 ) = 0.125,
12
e, quindi, per la funzione obiettivo (da massimizzare)
X X X
f (u, ω) =
p(ω1 , ω2 , ω3 )(−ru13 (ω1 , ω2 , ω3 ) + qu23 (ω1 , ω2 , ω3 )).
ω1 ∈Ω ω2 ∈Ω ω3 ∈Ω
Risolvendo il problema otteniamo la seguente soluzione ottima:
u10 = 41479.3,
ω1
up
up
up
up
down
down
down
down
4
ω2
up
up
down
down
up
up
down
down
ω3
up
down
up
down
up
down
up
down
u20 = 13520.7
u11
u21
65094.6
2168.14
36743.2
u12
u22
83839.9
0
0
71428.6
0
71428.6
64000
0
22368
u13
0
0
0
0
0
0
0
12160
u23
24799.9
8870.3
1428.57
0
1428.57
0
0
0
Programmazione Stocastica lineare a due stadi
Consideriamo nuovamente il caso nel venditore ambulante di giornali visto nella sezione precedente. È piuttosto facile ravvisare nel problema del venditore
un problema di programmazione stocastica a due stadi con una sola v.a. (la
domanda D di giornali).
x∈ℜ
D∈Ω
u 1 ∈ ℜ2
primo stadio: decisione sul numero di giornali da acquistare
osservazione: diventa noto il numero di giornali venduti
secondo stadio: decisioni di ricorsione
il vettore delle decisioni di ricorsione è composto di due elementi e precisamente
u11 ≡ giornali venduti al semaforo = min{x, D} e
u12 ≡ giornali rivenduti all’editore = max{0, x − D}.
Riassumendo, il venditore deve decidere “here-and-now” quanti giornali x ∈ ℜ
acquistare dall’editore non sapendo quanto varrà la domanda D nella mattinata.
Tuttavia, quale che sia la realizzazione della v.a. D, egli reagirà all’osservazione
della v.a. con una decisione di secondo stadio o funzione di ricorsione u1 (D) :
Ω → ℜ2 essendo (u1 (D))1 il numero di giornali venduti al semaforo e (u1 (D))2
il numero di giornali rivenduti all’editore.
Più on generale, in un modello stocastico a due stadi si indica con
- x ∈ ℜn1 il vettore delle decisioni di primo stadio;
13
- ω ∈ Ω la v.a. la cui realizzazione fa da margine tra primo e secondo stadio;
- y(ω) ∈ ℜn2 il vettore delle decisioni di secondo stadio.
Il problema è dunque quello di risolvere
min IEω [f (x, y(ω))],
x,y(ω)
potendo f (x, y(ω)) essere una generica funzione (non lineare) e a valori sullo
spazio esteso ℜ̄. Poniamoci ora in un contesto lineare. In questo assetto alla
decisione di primo stadio x corrisponde un costo lineare cT x e dei vincoli lineari
in forma standard Ax = b, x ≥ 0. Allo steso modo, alla decisione di secondo stadio y(ω) corrisponderà un costo lineare q(ω)T y(ω) e dei vincoli lineari
W (ω)y(ω) = h(ω) − T (ω)x, y(ω) ≥ 0. Il problema è dunque
min cT x + IEω [q(ω)T y(ω)]
Ax = b, x ≥ 0
W (ω)y(ω) = h(ω) − T (ω)x q.c.
y(ω) ≥ 0 q.c.
(6)
ove le matrici W (ω) e T (ω) sono dette, rispettivamente, matrice di ricorsione e
matrice della tecnologia. Se la matrice W non dipende dalla v.a. ω si parla di
problema con ricorsione fissa.
Notiamo che i vincoli di secondo stadio dipendono dalla realizzazione della v.a.
ω. Nel seguito supporremo che tali vincoli siano soddisfatti q.c. ovvero quasi
certamente cioè che siano soddisfatti per ogni ω ∈ Ω tranne che al più per ogni
ω ∈ A con Pω (A) = 0 (ovvero A insieme con probabilità nulla).
Se pensiamo ad x e ω come parametri, ovvero ci poniamo nel secondo stadio
quando la decisione iniziale è stata presa e la v.a. si è realizzata, il problema
che ci troviamo ad affrontare è il seguente:
min q(ω)T y
W (ω)y = h(ω) − T (ω)x
y ≥ 0.
Sia ora Q(x, ω) = inf{q(ω)T y : W (ω)y = h(ω) − T (ω)x, y ≥ 0} e Q(x) =
IEω [Q(x, ω)]. Definiamo problema proiettato il seguente
min cT x + Q(x)
Ax = b, x ≥ 0
ove, in pratica, è stata eliminata la dipendenza esplicita dalla v.a. ω. Abbiamo
ricorsione relativamente completa quando Q(x) assume valori finiti per ogni
x ∈ ℜn1 tale che Ax = b e x ≥ 0. Diciamo, invece, che si ha ricorsione completa
quando Q(x) assume valori finiti per ogni x ∈ ℜn1 .
14
5
Esempi di modelli stocastici lineari a due stadi
5.1
Il problema dell’azienda agricola
Una azienda agricola europea è specializzata nella coltivazione di grano, frumento e barbabietole da zucchero e nell’allevamento di mucche da latte. In totale
l’azienda possiede 500 acri di terra che possono essere utilizzati per i diversi tipi
di coltivazione. È noto che ogni anno sono necessari per l’allevamento del bestiame almeno 200 Ton. di farina e 240 Ton. di frumento. Naturalmente l’azienda
può fare fronte a queste necessita o con il proprio raccolto oppure acquistando
da un grossista della zona. La produzione agricola dell’azienda oltre a servire
per l’allevamento del bestiame può essere venduta sul mercato ad un prezzo di
170 euro per una Ton. di farina e 150 euro per una Ton. di frumento. I prezzi
di acquisto di farina e frumento dal grossista sono maggiorati del 40% per via
dei costi di trasporto che quest’ultimo deve sostenere. Per quanto riguarda la
barbabietola, il prezzo di vendita è di 36 euro/Ton; tuttavia, la commissione
europea ha imposto all’azienda una quota di produzione per le barbabietole di
6000 Ton. l’anno. La quantità di barbabietola eventualmente prodotta oltre
questa quota potrà essere venduta ad un prezzo ribassato e precisamente a 10
euro/Ton.
Basandosi sulla propria esperienza passata, l’azienda agricola sa che ogni acro
di terra coltivato a frumento, grano o barbabietola frutta rispettivamente 2.5,
3 e 20 Ton. Per finire, per ogni acro di terra l’azienda deve sostenere dei costi
di semina che sono di 150 euro, 230 euro e 260 euro rispettivamente per grano,
frumento e barbabietole. La tabella che segue riassume i dati fin qui esposti:
Raccolto (Ton./acri)
Costo di semina (euro/acri)
Prezzo di vendita (euro/Ton.)
Grano
2.5
150
170
Frumento
3
230
150
Prezzo di acquisto (euro/Ton.)
Richieste min. (Ton.)
238
200
210
240
Barbabietole
20
260
36 sotto 6000 Ton.
10 sopra 6000 Ton.
–
–
Per scegliere le migliore strategia di semina, l’azienda dovrebbe ricorrere al
seguente modello lineare.
Siano:
x1 = acri di terra piantati a grano;
x2 = acri di terra piantati a frumento;
x3 = acri di terra piantati a barbabietole;
w1 = ton. di grano venduto;
y1 = ton. di grano acquistato;
w2 = ton. di frumento venduto;
15
y2 = ton. di frumento acquistato;
w3 = ton. di barbabietole vendute a prezzo pieno;
w4 = ton. di barbabietole vendute a prezzo ribassato;
Il problema da risolvere è il seguente:
min 150x1 + 230x2 + 260x3 + 238y1 − 170w1 +
+210y2 − 150w2 − 36w3 − 10w4
c.v. x1 + x2 + x3 ≤ 500
2.5x1 + y1 − w1 ≥ 200
3x2 + y2 − w2 ≥ 240
20x3 − w3 − w4 ≥ 0
w3 ≤ 6000
xi , yj , wh ≥ 0 i = 1, 2, 3, j = 1, 2, h = 1, 2, 3, 4.
Risolvendo il problema precedente otteniamo la seguente soluzione ottima
Coltivazione
Grano Frumento
xi
120
80
Raccolto (Ton.)
300
240
wi
100
–
yi
–
–
Profitto complessivo: 118600 euro
Barbabietole
300
6000
6000 (w4 = 0)
L’azienda agricola pur soddisfatta da questa soluzione ottima sa perfettamente
che di anno in anno e a parità di seminato, i raccolti possono variare sensibilmente. in particolare non è usuale avere stagioni con raccolti che sono superiori
o inferiori del 20% rispetto alle stime usate nel problema precedente. Più precisamente, nel caso di una stagione particolarmente buona ogni acro seminato
a grano, frumento e barbabietole frutterà, rispettivamente, 3, 3.6 e 24 Ton. di
raccolto. Vice versa, nel caso di una stagione sotto la media ogni acro seminato
a grano, frumento e barbabietole frutterà, rispettivamente, 2, 2.4 e 16 Ton. di
raccolto. Possiamo a questo punto risolvere due ulteriori problemi di ottimizzazione corrispondenti agli scenari di stagione sopra e sotto la media e ottenere
le seguenti soluzioni ottime:
stagione sopra la media: raccolto + 20%
Coltivazione
Grano Frumento
xi
183.33
66.67
Raccolto (Ton.)
550
240
wi
350
–
yi
–
–
Profitto complessivo: 167667 euro
16
Barbabietole
250
6000
6000 (w4 = 0)
stagione sotto la media: raccolto - 20%
Coltivazione
Grano Frumento
xi
100
25
Raccolto (Ton.)
200
60
wi
–
–
yi
–
180
Profitto complessivo: 59950 euro
Barbabietole
375
6000
6000 (w4 = 0)
Ragionando sul problema in esame è facile convincersi del fatto che la decisione
sulle quantità di terra da seminare con le differenti colture (xi , i = 1, 2, 3) deve
essere presa prima di conoscere l’esito della stagione (se nella norma, sopra la
media o sotto la media). Al contrario le quantità da vendere e comprare dei
differenti prodotti (yi , i = 1, 2 e wj , j = 1, 2, 3, 4) dipendono dal raccolto.
Supponiamo di assegnare a ciascuno dei tre scenari disponibili (sotto la media,
in media e sopra la media) un indice s = 1, 2, 3 e definire le variabili wjs ,
j = 1, 2, 3, 4 e yis , i = 1, 2 dove, per esempio, w32 rappresenta la quantità di
barbabietole vendute a prezzo pieno nel caso di un raccolto nella media.
P
Se ipotiziamo che l’s-esimo scenario abbia probabilita ps = 1/3 con 3s=1 ps = 1
allora possiamo scrivere il seguente problema
min 150x1 + 230x2 + 260x3 +
+ 13 (238y11 − 170w11 + 210y21 − 150w21 − 36w31 − 10w41 )
+ 31 (238y12 − 170w12 + 210y22 − 150w22 − 36w32 − 10w42 )
+ 31 (238y13 − 170w13 + 210y23 − 150w23 − 36w33 − 10w43 )
c.v. x1 + x2 + x3 ≤ 500
2x1 + y11 − w11 ≥ 200
2.4x2 + y21 − w21 ≥ 240
16x3 − w31 − w41 ≥ 0
w31 ≤ 6000
2.5x1 + y12 − w12 ≥ 200
3x2 + y22 − w22 ≥ 240
20x3 − w32 − w42 ≥ 0
w32 ≤ 6000
3x1 + y13 − w13 ≥ 200
3.6x2 + y23 − w23 ≥ 240
24x3 − w33 − w43 ≥ 0
w33 ≤ 6000
xi , yjs , whs ≥ 0 i = 1, 2, 3, j = 1, 2, h = 1, 2, 3, 4, s = 1, 2, 3.
Risolvendo questo problema otteniamo la soluzione seguente:
17
s=1
xi
Raccolto (Ton.)
wi1
Grano
170
340
140
Frumento
80
192
0
s=2
yi1
Raccolto (Ton.)
wi2
0
425
225
48
240
0
s=3
yi2
Raccolto (Ton.)
wi3
0
510
310
0
288
48
yi3
0
Profitto complessivo: 108390 euro
Barbabietole
250
4000
w31 = 4000,
w41 = 0
5000
w32 = 5000,
w42 = 0
6000
w33 = 6000,
w43 = 0
0
Notiamo che, a differenza di quanto accadeva nell’esempio del venditore di giornali, qui la v.a. ω non è direttamente associata ad un valore numerico. La v.a.
ω ha valori sull’insieme Ω = {ω1 , ω2 , ω3 } con ω1 = “stagione sopra la media”,
ω2 = “stagione nella media” e ω3 = “stagione sotto la media”. Nel problema
precedente, gli unici elementi che dipendono dalla v.a. ω sono gli elementi della
matrice della tecnologia essendo


t11 (ω)
0
0
.
0
t22 (ω)
0
T (ω) = 
0
0
t33 (ω)
5.2
Gestione degli investimenti per un impianto di distribuzione elettrica
Un ente per la distribuzione di energia elettrica deve pianificare i propri investimenti in nuovi impianti di generazione per poter soddisfare la domanda
nazionale attuale e futura su un orizzonte temporale di 15 anni. I nuovi impianti
devono essere costruiti all’inizio del primo anno e devono essere operativi per
tutti i 15 anni dell’orizzonte di pianificazione. L’ente dispone di un budget b
di 10 miliardi di euro che possono essere allocati per la costruzione di 4 differenti tipi di impianti di produzione e precisamente: impianto con turbine a
gas, impianto a carbone, impianto nucleare e impianto idroelettrico. Il costo di
costruzione di ciascun tipo di impianto dipende dalla sua capacità ovvero dalla massima potenza elettrica (in GW) erogabile dell’impianto stesso ed infatti,
abbiamo i seguenti costi:
Impianto
Gas
Carbone
Nucleare
Idrico
mln. euro/GW di capacità
110
180
450
950
18
Inoltre, data la limitata disponibilità di corsi d’acqua utilizzabili per produrre
elettricità, un eventuale nuovo impianto idroelettrico non potrà avere una capacità superiore a 5 GW.
Oltre all’investimento iniziale per la costruzione, ciascun tipo di impianto comporta anche dei costi di esercizio (espressi in euro/KWh) come riportato nella
tabella seguente in cui si riporta anche il costo per l’acquisto di un KWh di
energia elettrica da un fornitore estero.
Impianto
Gas
Carbone
Nucleare
Idrico
Acquisto
costo di eser. euro/KWh
ωG
ωC
0.0140
0.0040
0.1500
I costi per la produzione di un KWh di energia elettrica con impianto a gas e a
carbone sono v.a. discrete indipendenti con le seguenti distribuzioni:
ωG
0.0310
0.0330
0.0390
0.0450
0.0490
P (ωG )
0.1
0.2
0.4
0.2
0.1
ωC
0.0170
0.0210
0.0240
0.0290
0.0310
P (ωC )
0.1
0.2
0.4
0.2
0.1
Per pianificare al meglio i propri investimenti, l’ente si basa sulla conoscenza
dell’andamento della richiesta di potenza elettrica per il primo anno. In particolare, sono stati individuati 5 blocchi di consumo e per ciascuno di essi è nota
la durata in ore (nell’anno) e la potenza richiesta (d◦j in GW) come riportato
nella tabella che segue
Blocco j
#1
#2
#3
#4
#5
d◦j (GW)
10.0
8.4
6.7
5.4
4.3
Durata (ore)
490
730
2190
3260
2090
Cosı̀, è noto che, nel corso del primo anno, ci sarà una richiesta di 10 GW
di potenza per un ammontare complessivo di 490 ore. La richiesta sarà di
8.4 GW per 730 ore, 6.7 GW per 2190 e cosı̀ via. Basandosi sui dati in suo
possesso e grazie a delle approfondite ricerche di mercato, l’ente ha stabilito che
nei prossimi 15 anni (periodo di pianificazione dell’investimento), si potranno
verificare i seguenti trend di incremento/decremento delle richieste di potenza
elettrica:
19
ωR
-0.01
0.01
0.03
0.05
0.07
P (ωR )
0.2
0.2
0.2
0.2
0.2
È inoltre noto che le tre v.a. ωG , ωC e ωR sono indipendenti tra di loro. Quindi,
se indichiamo con ΩG , ΩC e ΩR , rispettivamente, lo spazio di tutti i possibili
eventi associati alle tre v.a. ωG , ωC e ωR , avremo che Ω = ΩG × ΩC × ΩR è
lo spazio degli scenari possibili che sono complessivamente pari a |ΩG × ΩC ×
ΩR | = 5 × 5 × 5 = 125. Pertanto, ogni scenario ω = (ωG , ωC , ωR ) ∈ Ω ha una
probabilità data dal prodotto delle probabilità dei singoli eventi cioè: P (ω) =
P (ωG )P (ωC )P (ωR ).
Indichiamo con:
- xi , i = 1, 2, 3, 4, la capacità (in GW) dei nuovi impianti, rispettivamente,
a Gas, a Carbone, Nucleare e Idroelettrico.
- ci , i = 1, 2, 3, 4, i costi di costruzione (in milioni di euro/GW) associati ai
4 differenti tipi di impianto.
- yijk (ω), i = 1, . . . , 5, j = 1, . . . , 5 e k = 1, . . . , 5, la frazione di capacità
elettrica (in GW) utilizzata per produrre elettricità con impianto di tipo
i, per il blocco di richiesta j nell’anno k. Notiamo che, quando i = 5, y5jk
indicherà la quantità di elettricità acquistata da un fornitore estero per
soddisfare la domanda nel blocco j dell’anno k. Queste variabili sono, in
realtà funzioni della v.a. ω.
- q1 (ω), q2 (ω), q3 , q4 , q5 , i costi di esercizio (in euro/KWh), rispettivamente,
di impianti a Gas, Carbone, Nucleare, Idroelettrico e per l’acquisto da
un fornitore estero. Notiamo che i primi due costi dipendono dalle realizzazioni della v.a. ω ∈ Ω.
- hj , j = 1, 2, 3, 4, 5, la durata (in ore) di ciascun blocco di consumo.
- Djk (ω), j = 1, . . . , 5 e k = 1, . . . , 15, la potenza richiesta (in GW) nel
blocco j dell’anno k. In particolare, vale la relazione Djk (ω) = d◦j (1 + (k −
1)ωR ), per ogni j = 1, . . . , 5 e k = 1, . . . , 15, dove d◦j indica la richiesta nel
primo anno per i vari blocchi di domanda. Le informazioni sulle quantità
Djk (ω) possono essere cosı̀ riassunte:
20
j
1
1
1
1
1
..
.
5
5
5
5
5
ωR
-0.1
0.1
0.3
0.5
0.7
..
.
-0.1
0.1
0.3
0.5
0.7
Dj1 (ω)
10
10
10
10
10
..
.
4.3
4.3
4.3
4.3
4.3
Dj2 (ω)
9.9
10.1
10.3
10.5
10.7
..
.
4.257
4.343
4.429
4.515
4.601
...
...
...
...
...
...
···
...
...
...
...
...
Dj14 (ω)
8.7
11.3
13.9
16.5
19.1
..
.
3.741
4.859
5.977
7.095
8.213
Dj15 (ω)
8.6
11.4
14.2
17.0
19.8
..
.
3.698
4.902
6.106
7.31
8.514
Riassumendo, otteniamo il seguente problema di programmazione stocastica
lineare a due stadi.


5 X
5 X
15
X
min cT x + IEω 
qi (ω)hj yijk (ω)
x,y
i=1 j=1 k=1
cT x ≤ 10000
x4 ≤ 5.0
x≥0
yijk (ω) ≤ xi ,
5
X
yijk (ω) ≥ Djk (ω),
j = 1, . . . , 5, k = 1, . . . , 15, i = 1, . . . , 4,
j = 1, . . . , 5, k = 1, . . . , 15
i=1
yijk (ω) ≥ 0,
j = 1, . . . , 5, k = 1, . . . , 15, i = 1, . . . , 5.
Risolvendo il problema, otteniamo il seguente piano di investimenti, in termini
di capacità installate (variabili xi , i = 1, . . . , 4):
Impianto
Gas
Carbone
Nucleare
Idrico
i
1
2
3
4
Capacità xi
3.366 GW
3.953 GW
4.313 GW
5.000 GW
Il costo atteso complessivo del piano di investimento all’ottimo è pari a 15828.81856
milioni di euro. Il piano di investimenti ottimo prevede, tra l’altro, che, date le
installazioni ottime di capacità x⋆ come riportate nella tabella che precede, un
trend di crescita della domanda del 7% e i costi operativi delle centrali a gas e
carbone rispettivamente di 3.9 e 2.4 eurocent per KWh, le potenze erogate dai
singoli impianti nel 15 anno siano pari a (yij15 (ω) per i = 1, . . . , 5, j = 1, . . . , 5
e dove ω = (ωG , ωC , ωR ) con ωG = 0.039, ωC = 0.024 e ωR = 0.07).
21
Impianto
1
3.366
3.953
4.313
5
3.868
Gas
Carbone
Nucleare
Idrico
Fornitore
6
Blocco di domanda
2
3
4
3.366 0.469
0
3.953 3.953 1.757
4.313 4.313 4.313
5
5
5
0.588
0
0
5
0
0
3.815
5
0
Programmazione stocastica lineare a due stadi
con ricorsione fissa
Consideriamo nuovamente il problema (6) ove W (ω) = W ovvero
min
cT x + IEω [q(ω)T y(ω)]
Ax = b, x ≥ 0
W y(ω) = h(ω) − T (ω)x q.c.
y(ω) ≥ 0 q.c.
oppure, equivalentemente, il problema proiettato
min cT x + IEω [Q(x, ω)]
Ax = b, x ≥ 0
dove
Q(x, ω) = inf
qT y
Wy = h − Tx
y ≥ 0.
(7)
(8)
Per ogni realizzazione di ω, Il problema (8) è un problema lineare per il quale
possiamo scrivere il duale che è
max
π T (h − T x)
W T π ≤ q.
(9)
Consideriamo ora la funzione sq (χ) = inf{q T y : W y = χ, y ≥ 0}. Ovviamente
risulta Q(x, ω) = sq (h − T x). Indichiamo con Π(q) l’insieme ammissibile del
problema (9) ovvero
Π(q) = {π : W T π ≤ q}
e con Y (χ) l’insieme ammissibile del problema (8) cioè
Y (χ) = {y : W y = χ, y ≥ 0}.
Dalla teoria della dualità nella programmazione lineare abbiamo che, se Π(q) 6=
∅,
sq (χ) = sup π T χ,
π∈Π(q)
ovvero, secondo la definizione, sq (χ) è una funzione di supporto per l’insieme
Π(q).
22
Sia Q(x) = IEω [Q(x, ω)] la cosı̀ detta funzione di ricorsione mediante la quale
possiamo scrivere, come gia visto, il problema proiettato seguente
min
cT x + Q(x)
Ax = b, x ≥ 0.
(10)
Dal momento che, a seconda delle realizzazioni della v.a. ω, il problema di secondo stadio può essere illimitato inferiormente (Q(x, ω) = −∞), inammissibile
(Q(x, ω) = +∞) oppure avere un valore ottimo finito, il valore atteso mediante
il quale è definita la funzione proiettata Q(x) potrebbe non essere ben definito
se Pω ([Q(x, ω) = −∞] ∪ [Q(x, ω) = −∞]) > 0. Per questo motivo, assumiamo
convenzionalmente che +∞ + (−∞) = +∞. Nel caso in cui ω sia una v.a. discreta, questo equivale a ritenere inammissibile una decisione di primo di stadio
che potrebbe definire una azione di secondo stadio indefinita (+∞ + (−∞)).
Definiamo K1 = {x : Ax = b, x ≥ 0} l’insieme definito dai vincoli di primo
stadio e K2 = {x : Q(x) < +∞} l’insieme ammissibile del secondo stadio
implicitamente definito dalla funzione Q(x). Con questa notazione, è possibile
riscrivere il problema (10) come
min
cT x + Q(x)
x ∈ K1 ∩ K2 .
Più in particolare, nel caso in cui ω sia una v.a. discreta, definiamo gli insiemi
ammissibili elementari come
K2 (ω) = {x : Q(x, ω) < +∞} = Y (h(ω) − T (ω)x).
Sia inoltre
K2P
= {x : ∀ω ∈ Ω Y (h(ω) − T (ω)x) 6= ∅}
= ∩ω∈Ω K2 (ω)
l’insieme di possibilità degli insiemi ammissibili di secondo stadio.
Proposizione 1
a) Per ogni ω, K2 (ω) è un poliedro e K2P è un insieme chiuso e convesso.
b) Se ω è una v.a. discreta, allora K2P è un poliedro e coincide con K2 .
Dim.: Per provare il punto a) è sufficiente notare che l’insieme K2 (ω) è definito
da un insieme di uguaglianze e disuguaglianze lineari. Per quanto riguarda il
punto b) notiamo anzitutto che K2P essendo, in questo caso, l’intersezione di un
numero finito di poliedri, è a sua volta un poliedro. Proviamo ora che K2P = K2 .
Dimostriamo prima di tutto che K2 ⊆ K2P . Sia x ∈ K2 tale che Q(x) < +∞.
Per definizione sappiamo che Q(x) è una somma pesata con coefficienti positivi
di un numero finito di termini Q(x, ω). Dunque, siccome Q(x) < +∞ e grazie
alla convenzione +∞ + (−∞) = +∞, deve risultare Q(x, ω) < +∞ per ogni
23
ω ∈ Ω il che implica x ∈ K2 (ω) per ogni ω ovvero x ∈ K2P . Proviamo ora che
K2P ⊆ K2 . Sia, pertanto, x ∈ K2P . Per tale valore di x risulta Q(x, ω) < +∞
per ogni ω che, nuovamente, implica Q(x) < +∞ e dunque x ∈ K2 .
2
Dalla precedente proposizione segue che, quando ω è una v.a. discreta, l’insieme
ammissibile del problema (10) K1 ∩ K2 è un poliedro. Vediamo ora cosa è
possibile dire riguardo le proprietà della funzione di ricorsione Q(x).
Proposizione 2 Per il problema di programmazione stocastica (7), se Q(x, ω) >
−∞ e ω è una v.a. discreta, allora Q(·, ω) è una funzione convessa e lineare a
tratti sull’insieme K1 ∩ K2 per ogni ω ∈ Ω.
Dim.: Dimostriamo anzitutto che Q(x, ω) è una funzione convessa. Si considerino due vettori x1 , x2 ∈ K1 ∩ K2 e siano Q(x1 , ω) = min{q T y : W y =
h − T x1 , y ≥ 0} = q T y1⋆ e Q(x2 , ω) = min{q T y : W y = h − T x2 , y ≥ 0} = q T y2⋆ .
Sia ora x = λx1 + (1 − λ)x2 , λ ∈ (0, 1), un vettore combinazione convessa di
x1 e x2 . Per il fatto che K1 ∩ K2 è un poliedro e quindi un insieme convesso,
risulta certamente x ∈ K1 ∩ K2 . Sappiamo inoltre che Q(x, ω) = min{q T y :
W y = h − T x, y ≥ 0} = q T y ⋆ . Dal momento che
W y1⋆ = h − T x1 , y1⋆ ≥ 0,
W y2⋆ = h − T x2 , y2⋆ ≥ 0,
risulta che, per ogni λ ∈ (0, 1)
W (λy1⋆ + (1 − λ)y2⋆ ) = h − T x,
λy1⋆ + (1 − λ)y2⋆ ≥ 0,
cioè, il vettore λy1⋆ + (1 − λ)y2⋆ è ammissibile per il problema min{q T y : W y =
h − T x, y ≥ 0}. Quindi, si ha
Q(x, ω) = q T y ⋆
≤ q T (λy1⋆ + (1 − λ)y2⋆ )
= λq T Q(x1 , ω) + (1 − λ)Q(x2 , ω).
Per poter dimostrare che Q(x, ω) è anche una funzione lineare a tratti, ricordiamo anzitutto che dalla definizione dell’insieme K2 e per il fatto che ω è una
v.a. discreta, segue che, per ogni x ∈ K1 ∩ K2 , Q(x, ω) < +∞. Ovvero, per
ogni valore di ω e x ∈ K1 ∩ K2 , il problema min{q T y : W y = h − T x, y ≥ 0} è
sempre ammissibile. Per ipotesi sappiamo che Q(x, ω) > −∞, il che vuol dire
che il problema di secondo stadio non è illimitato inferiormente. Quindi, per
il teorema fondamentale della programmazione lineare, il problema di secondo stadio ammette soluzioni ottime e, almeno una di queste, è un vertice del
poliedro {y : W y = h − T x, y ≥ 0}. Dal momento che un poliedro ammette
un numero finito di vertici e quindi un numero finito di basi ottime, segue che
Q(x, ω) è lineare a tratti. Infatti, sia y ⋆ la soluzione ottima del problema di
secondo stadio e sia B la base associata a detta soluzione ovvero
−1
B (h − T x)
⋆
y =
≥ 0.
0
24
T −1
Abbiamo che Q(x, ω) = q T y ⋆ = qB
B (h − T x), per ogni x tale che B è una
base ottima del problema di secondo stadio da cui segue che Q(x, ω) è lineare a
tratti.
2
A questo punto abbiamo dimostrato che la funzione Q(x, ω), sotto opportune
ipotesi, è convessa e lineare a tratti. Essa ammette dunque un subdifferenziale
ovvero, secondo la definizione, esiste non vuoto l’insieme
∂Q(x, ω) = {ξ ∈ ℜn : Q(y, ω) ≥ Q(x, ω) + ξ T (y − x), ∀ y ∈ K1 ∩ K2 },
per ogni x ∈ K1 ∩ K2 .
Dalla teoria della dualità, sappiamo che se per x ∈ K1 ∩ K2 e ω ∈ Ω il problema
di secondo stadio ammette soluzione ottima allora ammette soluzione ottima
anche il problema (9) e, per il teorema della dualità forte, sappiamo che
Q(x, ω) = (π ⋆ )T (h − T x).
Dal teorema della dualità debole risulta invece
Q(y, ω) ≥ (π ⋆ )T (h − T y),
da cui otteniamo
Q(y, ω) − Q(x, ω) ≥ −(π ⋆ )T T (y − x).
È quindi evidente, in base alla definizione, che il vettore −T T π ⋆ è un subgradiente della funzione Q(x, ω) ovvero un elemento del subdifferenziale ∂Q(x, ω).
Più precisamente, è possibile dimostrare la seguente proposizione
Proposizione 3 Supponiamo che, per un certo x = x0 , il valore di Q(x0 , ω)
sia finito. Allora Q(x, ω) ammette un subdifferenziabile in x0 e risulta
∂Q(x0 , ω) = conv{−T T D(x0 , ω)},
dove D(x, ω) = arg maxπ∈Π(q) π T (h − T x) è l’insieme delle soluzioni ottime del
problema duale (9).
Le proprietà fin qui viste per la funzione Q(x, ω), nel caso in cui ω sia una v.a.
discreta, si estendono senza troppiPproblemi alla funzione di ricorsione. Infatti,
N
dato che Q(x) = IEω [Q(x, ω)] = k=1 pk Q(x, ωk ) ovvero una sommatoria con
coefficienti positivi delle funzioni Q(x, ωk ), k = 1, . . . , N , risulta:
1) Q(x) è una funzione convessa e lineare a tratti per ogni x ∈ K1 ∩ K2 ;
2) Q(x) ammette un subdifferenziale per ogni x ∈ K1 ∩ K2 e risulta
∂Q(x) =
N
X
pk ∂Q(x, ωk ).
k=1
In particolare, se ξk ∈ ∂Q(x, ωk ), k = 1, . . . , N , allora ξ =
∂Q(x).
25
PN
k=1
pk ξk ∈
6.1
Esempio: il problema del venditore di giornali
Consideriamo nuovamente l’esempio visto in sezione 2.4. Come gia abbiamo
visto, il problema del venditore di giornali può essere formulato come problema di programmazione stocastica lineare a due stadi con risorsione fissa. La
situazione del venditore è infatti la seguente
x∈ℜ
D∈Ω
primo stadio: decisione sul numero di giornali da acquistare
osservazione: diventa noto il numero di giornali venduti
y ∈ ℜ2
secondo stadio: decisioni di ricorsione
In particolare, una volta acquistati gli x giornali dall’editore e avvenuta la realizzazione della v.a. D, il venditore adotterà le decisioni di secondo stadio ovvero
deciderà quanti giornali vendere al semaforo e quanti rivenderne all’editore a
fine mattinata. Il problema è quindi
minx
cT x + Q(x)
x ≥ 0,
con Q(x) = IED [Q(x, D)] e
Q(x, D)
= inf −sy1 − ry2
y
y1 + y2 = x
y1 ≤ D
y1 , y2 ≥ 0.
Il problema duale del precedente problema (primale) di secondo stadio assume
la forma seguente
sup xπ1 − Dπ2
π
π1 − π2 ≤ −s
π1 ≤ −r
π2 ≥ 0,
la cui regione ammissibile, non vuota, è riportata in figura. Notiamo che la
regione ammissibile del problema duale, ovvero l’insieme Π(q), non dipende in
alcun modo ne dal valore delle variabili di decisione di primo stadio (x) ne dal
valore della v.a. D. Dal momento che la regione ammissibile del problema
primale, vale a dire l’insieme Y (h − T x), è sempre non vuota e anzi il problema
primale ammette sempre soluzione ottima (qualunque sia il valore di x e D)
pari a
 x


se x ≤ D

⋆ 
D−x

y1
=
y2⋆
 
D


se x > D

x−D
26
@
@@
@@
π2
@@@
6
@@
@
@@@
@@
@
@@@
Regione ammissibile
@@
@
@@
@
@@@
del Duale
@@@
@@
@
p s
@@
@
pp
@@@
p
@
@ @@ p p p
@
@ @ pp
@
@
@
@
@
@rpp p p p p p p p p p
@
@
@
@
@
pp
@
@ @
@ @
@
s−r
@@@@@@@@@@@@@@@@@@@@@@@
pp
@
@
@
@
pp
@@@@@@@@@@@@@@@@@@@@@@@ @
@
pp
@r
@@@@@@@@@@@@@@@@@@@@@@@
p
π1
−s
−r
segue che anche il problema duale ammette sempre soluzione ottima e che almeno una di queste coinciderà con uno dei due vertici della regione ammissibile.
Il valore della soluzione ottima del duale sarà pertanto pari a

⋆ π1
−s


se
=
 −sx
⋆
π2⋆ 0
π1
−r


=
 −rx − D(s − r) se
π2⋆
s−r
Ricordando che, nell’esempio in esame T = (−1, 0)T , abbiamo che il subdifferenziale ∂Q(x, D) contiene il subgradiente −s se −sx = −sy1⋆ − ry2⋆ e −r
se −rx − D(s − r) = −sy1⋆ − ry2⋆ o entrambi se −rx − D(s − r) = −sx.
Percui, se vogliamo stabilire come è fatto il subdifferenziale ∂Q(x, D) quando x = D = 30 e s = 1.8, c = 1.6 e r = 0.8, allora dobbiamo confrontare il
valore Q(30, 30) = −30s con i valori all’ottimo della funzione obiettivo duale
ovvero −30s e −30r − 30(s − r) = −30s. Pertanto
∂Q(30, 30) = conv{−s, −r},
ovvero il subdifferenziale contiene i subgradienti −s e −r e tutte le loro combinazioni convesse.
Ovviamente, potranno esistere valori di x e D per cui la funzione Q(x, D) è
differenziabile ovvero per cui il subdifferenziale si riduce ad un singleton contenente come suo unico elemento il gradiente ∇Q(x, D). Per esempio, se consideriamo per D = 30 il punto x = 40, abbiamo che Q(40, 30) = −30s − 10r =
−40r − 30(s − r), pertanto,
∂Q(40, 30) = {∇Q(40, 30)} = {−r}.
27
7
L-shaped method
Nel caso in cui ω sia una v.a. discreta è possibile, aumentando opportunamente
il numero delle variabili, riscrivere il problema (7) come un unico problema di
programmazione lineare. Supponiamo infatti di associare ad ogni realizzazione
ωk , k = 1, . . . , N , della v.a. le quantità qk , hk , Tk e le variabili di secondo stadio
yk . Possiamo riscrivere il problema (7) come segue.
min
x,y1 ,...,yN
T
c x+
N
X
pk qkT yk
k=1
Ax = b
Tk x + W yk = hk , k = 1, . . . , N
x ≥ 0, yk ≥ 0, k = 1, . . . , N,
(11)
la cui struttura a blocchi, come emerge dalla figura, richiama quella della lettera
“L”, da cui il nome L-shaped per il metodo usato per risolvere problemi di
programmazione stocastica lineare a due stadi con ricorsione fissa.
A
T1
T2
p
p
p
p
p
W
W
p p
p p
p p
p
TN
W
Figura 1: Struttura a blocchi del problema (11).
L’inconveniente del problema (11) è che il numero di variabili e/o vincoli cresce
molto rapidamente con il numero di scenari ovvero con il numero di realizzazioni
possibili per la v.a. ω. Per questo motivo è preferibile affrontare questi problemi sfruttandone maggiormente la struttura e non trattandoli come dei comuni
problemi di programmazione lineare.
Per questa ragione, è preferibile considerare direttamente il problema (7). Per
questo problema, notiamo, anzitutto, che non si tratta di un problema di programmazione lineare per il semplice fatto che la funzione obiettivo cT x + Q(x)
non è una funzione lineare ma una funzione convessa lineare a tratti. Tuttavia,
28
è possibile riscrivere il problema (7) nella seguente forma equivalente
min
cT x + θ
Q(x) ≤ θ
x ∈ K1 ∩ K2 ,
(12)
dove, mediante l’aggiunta di una variabile (θ), ci si è ricondotti ad un problema
di programmazione lineare. Il vincolo Q(x) ≤ θ, infatti, può essere riscritto
come (x, θ) ∈ epi Q(x) ove, lo ricordiamo, Q(x) è una funzione poliedrale per
la quale risulta epi Q(x) essere un poliedro.
Come rappresentare il poliedro K2 .
Notiamo che imporre il vincolo x ∈ K2 è equivalente al richiedere che, per ogni
k = 1, . . . , N ,
hk − Tk x ∈ pos W,
ovvero che i vettori hk − Tk x appartengano all’involucro affine di W . Se introduciamo, per k = 1, . . . , N i problemi
min wk = eT v + eT v −
W y + v + − v − = hk − Tk x,
y, v + , v − ≥ 0,
(Pfk )
ove eT = (1, . . . , 1), abbiamo che hk − Tk x ∈ pos W se e solo se wk⋆ = 0. I
problemi duali dei precedenti problemi (Pfk ), per k = 1, . . . , N , sono i seguenti
max
σ T (hk − Tk x)
σ T W ≤ 0,
−e ≤ σ ≤ e.
(Dfk )
Indicando con σi , i = 1, . . . , I, i vertici della regione ammissibile del problema
duale (Dfk ), abbiamo che condizione necessaria e sufficiente affinché wk⋆ = 0, è
che
σiT (hk − Tk x) ≤ 0, ∀ i = 1, . . . , I.
Pertanto, condizione necessaria e sufficiente affinché x ∈ K2 è che
σiT (hk − Tk x) ≤ 0, ∀ i = 1, . . . , I, k = 1, . . . , N.
(13)
Come rappresentare il poliedro epi Q(x).
Sia xν ∈ K1 ∩ K2 e ricordiamo che, per definizione, Q(xν ) = IEω [Q(xν , ω)] ove
Q(xν , ωk ) = inf
qkT y
W y = hk − Tk xν
y ≥ 0.
29
(Pqk )
Il duale del precedente problema (Pqk ), per k = 1, . . . , N , può essere scritto nella
forma
max π T (hk − Tk xν )
(Dqk )
π T W ≤ qkT .
Siccome xν ∈ K1 ∩ K2 , Q(xν , ωk ) < +∞, per ogni k = 1, . . . , N . Supponiamo
che Q(xν , ωk ) > −∞ allora risulta
Q(xν , ωk ) = max {πjT (hk − Tk xν )},
j=1,...,J
dove πj , j = 1, . . . , J, indicano tutti i vertici della regione ammissibile del
problema duale (Dqk ). Per la funzione di ricorsione Q(x) risulta
Q(xν ) =
N
X
pk max {πjT (hk − Tk xν )}.
k=1
(xν )
Inoltre, indicato con πjk
j=1,...,J
il vettore tale che
(xν )
(πjk )T (hk − Tk xν ) = max {πjT (hk − Tk xν )},
j=1,...,J
ν
(x )
risulta η = −(πjk )T Tk ∈ ∂Q(x̄) un subgradiente di Q(x) in x = xν e, ricordando la convessità di Q(x), per ogni x ∈ K1 ∩ K2 , vale la disuguaglianza del
subgradiente
Q(xν ) + η T (x − xν ) =
N
X
(xν )
pk (πjk )T (hk − Tk x) ≤ Q(x).
k=1
Il vincolo Q(xν ) + η T (x − xν ) ≤ θ definisce una faccia del poliedro epi Q(x)
perchè soddisfa le seguenti proprietà:
1. il punto (xν , θ̄) = (xν , Q(xν )) giace sull’iperpiano Q(xν ) + η T (x − xν ) = θ;
2. per ogni punto (x, θ) ∈ epi Q(x) risulta Q(xν ) + η T (x − xν ) ≤ θ;
3. per ogni punto (x, θ) tale che Q(xν ) + η T (x − xν ) > θ risulta (x, θ) 6∈
epi Q(x).
Dal momento che i vertici del poliedro K1 ∩ K2 nonché i vertici del poliedro del
problema duale sono in numero finito cosı̀ pure sarà il numero dei vincoli
Q(xν ) + η T (x − xν ) ≤ θ,
(14)
mediante i quali si può imporre l’appartenenza del generico punto (x, θ) all’epigrafo della funzione Q(x).
Il metodo L-shaped risolve il problema (7) posto nella forma (12) ma con una
strategia di generazione dei vincoli di tipo (13) o (14). Questo vuol dire che
30
inizialmente si tenta di risolvere il problema senza alcuno dei vincoli (13) o (14).
Se la soluzione di questo problema rilassato (in cui mancano cioè alcuni vincoli)
non viola alcuno dei vincoli mancanti allora questa viene dichiarata ottima ed il
metodo si arresta. Altrimenti, si determina un vincolo che risulta essere violato e
tale vincolo viene aggiunto al problema rilassato che sarà risolto nella successiva
iterazione. Il processo ha convergenza finita perchè al più un numero finito di
vincoli possono essere aggiunti al problema.
L-shaped method
Step 0. Poni r = s = ν = 0.
Step 1. Poni ν = ν + 1.
Se s > 0 risolvi il problema
Se s = 0 risolvi il problema
min cT x + θ
Ax = b
Dl x ≥ dl l = 1, . . . , r
El x + θ ≥ el l = 1, . . . , s
x ≥ 0, θ ∈ ℜ.
min
cT x
Ax = b
Dl x ≥ dl
x ≥ 0.
l = 1, . . . , r
e poni θν = −∞
Sia (xν , θν ) una sua soluzione ottima.
Step 2. If xν 6∈ K2 then
aggiungi un Feasibility cut
poni r = r + 1 e torna allo Step 1.
Endif
Step 3. Do k = 1, . . . , N
risolvi il problema
min wk = qkT y
W y = hk − Tk xν
y ≥ 0.
If wk⋆ = −∞ then STOP Endif.
End Do
definisci un Optimality cut θ ≥ es+1 − Es+1 x.
If θν ≥ es+1 − Es+1 xν then STOP.
else
poni s = s + 1 e torna allo Step 1.
Endif
31
(15)
Feasibility cut
Do k = 1, . . . , N
risolvi il problema
min
wk = eT v + eT v −
W y + v + − v − = hk − Tk xν ,
y, v + , v − ≥ 0,
dove eT = (1, . . . , 1).
Se wk⋆ > 0, sia σkν il vettore delle variabili duali all’ottimo, poni
Dr+1 = (σkν )T Tk e dr+1 = (σkν )T hk ,
e RETURN.
End Do
Optimality cut
Sia πkν il vettore delle variabili duali all’ottimo del problema (15). Poni
Es+1 =
N
X
pk (πkν )T Tk e es+1 =
k=1
N
X
pk (πkν )T hk ,
k=1
e RETURN
Ricordiamo che si definisce involucro affine di una matrice W l’insieme
posW = {χ : W y = χ, y ≥ 0},
ovvero l’insieme delle combinazioni con coefficienti positivi o nulli delle colonne
della matrice W . In base a questa definizione, risulta Y (χ) = posW . Sappiamo
inoltre che dato il poliedro Π(q) = {π : W T π ≤ q}, il suo cono di recessione è il
seguente insieme Π0 = Π(0) = {π : W T π ≤ 0}.
Nella prima iterazione ν = 1 del metodo L-shaped, si risolve come prima cosa
il problema
min cT x
Ax = b, x ≥ 0
determinandone la soluzione ottima xν e si pone θν = −∞. Fatto questo, nello
Step 2 si controlla se il punto xν appartiene all’insieme K2 = {x : Q(x) < +∞}.
Per verificare se xν ∈ K2 , nella procedura Feasibility cut, si risolve, per k =
32
1, . . . , N , il problema
min wk = eT v + + eT v −
W y + v + − v − = hk − Tk xν ,
y, v + , v − ≥ 0,
che ammette sempre soluzione ottima con wk⋆ ≥ 0. Se wk⋆ > 0, risulta Y (hk −
Tk xν ) = ∅ ovvero hk − Tk xν 6∈ posW e quindi xν 6∈ K2 e deve esistere un
iperpiano che separa hk − Tk xν da posW cioè un iperpiano σ T χ = 0 tale che
σT χ > 0
σT χ ≤ 0
per χ = hk − Tk xν
per ogni χ ∈ pos W.
Per determinare un simile iperpipano risolviamo il problema
max
σ T (hk − Tk x)
σ T W ≤ 0,
−e ≤ σ ≤ e.
duale del precedente problema usato per il controllo dell’ammissibilità. All’ottimo avremo determinato un vettore σ̄ tale che
σ̄ T W ≤
σ̄ T (hk − Tk xν ) >
0
0.
(16)
(17)
Dimostriamo che σ̄ T χ = σ̄ T (hk − Tk x) = 0 è l’iperpiano che cerchiamo. Infatti:
1. dalla relazione (17) segue che σ̄ T (hk − Tk xν ) > 0;
2. al contrario, per ogni χ ∈ posW ovvero tale per cui esiste y ≥ 0 e χ = W y,
risulta σ̄ T χ = σ̄ T W y e dalla relazione (16) otteniamo che σ̄ T W y ≤ 0.
Se si arriva ad eseguire lo Step 3, siamo certi che il punto xν appartiene all’insieme K1 ∩ K2 . Nello Step 3, per ogni k = 1, . . . , N , si risolve il problema (15)
determinandone la soluzione ottima wk⋆ . Ovviamente, se il k-esimo problema
(15) è illimitato inferiormente il metodo termina con l’indicazione che il problema complessivo è illimitato inferiormente nel senso che esiste una scelta delle
variabili x = xν tale che, con probabilità pk il costo cT xν + Q(xν ) = −∞.
Viceversa, se per ogni k il problema (15) ammette soluzione ottima (ricordiamo
che non può essere inammissibile), indicate con πkν le variabili duali all’ottimo, possiamo definire, per ogni k = 1, . . . , N , un subgradiente della funzione
Q(x, ωk ) in xν , percui, vista la convessità di Q(x, ω), vale la seguente relazione
Q(x, ωk ) ≥ (πkν )T hk − (πkν )T Tk x,
∀ x ∈ K1 ∩ K2 .
Sommando, per k = 1, . . . , N , le precedenti relazioni, otteniamo
Q(x) ≥
N
X
k=1
pk (πkν )T hk −
N
X
pk (πkν )T Tk x,
k=1
33
∀ x ∈ K1 ∩ K2
Proposizione 4 Sia
min cT x + θ
Ax = b
Dl x ≥ dl l = 1, . . . , r
PN
(xl ) T
k=1 pk (πjk ) (hk − Tk x) ≤ θ
x ≥ 0, θ ∈ ℜ.
(Pν )
l = 1, . . . , s
il problema che si è risolto allo Step 1 della ν-esima iterazione del metodo Lshaped. Supponiamo che la regione ammissibile del problema (Pν ) ammetta
almeno un vertice e che (xν , θν ) sia soluzione ottima del problema tale che:
1. xν ∈ K2 ,
PN
(xν ) T
ν
ν
2.
k=1 pk (πjk ) (hk − Tk x ) ≤ θ .
Allora (xν , θν ) è soluzione ottima del problema (12).
Dim.: Dimostriamo anzitutto che risulta
N
X
(xν )
pk (πjk )T (hk − Tk xν ) = θν .
k=1
PN
(xν )
Supponiamo, per assurdo, che k=1 pk (πjk )T (hk − Tk xν ) < θν . Visto che
(xν , θν ) è un vertice del problema (Pν ), devono esistere n + 1 vincoli attivi in
(xν , θν ) linearmente indipendenti. Inoltre almeno uno dei vincoli che coinvolgono θ deve essere attivo. Sia esso
N
X
(xl̄ )
pk (πjk )T (hk − Tk xν ) = θν .
k=1
Risulterà quindi
N
X
(xl̄ )
pk (πjk )T (hk − Tk xν ) >
k=1
N
X
(xν )
pk (πjk )T (hk − Tk xν )
k=1
(xν )
che è, ovviamente, in contraddizione con il fatto che πjk è soluzione del
problema (Dqk ).
Pertanto, θν = Q(xν ) e quindi (xν , θν ) ∈ epi Q(x). Quindi, possiamo concludere
che (xν , θν ) è anche soluzione del problema (12).
2
La proposizione che precede, giustifica il criterio di arresto presente nello Step
3 del metodo secondo cui il metodo si interrompe producendo una soluzione
ottima del problema di programmazione stocastica non appena viene generata
una coppia (xν , θν ) tale che θν ≥ es+1 − Es+1 xν . Se, al contrario, la coppia
(xν , θν) è tale che θν < es+1 − Es+1 xν allora il vincolo
θ ≥ es+1 − Es+1 x
viene aggiunto al problema dello Step 1 e si riparte.
34
7.1
Esempio: soluzione del problema del venditore di giornali
Il venditore acquista i giornali dall’editore a c = 1.6 euro a copia e li rivende al
semaforo a s = 1.8 euro l’uno. I giornali non venduti possono essere rivenduti
all’editore a r = 0.8 euro a copia. Il venditore sa che per la domanda di giornali
al semaforo, sono possibili 5 esiti distinti come riportato in tabella.
k
1
2
3
4
5
pk
1/7
2/7
2/7
1/7
1/7
Dk
30
40
50
60
100
Nel seguito supporremo, senza perdita di generalità, che il venditore non possa
acquistare più di 200 giornali dall’editore. Vediamo come risolvere il problema
di programmazione stocastica del venditore di giornali usando il metodo “Lshaped”.
Iterazione 1
: ν = 1, r = s = 0.
Step 1. Dobbiamo risolvere il problema
min cx + θ
0 ≤ x ≤ 200,
ottenendo quindi la soluzione (xν , θν ) = (0, −∞).
Step 2. Notiamo che per qualunque valore x ∈ K1 , risulta sempre Y (hk −
T x) 6= ∅. Infatti, per ogni k = 1, . . . , 5, il problema
min
α1 + α2
y1 + y2 + α1 = x
y1 + y3 + α2 = Dk
yi , αj ≥ 0, i = 1, 2, 3, j = 1, 2
ammette sempre soluzione ottima
(y1 , y2 , y3 , α1 , α2 )⋆ = (min{x, Dk }, max{0, x−Dk }, Dk −min{x, Dk }, 0, 0),
quindi nessun feasibility cut verrà generato dall’algoritmo. Questo è equivalente al fatto che, per ogni x ∈ K1 risulta x ∈ K2 e cioè, K1 ⊂ K2 .
Step 3. Per k = 1, . . . , 5, risolviamo le coppie di problemi primale-duale
min −sy1 − ry2
y1 + y2 = xν = 0
(P )
y1 + y3 = Dk
yi ≥ 0, i = 1, 2, 3
35
max
(D)
0π1 + Dk π2
π1 + π2 ≤ −s
π1 ≤ −r
π2 ≤ 0.
Notiamo che la regione ammissibile duale ha solo due vertici e precisamente i punti A = (−s, 0)T e B = (−r, r − s)T . Risolvendo la coppia di problemi di sopra troviamo le soluzioni (y1k , y2k , y3k ) = (0, 0, Dk ) e
(π1k , π2k ) = (−s, 0), per k = 1, . . . , 5. Pertanto, la funzione Q(x) ammette,
nel punto xν , il subgradiente
η=−
5
X
pk (π1k
π2k )T Tk
T
= −(−s 0)
k=1
−1
0
= −s ∈ ∂Q(xν ).
Quindi, deve valere la relazione
θ ≥ Q(x) ≥ Q(xν ) + η(x − xν ) = es+1 − Es+1 x
dove
es+1 =
5
X
pk (π1k
π2k )T hk
=
k=1
e
Es+1 =
5
X
pk (−s 0)
k=1
5
X
T
0
Dk
=0
pk (π1k π2k )T Tk = −η = s.
k=1
Dal momento che θν = −∞ < es+1 − Es+1 xν = 0, imponiamo il vincolo
θ + Es+1 x ≥ es+1 , ovvero
θ + sx ≥ 0.
A questo punto, si pone s = s + 1, ν = ν + 1 e si torna allo Step 1.
Iterazione 2
: ν = 2, r = 0, s = 1.
Step 1. Dobbiamo risolvere il problema
min cx + θ
0 ≤ x ≤ 200
θ + sx ≥ 0
la cui regione ammissibile con le curve di livello della funzione obiettivo e
36
il punto di ottimo (200, −200s) sono riportati nella figura che segue.
-200
-250
-300
-350
160
170
180
190
200
210
In rosso è indicato l’optimality cut aggiunto allo Step 3 dell’iterazione
precedente.
Step 2. Come gia detto per l’iterazione 1, non è necessario aggiungere
alcun feasibility cut.
Step 3. Per k = 1, . . . , 5, risolviamo le coppie di problemi primale-duale
min
(P )
−sy1 − ry2
y1 + y2 = xν = 200
y1 + y3 = Dk
yi ≥ 0, i = 1, 2, 3
max 200π1 + Dk π2
π1 + π2 ≤ −s
(D)
π1 ≤ −r
π2 ≤ 0.
Le soluzioni ottime del primale e, rispettivamente, del duale, sono in questo
caso, per k = 1, . . . , 5: (y1k , y2k , y3k ) = (Dk , 200 − Dk , 0) e (π1k , π2k ) =
(−r, r − s). Pertanto, la funzione Q(x) ammette, nel punto xν = 200,
il subgradiente
η=−
5
X
pk (π1k π2k )T Tk = −(−r r − s)T
k=1
−1
0
= −r ∈ ∂Q(xν ).
Quindi, deve valere la relazione
θ ≥ Q(x) ≥ Q(xν ) + η(x − xν ) = es+1 − Es+1 x
dove
es+1 =
5
X
k=1
pk (π1k π2k )T hk =
5
X
pk (−r r − s)T
k=1
37
0
Dk
=
370
(r − s)
7
e
Es+1 =
5
X
pk (π1k π2k )T Tk = −η = r.
k=1
ν
Dal momento che θ = −200s < es+1 − Es+1 xν =
imponiamo il vincolo θ + Es+1 x ≥ es+1 , ovvero
θ + rx ≥
370
7 (r
− s) − 200r,
370
(r − s).
7
Poniamo s = s + 1, ν = ν + 1 e torniamo allo Step 1.
Iterazione 3
: ν = 3, r = 0, s = 2.
Step 1. Dobbiamo risolvere il problema
min cx + θ
0 ≤ x ≤ 200
θ + sx ≥ 0
θ + rx ≥ 370
7 (r − s)
la cui regione ammissibile con le curve di livello della funzione obiettivo e
il punto di ottimo (370/7, −s370/7) sono riportati nella figura che segue.
0
50
100
150
200
-50
-100
-150
-200
In rosso è indicato l’optimality cut aggiunto nell’iterazione 1 e in blu quello
aggiunto nell’iterazione 2. Notiamo che, coerentemente con il loro nome,
optimality cut, i vincoli che via via andiamo ad aggiungere allo Step 3, non
fanno altro che tagliar fuori dalla regione ammissibile i punti di ottimo
trovati allo Step 1 dell’iterazione precedente.
Step 2. Non viene aggiunto nessun feasibility cut.
38
Step 3. Per k = 1, . . . , 5, risolviamo le coppie di problemi primale-duale
min
(P )
−sy1 − ry2
y1 + y2 = xν = 370/7
y1 + y3 = Dk
yi ≥ 0, i = 1, 2, 3
max π1 370/7 + Dk π2
π1 + π2 ≤ −s
(D)
π1 ≤ −r
π2 ≤ 0.
Le soluzioni ottime del primale e, rispettivamente, del duale, sono riportate
nella tabella che segue:
k
1
2
3
4
5
pk
1/7
2/7
2/7
1/7
1/7
yk
(30, 370/7 − 30, 0)
(40, 370/7 − 40, 0)
(50, 370/7 − 50, 0)
(370/7, 0, 60 − 370/7)
(370/7, 0, 100 − 370/7)
Dk
30
40
50
60
100
πk
(−r, r − s)
(−r, r − s)
(−r, r − s)
(−s, 0)
(−s, 0)
Pertanto, la funzione Q(x) ammette, nel punto xν = 370/7, il subgradiente
η=−
5
X
k=1
5
2
pk (π1k π2k )T Tk = − r − s ∈ ∂Q(xν ).
7
7
Quindi, deve valere la relazione
θ ≥ Q(x) ≥ Q(xν ) + η(x − xν ) = es+1 − Es+1 x
dove
es+1 =
5
X
pk (π1k π2k )T hk = 30(r − s)
k=1
e
Es+1 =
5
X
pk (π1k π2k )T Tk = −η =
k=1
5
2
r + s.
7
7
370
Dal momento che θ = −s370/7 < es+1 − Es+1 xν = 30(r − s) − 5r+2s
7
7 ,
imponiamo il vincolo θ + Es+1 x ≥ es+1 , ovvero
ν
θ+
5r + 2s
x ≥ 30(r − s).
7
Poniamo s = s + 1, ν = ν + 1 e torniamo allo Step 1.
Iterazione 4
: ν = 4, r = 0, s = 3.
39
Step 1. Dobbiamo risolvere il problema
min
cx + θ
0 ≤ x ≤ 200
θ + sx ≥ 0
θ + rx ≥ 370
7 (r − s)
θ + 5r+2s
x
≥ 30(r − s)
7
la cui regione ammissibile con le curve di livello della funzione obiettivo e
il punto di ottimo (42, −75.6) sono riportati nella figura che segue.
50
0
100
150
200
-50
-100
-150
-200
In rosso, blu e verde sono indicati gli optimality cut aggiunti, rispettivamente, nelle iterazioni 1, 2 e 3.
Step 2. Non viene aggiunto nessun feasibility cut.
Step 3. Per k = 1, . . . , 5, risolviamo le coppie di problemi primale-duale
min −sy1 − ry2
y1 + y2 = xν = 42
(P )
y1 + y3 = Dk
yi ≥ 0, i = 1, 2, 3
max
(D)
42π1 + Dk π2
π1 + π2 ≤ −s
π1 ≤ −r
π2 ≤ 0.
Le soluzioni ottime del primale e, rispettivamente, del duale, sono riportate
nella tabella che segue:
k
1
2
3
4
5
pk
1/7
2/7
2/7
1/7
1/7
Dk
30
40
50
60
100
yk
(30, 12, 0)
(40, 2, 0)
(42, 0, 8)
(42, 0, 18)
(42, 0, 58)
40
πk
(−r, r − s)
(−r, r − s)
(−s, 0)
(−s, 0)
(−s, 0)
Pertanto, la funzione Q(x) ammette, nel punto xν = 42, il subgradiente
η=−
5
X
k=1
3
4
pk (π1k π2k )T Tk = − r − s ∈ ∂Q(xν ).
7
7
Quindi, deve valere la relazione
θ ≥ Q(x) ≥ Q(xν ) + η(x − xν ) = es+1 − Es+1 x
dove
es+1 =
5
X
pk (π1k π2k )T hk =
k=1
e
Es+1 =
5
X
110
(r − s)
7
pk (π1k π2k )T Tk = −η =
k=1
3
4
r + s.
7
7
ν
Dal momento che θ = −75.6 < es+1 − Es+1 xν =
imponiamo il vincolo θ + Es+1 x ≥ es+1 , ovvero
θ+
110
7 (r
− s) −
3r+4s
7 42,
110
3r + 4s
x≥
(r − s).
7
7
Poniamo s = s + 1, ν = ν + 1 e torniamo allo Step 1.
Iterazione 5
: ν = 5, r = 0, s = 4.
Step 1. Dobbiamo risolvere il problema
min cx + θ
0 ≤ x ≤ 200
θ + sx ≥ 0
θ + rx ≥ 370
7 (r − s)
θ + 5r+2s
x
≥ 30(r − s)
7
θ + 3r+4s
x
≥ 110
7
7 (r − s)
la cui regione ammissibile con le curve di livello della funzione obiettivo e
41
il punto di ottimo (110/3, −66) sono riportati nella figura che segue.
-40
-60
-80
-100
-120
20
40
60
80
100
Step 2. Non viene aggiunto nessun feasibility cut.
Step 3. Per k = 1, . . . , 5, risolviamo le coppie di problemi primale-duale
min
(P )
−sy1 − ry2
y1 + y2 = xν = 110/3
y1 + y3 = Dk
yi ≥ 0, i = 1, 2, 3
max π1 110/3 + Dk π2
π1 + π2 ≤ −s
(D)
π1 ≤ −r
π2 ≤ 0.
Le soluzioni ottime del primale e, rispettivamente, del duale, sono riportate
nella tabella che segue:
k
1
2
3
4
5
pk
1/7
2/7
2/7
1/7
1/7
Dk
30
40
50
60
100
yk
(30, 20/3, 0)
(110/3, 0, 10/3)
(110/3, 0, 40/3)
(110/3, 0, 70/3)
(110/3, 0, 190/3)
πk
(−r, r − s)
(−s, 0)
(−s, 0)
(−s, 0)
(−s, 0)
Pertanto, la funzione Q(x) ammette, nel punto xν = 110/3, il subgradiente
η=−
5
X
k=1
1
6
pk (π1k π2k )T Tk = − r − s ∈ ∂Q(xν ).
7
7
Quindi, deve valere la relazione
θ ≥ Q(x) ≥ Q(xν ) + η(x − xν ) = es+1 − Es+1 x
42
dove
es+1 =
5
X
pk (π1k π2k )T hk =
k=1
e
Es+1 =
5
X
30
(r − s)
7
pk (π1k π2k )T Tk = −η =
k=1
1
6
r + s.
7
7
Dal momento che θν = −66 < es+1 − Es+1 xν = −30/7 − (11.6/7)110/3 =
−65, imponiamo il vincolo θ + Es+1 x ≥ es+1 , ovvero
30
r + 6s
x≥
(r − s).
7
7
Poniamo s = s + 1, ν = ν + 1 e torniamo allo Step 1.
θ+
Iterazione 6
: ν = 6, r = 0, s = 5.
Step 1. Dobbiamo risolvere il problema
min cx + θ
0 ≤ x ≤ 200
θ + sx ≥ 0
θ + rx ≥ 370
7 (r − s)
θ + 5r+2s
x
≥ 30(r − s)
7
θ + 3r+4s
x
≥ 110
7
7 (r − s)
r+6s
θ + 7 x ≥ 30
7 (r − s)
la cui regione ammissibile con le curve di livello della funzione obiettivo e
il punto di ottimo (40, −494/7) sono riportati nella figura che segue.
-40
-45
-50
-55
-60
-65
-70
-75
25
30
35
43
40
Step 2. Non viene aggiunto nessun feasibility cut.
Step 3. Per k = 1, . . . , 5, risolviamo le coppie di problemi primale-duale
min −sy1 − ry2
y1 + y2 = xν = 40
(P )
y1 + y3 = Dk
yi ≥ 0, i = 1, 2, 3
max
(D)
40π1 + Dk π2
π1 + π2 ≤ −s
π1 ≤ −r
π2 ≤ 0.
Le soluzioni ottime del primale e, rispettivamente, del duale, sono riportate
nella tabella che segue:
k
1
2
3
4
5
pk
1/7
2/7
2/7
1/7
1/7
Dk
30
40
50
60
100
yk
(30, 10, 0)
(40, 0, 0)
(40, 0, 10)
(40, 0, 20)
(40, 0, 60)
πk
(−r, r − s)
(−s, 0)
(−s, 0)
(−s, 0)
(−s, 0)
Pertanto, la funzione Q(x) ammette, nel punto xν = 40, il subgradiente
η=−
5
X
k=1
1
6
pk (π1k π2k )T Tk = − r − s ∈ ∂Q(xν ).
7
7
Quindi, deve valere la relazione
θ ≥ Q(x) ≥ Q(xν ) + η(x − xν ) = es+1 − Es+1 x
dove
es+1 =
5
X
pk (π1k π2k )T hk =
k=1
e
Es+1 =
5
X
30
(r − s)
7
pk (π1k π2k )T Tk = −η =
k=1
6
1
r + s.
7
7
ν
Dal momento che θ = −494/7 = es+1 − Es+1 xν = −30/7 − (11.6/7)40,
concludiamo che la soluzione corrente (xν , θν ) = (40, −494/7) è ottima
per il problema di programmazione stocastica.
8
Indicatori di validità e attendibilità della Programmazione Stocastica
Come gia abbiamo avuto modo di osservare in precedenza, risolvere un problema di programmazione stocastica è, in taluni casi, equivalente alla soluzione di
44
un problema di grandi dimensioni cioè con un elevato numero di variabili e/o
di vincoli. Pertanto, in generale, determinare la soluzione di un problema di
programmazione stocastica lineare è molto costoso. Ha quindi senso chiedersi
quando e quanto sia conveniente abbandonare il problema stocastico e risolvere,
invece, un problema più semplice. In particolare, è lecito chiedersi in che misura
approcci più semplici al problema, come per esempio quello che prevede di sostituire alla v.a. il suo valor medio, forniscano soluzioni che si discostano dall’ottimo o fine a che punto questi approcci “alternativi” non siano completamente
inaccurati.
Una risposta a queste domande può essere fornita da due importanti indicatori
di bontà dell’approccio stocastico che sono:
1. l’EVPI (Expected Value of Perfect Information) ovvero il valore atteso in
condizioni di informazione completa e,
2. l’VSS (Value of Stochastic Solution) ovvero il valore della soluzione stocastica.
Al fine di definire chiaramente questi due importanti indicatori, introduciamo
brevemente alcune definizioni che ci faranno comodo nel seguito.
Sia ω ∈ Ω una v.a. discreta ovvero che può avere solo un numero finito N di realizzazioni ω1 , . . . , ωN . Sia f (x⋆ (ωi ), ωi ) il valore ottimo del problema di programmazione stocastica quando si fissa la v.a. ω al valore ωi , i = 1, . . . , N . Percui,
data una realizzazione ωi della v.a. ω ∈ Ω, associamo ad essa il valore all’ottimo del corrispondente problema di programmazione matematica f (x⋆ (ωi ), ωi ).
f (x⋆ (ω), ω) è pertanto una v.a. composta per cui è possibile definire il valore
atteso. Infatti, definiamo soluzione wait-and-see
W S = IEω [f (x⋆ (ω), ω)],
il valore atteso della v.a. f (x⋆ (ω), ω). Definiamo, invece, soluzione here-and-now
come il valore della soluzione del problema di programmazione stocastica
HN = min IEω [f (x, ω)].
x
A questo punto, definiamo l’EVPI (valore atteso in condizioni di informazione
completa) come
EV P I = HN − W S.
Nel caso dell’esempio del venditore di giornali, abbiamo che, come risulta dalla
tebella che segue, W S = −10.57142
i
1
2
3
4
5
pi
1/7
2/7
2/7
1/7
1/7
Di
30
40
50
60
100
x⋆ (Di )
30
40
50
60
100
45
f (x⋆ (Di ), Di )
-6
-8
-10
-12
-20
Risulta, inoltre, x⋆ = 40,
Q(x⋆ ) =
−1.8D1 − 0.8(x⋆ − D1 )
1.8x⋆
−6
7
7
e quindi HN = 1.6x⋆ + Q(x⋆ ) = −6.57142. Otteniamo cosı̀ EV P I = 4. Questo
valore è pari alla quantità di euro che il venditore sarebbe disposto a pagare
ogni mattina pur di sapere quanti giornali potrà vendere al semaforo.
Quando ci si trova a dover risolvere un problema di ottimizzazione con incertezza, molto spesso si ritiene di poter ottenere una buona approssimazione della
soluzione ottima semplicemente risolvendo il problema
EV = min f (x, ω̄),
x
dove si è sostituita la v.a. ω con il suo valore atteso ω̄ = IEω [ω]. La soluzione
di questo problema è nota come soluzione EV (expected value) ovvero soluzione
di valor medio. L’indicatore VSS è quello che ci dice quanto la soluzione EV
sia lontana dall’essere ottima per il problema di programmazione stocastica. In
particolare, una volta nota la soluzione ottima x⋆ (ω̄) del problema EV, definiamo
EEV = IEω [f (x⋆ (ω̄), ω)] = cT x⋆ (ω̄) + Q(x⋆ (ω̄))
ovvero il valore atteso quando si usa la soluzione di valor medio. Il VSS è definito
come
V SS = EEV − HN.
Di nuovo, nel caso del venditore di giornali, abbiamo che x⋆ (D̄) = 370/7 e
EV = −10.57142. Come risulta dalla tabella che segue abbiamo
i
1
2
3
4
5
pi
1/7
2/7
2/7
1/7
1/7
Di
30
40
50
60
100
cx⋆ (D̄) + Q(x⋆ (D̄), Di )
1.6x⋆ (D̄) − 1.8D1 − 0.8(x⋆ (D̄) − D1 )
1.6x⋆ (D̄) − 1.8D2 − 0.8(x⋆ (D̄) − D2 )
1.6x⋆ (D̄) − 1.8D3 − 0.8(x⋆ (D̄) − D3 )
1.6x⋆ (D̄) − 1.8x⋆ (D̄)
1.6x⋆ (D̄) − 1.8x⋆ (D̄)
f (x⋆ (D̄), Di )
12.28571
2.28571
-7.71429
-10.57143
-10.57143
EEV = IED [f (x⋆ (D̄), D)] = −2.81634 e quindi V SS = 3.75508 che ci dice
esattamente di quanto la soluzione del problema EV si discosta dalla soluzione
HN.
Consideriamo ora l’esempio dell’azienda agricola. Come gia abbiamo visto, in
questo caso HN = −108390. Calcoliamo ora WS. Nel caso in cui ω = ω1 (stagione sotto la media), risulta z(x(ω1 ), ω1 ) = −59950; quando ω = ω2 (stagione
nella media), z(x(ω2 ), ω2 ) = −118600; quando ω = ω3 (stagione sopra la media), z(x(ω3 ), ω3 ) = −167667; cosicche W S = z(x(ω1 ), ω1 )/3 + z(x(ω2), ω2 )/3 +
z(x(ω3 ), ω3 )/3 = −115405.56. Quindi otteniamo
EV P I = HN − W S = 7015.6
46
il che significa che l’azienda sarebbe disposta a pagare 7015.6 euro ogni anno
per poter conoscere quale sarà l’andamento della stagione.
Calcoliamo ora l’EV dell’esempio dell’azienda agricola che otteniamo semplicemente risolvendo il problema in cui si è fissata la v.a. al suo valor medio il che,
ovviamente, equivale a risolvere il problema nel caso di stagione nella media per
cui otteniamo x(ω̄) = (120, 80, 300)T . Possiamo ora calcolare l’EEV che risulta
EEV = −107240 per cui otteniamo
V SS = EEV − HN = 1150.
9
Richiami sulle funzioni convesse
¯ gli insiemi estesi ℜ ∪ {+∞} e ℜ ∪ {−∞, +∞}, rispettivaIndichiamo con ℜ̄ e ℜ̄
¯ con S ⊆ ℜn . L’insieme
mente. Sia f : S → ℜ̄
{(x, µ) : x ∈ S, µ ∈ ℜ, f (x) ≤ µ}
è detto epigrafo di f ed è indicato epi f . Diciamo che f è una funzione convessa
su S se epi f è un sottoinsieme convesso di ℜn+1 . Il dominio effettivo della
funzione f su S è la proiezione su ℜn dell’epigrafo di f ovvero
dom f = {x : ∃µ, (x, µ) ∈ epi f } = {x : f (x) < +∞}.
Notiamo che data una funzione f su S convessa, è sempre possibile ottenere
una funzione f ′ definita su ℜn e ancora convessa. Per fare questo è sufficiente
considerare la funzione
f (x) se x ∈ S,
′
f (x) =
+∞ se x 6∈ S.
Dal momento che le funzioni che stiamo calcolando hanno valori nello spazio
¯ , è necessario fornire delle regole che specifichino il risultato delle operesteso ℜ̄
azioni aritmetiche fondamentali quando sono coinvolti i simboli +∞ o −∞. In
particolare, adotteremo le seguenti regole:
1. α + ∞ = ∞ + α = ∞ per ogni −∞ < α ≤ ∞;
2. α − ∞ = −∞ + α = −∞ per ogni −∞ ≤ α < ∞;
3. α∞ = ∞α = ∞, α(−∞) = (−∞)α = −∞ per ogni 0 < α ≤ ∞;
4. α∞ = ∞α = −∞, α(−∞) = (−∞)α = ∞ per ogni −∞ ≤ α < 0;
5. 0∞ = ∞0 = 0 = 0(−∞) = (−∞)0, −(−∞) = ∞;
6. ∞ − ∞ = −∞ + ∞ = ∞;
7. inf ∅ = +∞, sup ∅ = −∞.
47
Una funzione convessa f è detta propria se f (x) < +∞ per almeno un x e
f (x) > −∞ per ogni x. Una funzione convessa che non è propria è detta
impropria. In particolare se f è una funzione convessa finita definita sull’insieme
convesso C, la funzione f ′ definita su tutto ℜn
f (x) se x ∈ C,
′
f (x) =
+∞ se x 6∈ C.
è una funzione convessa e propria.
Proposizione 5 Sia f : C → ℜ̄, dove C ⊆ ℜn . f è una funzione convessa su
C se e solo se
1. C è un insieme convesso;
2. comunque scelti x, y ∈ C risulta
f ((1 − λ)x + λy) ≤ (1 − λ)f (x) + λf (y),
0 < λ < 1.
¯ . f è una funzione convessa se e solo se
Proposizione 6 Sia f : ℜn → ℜ̄
n
comunque presi due punti x, y ∈ ℜ per cui esistono due scalari α, β < +∞ tali
che f (x) < α e f (y) < β allora
f ((1 − λ)x + λy) < (1 − λ)α + λβ,
0 < λ < 1.
Proposizione 7 (Disuguaglianza di Jensen) Sia f : ℜn → ℜ̄. Allora f è
convessa se e solo se
f (λ1 x1 + . . . + λm xm ) ≤ λ1 f (x1 ) + . . . + λm f (xm ),
comunque scelti i vettori xi ∈ ℜn , i = 1, . . . , m, dove gli scalari λi sono tali che
λi ≥ 0, i = 1, . . . , m e
m
X
λi = 1.
i=1
10
Definizioni
Dato S ⊂ ℜn , diciamo che S è un insieme convesso quando, comunque presi
due punti x, y ∈ S, risulta, per ogni scalare λ ∈ [0, 1]
λx + (1 − λ)y ∈ S.
Sia S ⊂ ℜn un insieme convesso. Diciamo che la funzione f : ℜn → ℜ è
convessa su S quando, comunque presi due punti x, y ∈ S, risulta, per ogni
scalare λ ∈ [0, 1]
f (λx + (1 − λ)y) ≤ λf (x) + (1 − λ)f (y).
48
Data una funzione f : ℜn → ℜ convessa sul’insieme convesso X e un punto
xo ∈ X, diciamo che un vettore η ∈ ℜn è un subgradiente di f (x) in xo
quando:
f (x) ≥ f (xo ) + η T (x − xo ), ∀ x ∈ X.
Data una funzione f : ℜn → ℜ convessa sul’insieme convesso X e un punto
xo ∈ X, definiamo il subdifferenziale di f (x) in xo come segue:
∂f (xo ) = conv{η ∈ ℜn : f (x) ≥ f (xo ) + η T (x − xo ),
∀ x ∈ X}.
Una funzione f (x) si dice positivamente omogenea quando, comunque preso
uno scalare λ ≥ 0, risulta f (λx) = λf (x).
Sia S ⊂ ℜn un insieme convesso.
dell’insieme S la funzione σS : ℜn → ℜ
Definiamo funzione di supporto
σS (x) = sup xT y.
y∈S
Indice
1 Premessa
1
2 Nozioni preliminari
2.1 Formulazione deterministica . . . .
2.2 Osservazione e ottimizzazione . . .
2.3 Ottimizzazione e osservazione . . .
2.4 Esempio: il problema del venditore
.
.
.
.
1
2
2
3
3
3 Fondamenti di Programmazione Stocastica
3.1 Esempio: financial planning and control . . . . . . . . . . . . . .
8
10
4 Programmazione Stocastica lineare a due stadi
13
. . . . . .
. . . . . .
. . . . . .
di giornali
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
5 Esempi di modelli stocastici lineari a due stadi
15
5.1 Il problema dell’azienda agricola . . . . . . . . . . . . . . . . . . 15
5.2 Gestione degli investimenti per un impianto di distribuzione elettrica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
6 Programmazione stocastica lineare a due stadi con ricorsione
fissa
22
6.1 Esempio: il problema del venditore di giornali . . . . . . . . . . . 26
7 L-shaped method
28
7.1 Esempio: soluzione del problema del venditore di giornali . . . . 35
49
8 Indicatori di validità e attendibilità della Programmazione Stocastica
44
9 Richiami sulle funzioni convesse
47
10 Definizioni
48
Riferimenti bibliografici
[1] J.R. Birge and F. Louveaux “Introduction to Stochastic Programming”,
Springer Series in Operations Research and Financial Engineering. Springer
(1997)
[2] R.T. Rockafellar “Optimization under Uncertainty”, Lecture notes,
University of Washington.
50