10 Ein Dekompositionsverfahren für stochastische
Optimierungsprobleme mit rekombinierenden
Szenariobäumen
Küchler, C. ; Vigerske, S.
Kurzfassung: Dieses Kapitel stellt ein Dekompositionsverfahren für lineare mehrstufige stochastische Programme vor, welches auf rekombinierenden Szenariobäumen beruht. Rekombinierende
Szenariobäume sind numerisch vorteilhaft, da sie es erlauben, ein exponentielles Wachstum derKnotenzahl des Optimierungsproblems mit der Anzahl der Zeitstufen zu verhindern. Es wird gezeigt, wie sich mittels eines geeigneten Dekompositionsansatzes dieser Vorteil auch unter Gegenwart zeitkoppelnder Restriktionen bewahren lässt. Das vorgestellte Verfahren ist verwandt mit der
Nested Benders Dekomposition und nutzt die spezielle Struktur rekombinierender Szenariobäume
zur simultanen Generierung von Schnittebenen und zum “dynamischen Rekombinieren” von Szenarien. Weiter wird ein Verfahren zur Erstellung von rekombinierenden Bäumen vorgestellt.
Schlagworte: Mehrstufige Stochastische Optimierung, Nested Benders Dekomposition, Rekombinierende Szenariobäume

10.1

Einleitung

In vielen Anwendungen ist es notwendig, Entscheidungen zu treffen, deren Konsequenzen von unbekannten, erst in der Zukunft beobachtbaren Größen abhängen. Die Suchenach einer optimalen Entscheidung kann als ein Optimierungsproblem betrachtet werden, das Unsicherheit hinsichtlich eines oder mehrerer Modellparameter beinhaltet. Oft
sind gewisse statistische Informationen über diese “stochastischen” Parameter verfügbar,
beispielsweise aus historischen oder simulierten Daten gewonnene empirische Wahrscheinlichkeitsverteilungen. Dann ist es möglich, diese Parameter durch Zufallsvariablenoder stochastische Prozesse zu ersetzen. Auf diese Weise erhält man ein stochastisches
Optimierungproblem [19], dessen Lösungen Unterstützung bei der konkreten Entscheidungsfindung bieten können. Ein Beispiel hierfür ist das in Kapitel 8 betrachtete Modell,bei dem die eingespeiste Windenergie, die Strompreise und die Stromnachfrage als stochastische Größen gegeben sind.
Zwei- und mehrstufige stochastische Optimierungsprobleme zeichnen sich dadurch aus,
dass die in der ersten (Zeit-)Stufe zu treffenden Entscheidungen ohne Kenntnis der Realisierung der stochastischen Größen zu treffen sind, während die in der zweiten Stufe zu
treffenden Entscheidungen von den bis zu diesem Zeitpunkt beobachteten stochastischenGrößen abhängen dürfen. Mehrstufige Probleme beinhalten weitere Zeitstufen mit erst
dann beobachtbaren zufälligen Größen. Ihre Lösungen bestehen aus einer Abfolge vonEntscheidungen, die jeweils wieder nur von den bis zum jeweiligen Entscheidungszeitpunkt beobachtbaren Größen abhängen dürfen.

Mehrstufige stochastische Optimierungsprobleme können im Allgemeinen nur durch numerische Methoden gelöst werden. Dazu wird der zugrunde liegende stochastische Prozess, welcher beispielsweise durch ein Zeitreihenmodell gegeben sein kann, durch einen
Prozess approximiert, der nur endlich viele Werte annimmt. Letzterer kann durch einen
Szenariobaum repräsentiert werden (vgl. Kapitel 11). Da die Anzahl der Szenarien in
einem solchen Baum exponentiell mit der Länge des betrachteten Zeithorizonts wachsen
kann, vgl. Abbildung 10.1, werden zumeist nur Modelle mit einer geringen Anzahl vonZeitstufen oder einer deutlich reduzierten Anzahl von Szenarien betrachtet.
Eine Reduktion der Zeitstufen führt dazu, dass Entscheidungen in diesem Modell nur
zu sehr wenigen Zeitpunkten getroffen werden können. Dies wird, insbesondere für dynamische Probleme, dem zu modellierenden Problem jedoch oft nicht gerecht. Daher
besteht ein weitverbreiteter Ansatz zur Bestimmung akzeptabler Lösungen darin, einenrollierenden Planungshorizont zu verwenden. Dazu wird das Optimierungsproblem ineine Menge von Teilproblemen mit kürzeren Zeithorizonten zerlegt und diese sukzessivegelöst. Ist das mehrstufige Problem zeitlich zerlegbar, d.h. besteht beispielsweise keine
Kopplung zwischen den Entscheidungen eines Tages und jenen des Folgetages, so liefert
die Optimierung mit rollierendem Horizont tatsächlich optimale Lösungen. Allerdings
ist die genaue Berücksichtigung zeitkoppelnder Restriktionen, wie beispielsweise Speicherstände oder Mindestbetriebs- und -stillstandszeiten von Anlagen, problematisch.
Um den zugrunde liegenden stochastischen Prozess möglichst gut durch eine deutlich
reduzierte Anzahl von Szenarien zu approximieren, wurden Methoden zur Szenarioreduktion entwickelt, vgl. [10]. Diese erlauben es, durch Verwendung einer geringeren Anzahl von ausgewählten Szenarien das Optimierungsproblem auch für längere Zeiträume
zu lösen. Hierbei ist jedoch zu berücksichtigen, dass bei einer deutlichen Reduktion der
Szenarienanzahl das Modell die Unsicherheit bezüglich der stochastischen Daten zunehmend unterschätzt.
Ein weiterer Ansatz zur Reduzierung der Dimension des Optimierungsproblems besteht
in der Verwendung rekombinierender Szenariobäume. Das hierfür bestbekannte Beispiel
ist das Binomialmodell von Cox, Ross und Rubinstein [2] zur Berechnung von Optionspreisen. Für ein Problem mit T Zeitstufen verringert hier das Rekombinieren von Szenarien die Anzahl der Knoten von 2T −1 auf T (T +1)/2 und ermöglicht so erst die Lösungvon Problemen mit längeren Zeithorizonten. Rekombinierende Szenariobäume wurden
des weiteren bisher unter anderem zur Bewertung der Kosten von fluktuierender Windeinspeisung [21] und zur Bewertung von Gasspeichern [8] eingesetzt. Allerdings ist bei
der stochastischen Optimierung auf einem rekombinierenden Baum die präzise Abbildung von Nebenbedingungen, welche verschiedene Zeitstufen koppeln, nicht möglich. In[21] werden zeitkoppelnden Restriktionen daher nur näherungsweise erfüllt. In [8] wird
der zulässige Bereich der gekoppelten Entscheidungsvariablen diskretisiert und ein auf

dynamischer Programmierung beruhender Lösungsansatz verwendet. Im vorliegenden
Kapitel wird ein anderer Ansatz vorgestellt, der auf der folgenden Beobachtung beruht.
Die Kostenfunktion in einem Knoten des Szenariobaums beschreibt die erwarteten, minimal erreichbaren zukünftigen Kosten in Abhängigkeit der Entscheidungen aus vorhergehenden Zeitstufen und des Wertes der stochastischen Parameter in diesem Knoten.
Die erwarteten zukünftigen Kosten sind von dem im jeweiligen Knoten beginnenden
Teilbaum bestimmt. Gibt es zu einem Zeitpunkt Knoten im Szenariobaum, deren Teilbäume übereinstimmen, so sind die jeweiligen Kostenfunktionen identisch. Zerlegt man
das stochastische Problem zu diesen Zeitpunkten im Sinne der Nested Benders Dekomposition [1, 18] und approximiert die dort auftretenden Kostenfunktionen, so muss fürübereinstimmende Kostenfunktionen nur eine einzige Approximation bestimmt werden.
Dieses als “cut sharing” bezeichnete Prinzip wurde bereits von Infanger und Morton[13] und Pereira und Pinto [16] angewandt und führt durch die Ausnutzung von übereinstimmenden Kostenfunktionen zu einer beschleunigten Konvergenz des Lösungsalgorithmus. Allerdings gibt es mit wachsendem Zeithorizont eine exponentiell wachsende Anzahl von Knoten mit übereinstimmenden Kostenfunktionen, welche jeweils in den knotenabhängigen Entscheidungsvariablen evaluiert werden müssen. Während diese Knotenmenge in [16] auf eine zufällige Auswahl reduziert wird, wird in dem in diesem Kapitel vorgestellten Algorithmus folgende Eigenschaft genutzt: Findet man Knoten in denen die entsprechenden Entscheidungsvariablen ähnliche Werte besitzen, so genügt esdie (Lipschitz-stetige) Kostenfunktion in nur einer einzigen dieser Kontrollen auszuwerten. Damit können im Verlauf des Lösungsverfahrens Knoten mit gleichem Teilbaum“rekombiniert” werden, so lange die (jeweils besten bekannten) Entscheidungen in diesen Knoten ähnlich sind. Dieses Vorgehen kann als “dynamisches Rekombinieren” von
Szenarien interpretiert werden und ermöglicht, ein exponentielles Wachstum der Anzahl
von Kostenfunktionsauswertungen zu verhindern. Somit lassen sich stochastische Optimierungsprobleme mit zeitkoppelnden Restriktionen, vielen Zeitstufen und einer großen
Anzahl von Szenarien lösen.
Das vorliegende Kapitel gliedert sich wie folgt. In Abschnitt 10.2 werden die betrachtete Klasse von linearen mehrstufigen stochastischen Optimierungsproblemen vorgestellt
und das zeitliche Dekompositionsverfahren sowie das Konzept rekombinierender Szenariobäume skizziert. In Abschnitt 10.3 wird eine Modifikation der Benders Dekompositionfür Optimierungsprobleme entwickelt, die einen Szenariobaum beinhalten, in welchem
Szenarien in nur einer einzelnen Zeitstufe rekombiniert werden. Diese Technik bildet die
Grundlage für die Erweiterung auf den in Abschnitt 10.4 vorgestellten allgemeinen Fall
mit mehrfachem Rekombinieren. Ein Algorithmus zur Konstruktion von rekombinierenden Szenariobäumen aus einer gegebenen Menge von Szenarien wird in Abschnitt 10.5
diskutiert.

10.2

Problemformulierung und zeitliche Dekomposition

10.2.1

Lineare mehrstufige stochastische Probleme

Es sei ξ = (ξt )t=1,...,T ein Rs -wertiger zeitdiskreter stochastischen Prozess1 mit Zeithorizont T ∈ N auf einem Wahrscheinlichkeitsraum (Ω, F, P). Für t = 1, . . . , T sei
mit ξ[t] der Vektor (ξ1 , . . . , ξt ) bezeichnet. Problem (P) stellt ein allgemeines lineares
mehrstufiges stochastisches Optimierungsproblem (“MSP”) dar:
" T
#
X
min E
bt (ξt ), xt (ξ[t] )

(P)

t=1

so dass xt : Rs·t → Xt messbar,t = 1, . . . , T,

At,0 (ξt )xt (ξ[t] ) + At,1 (ξt )xt−1 (ξ[t−1] ) = ht (ξt ),t = 2, . . . , T.

Dabei bezeichnet xt die im Zeitpunkt t zu treffenden Entscheidungen; die entsprechenden zulässigen Mengen Xt ⊂ Rr seien abgeschlossene beschränkte Polyeder. Die Entscheidung xt = xt (ξ[t] ) darf von den (bis zum Zeitpunkt t) aus der Beobachtung von ξ[t]erhaltenen Informationen abhängen. Die Abbildungen ht (·), bt (·) und At,0 (·), At,1 (·)
sind vektor- bzw. matrix-wertig von geeigneter Dimension. Ohne Beschränkung der Allgemeinheit sei angenommen dass Kopplungen nur zwischen aufeinanderfolgenden Zeitschritten auftreten2 .
Um eine numerische Lösung von (P) zu ermöglichen wird angenommen, dass der Prozess
ξ Werte in einer endlichen Menge von Szenarien annimmt und damit durch einen Szenariobaum repräsentiert werden kann, vgl. [3]. Eine Baumstruktur entsteht dabei dadurch,
dass jeweils gewisse Szenarien bis zu einem Zeitpunkt t ∈ {2, . . . , T − 1} übereinstimmen. Die Menge aller Szenarien die jeweils bis zum Zeitpunkt t übereinstimmen, bilden
einen Knoten des Szenariobaumes. Es sei nt die Anzahl der Knoten zum Zeitpunkt t;diese seien mit i = 1, . . . , nt indiziert. Der Wert des Prozesses ξ bis zum Knoten i zur(i)Zeit t wird mit ξ[t] bezeichnet. Die Szenarien entsprechen den Werten des Prozesses bis(i)zu den Endknoten, d.h. sie sind gegeben durch ξ[T ] , i = 1, . . . , nT . Die Entscheidungxt (·) zum Zeitpunkt t ist konstant auf jedem Knoten in t, insbesondere ist xt (·) durch die
(i)

(i)

Werte xt := xt (ξ[t] ), i = 1, . . . , nt , definiert.
10.2.2

Dynamische Formulierung

Die diskrete Struktur des Prozesses ξ erlaubt es, Problem (P) als ein (großes) deterministisches lineares Optimierungsproblem zu formulieren. Allerdings ist eine direkte nume1 Im folgenden bezeichnen fett gedruckte Variablen (ξ, ξ , . . .) zufällige Größen, während die Variablent(ξ, ξt , . . .) für mögliche Werte dieser Größen stehen.
2 Gibt es Kopplungen von nicht aufeinanderfolgenden Zeitstufen, so kann das Optimierungsproblem durch
die Einführung zusätzlicher Entscheidungsvariablen auf die Form (P) gebracht werden.

rische Lösung dieses deterministischen Äquivalents durch verfügbare Löser für lineare
Programme in der Regel nicht ohne Weiteres möglich, da mit wachsendem Zeithorizont
die Anzahl der Knoten im Szenariobaum und damit die Größe des Problems (P) exponentiell wächst. Daher sind Dekompositionsverfahren notwendig, welche die besondere
Struktur von (P) ausnutzen können. Verschiedene solche an stochastische Optimierungsprobleme angepasste Lösungsmethoden werden in [18] vorgestellt.
Im Folgenden wird kurz eine einfache zeitliche Dekomposition vorgestellt. Dazu betrachtet man eine Zeitstufe R ∈ {2, . . . , T − 1} und definiert die Kostenfunktion zur Zeit R
(i)

(i)

und im Zustand (xR , ξ[R] ) ∈ Rr × Rs·R als
(i)

(i)

QR (xR , ξ[R] ) :=
" T
#
X
(i)
min E
hbt (ξt ), xt (ξ[t] )i ξ[R] = ξ[R]

(QR )

t=R+1

so dass xt : Rs·t → Xt messbar,t = R + 1, . . . , T,(i)

AR+1,0 (ξR+1 )xR+1 (ξ[R+1] ) = hR+1 (ξR+1 ) − AR+1,1 (ξR+1 )xR ,At,0 (ξt )xt (ξ[t] ) + At,1 (ξt )xt−1 (ξ[t−1] ) = ht (ξt ),(10.1)

t = R + 2, . . . , T.
(i)

Dabei bezeichnet für eine Funktion f der Term E[f (ξt+1 , . . . , ξT )|ξ[R] = ξ[R] ] die bedingte Erwartung von f (ξt+1 , . . . , ξT ), gegeben dass der Zufallsvektor ξ[R] den Wert
(i)

(i)

(i)

ξ[R] annimmt. Die Kostenfunktion QR (xR , ξ[R] ) gibt die erwarteten, minimal erreichbaren Kosten an, die nach dem Zeitpunkt R anfallen, wenn bis zum Zeitpunkt R der Wert(i)(i)
ξ[R] beobachtet und in R die Entscheidung xR getroffen wurde.
Mit Hilfe dieser Kostenfunktion kann das Problem (P) in der folgenden dynamischen
Form geschrieben werden:

min E

" R
X

#
hbt (ξt ), xt (ξ[t] )i + QR (xR (ξ[R] ), ξ[R] )

(Q0 )

t=1

so dass xt : Rs·t → Xt messbar,
At,0 (ξt )xt (ξ[t] ) + At,1 (ξt )xt−1 (ξ[t−1] ) = ht (ξt ),t = 1, . . . , R,t = 2, . . . , R.

Diese Formulierung hat den Vorteil, dass das Originalproblem (P) in ein Problem für die(i)ersten R Zeitstufen (Q0 ) und je ein Teilproblem (QR ) für jede Realisierung ξ[R] von ξ[R]zerfällt3 , vgl. Abbildung 10.1.
(i)
3 Der Zufallsvektor ξ[R] kann die Werte ξ[R] , i = 1, . . . , nR , annehmen.

(1)xR (ξ[R] )(QR )(1) (1)QR (x[R] , ξ[R] )

(P)

(Q0)

(2)
xR (ξ[R] )
(2)

(QR )

(2)

QR (x[R] , ξ[R] )(3)xR (ξ[R] )(QR )(3) (3)QR (x[R] , ξ[R] )(4)xR (ξ[R] )
(4)

(QR )

(4)

QR (x[R] , ξ[R] )(5)xR (ξ[R] )
(5)

(QR )

(5)

QR (x[R] , ξ[R] )(6)xR (ξ[R] )
(6)

(QR )

(6)

QR (x[R] , ξ[R] )(7)xR (ξ[R] )
(7)

(QR )

(7)

QR (x[R] , ξ[R] )Abbildung 10.1: Zeitliche Dekomposition eines stochastischen Programmes

10.2.3

Benders Dekomposition

Innerhalb des Benders Dekompositionsverfahrens [1, 18, 22] werden nun die unbekannten stückweise linearen konvexen Funktionen(i)QR (·, ξ[R] ) : Rr → R ∪ {+∞},

i = 1, . . . , nR ,in Problem (Q0 ) sukzessive durch unterschätzende lineare Approximationen in Formvon Schnittebenen ersetzt. Die so erhaltene Approximation des Problems (Q0 ) ist ein
lineares Optimierungsproblem und wird als Masterproblem bezeichnet. Die Lösung des(i)

Masterproblems liefert Entscheidungen xR , i = 1, . . . , nR in den Endknoten des Masterproblems. Löst man nun die linearen Teilprobleme (QR ) in diesen Punkten, d.h. be(i)(i)

rechnet man die Funktionswerte QR (xR , ξ[R] ), so erhält man aus den dualen Lösungen(i)

Schnittebenen, die zur Verbesserung der Approximation von QR (·, ξ[R] ) und damit zur
Aktualisierung des Masterproblems genutzt werden können.
Das Masterproblem und die Probleme (QR ) werden nun abwechselnd gelöst, bis allenötigen Schnittebenen zur Approximation von (Q0 ) bestimmt wurden. Die Entscheidun(i)(i)
gen xR , in denen die Funktionen QR (·, ξ[R] ) ausgewertet werden, stammen aus einerLösung des Masterproblems. Durch dieses Vorgehen wird erreicht, dass die Funktionen(i)QR (·, ξ[R] ) vornehmlich dort gut approximiert werden, wo ihre Werte für die Lösung von(Q0 ) tatsächlich von Interesse sind. Für eine detaillierte Darstellung dieser Vorgehensweise sei auf [18, 22] verwiesen, die Modifikation dieses Algorithmus für rekombinierende Szenariobäume wird in Abschnitt 10.3 ausführlicher erklärt.
Für lange Zeithorizonte T können die Teilprobleme (QR ) weiter zerlegt und selbst mittels des oben beschriebenen Approximationsverfahrens gelöst werden. Dieses Vorgehenwird als Nested Benders Dekomposition [1] bezeichnet. Allerdings wächst mit der Anzahl der Knoten zum Zeitpunkt R auch die Anzahl der zu approximierenden Funktionen(i)QR (·, ξ[R] ), i = 1, . . . , nR .
10.2.4

Rekombinierende Szenariobäume

Eine besondere Situation ergibt sich, wenn die in zwei Knoten i und j zu einem Zeitpunkt
R ansetzenden Teilbäume übereinstimmen, d.h. die entsprechenden bedingten Wahrscheinlichkeitsverteilungen von (ξR+1 , . . . , ξT ) identisch sind:
h
i
h
i
(i)
(k)
P (ξR+1 , . . . , ξT ) = a ξ[R] = ξ[R] = P (ξR+1 , . . . , ξT ) = a ξ[R] = ξ[R](10.2)
für alle a ∈ Rs·(T −R) , vgl. linke Grafik in Abbildung 10.2. Ein Ansatz zur Konstruktion
von Szenariobäumen mit der Eigenschaft (10.2) wird in Abschnitt 10.5 vorgestellt.

12x3H Ξ iLx3H Ξ iLx3H Ξ jLx3H Ξ jL

3

4

1

5

2

3

4

5

Abbildung 10.2: Szenariobaum (schematisch) mit Eigenschaft (10.2), R = 3 und mR = 2, d.h. zwei verschiedene Teilbäume beginnen zur Zeit 3; die schwarzen und die grauen Bäume stimmen
jeweils überein.
(i)

(i)

Ohne die zeitkoppelnde Restriktion (10.1) ist die Kostenfunktion QR (xR , ξ[R] ) unab(i)hängig vom Wert xR . Dann genügt es, innerhalb der Benders Dekomposition nur eines
(i)

(i)

(k)

(k)

der beiden identischen Teilprobleme QR (xR , ξ[R] ) und QR (xR , ξ[R] ) zu lösen4 . Insbesondere kann man dann den Prozess ξ durch einen Szenariobaum darstellen, in welchem
die Knoten i und k zusammengelegt, d.h. rekombiniert werden, vgl. Abbildung 10.2.
Durch wiederholtes Rekombinieren in verschiedenen Zeitstufen ergäbe sich ein Szenariobaum, dessen Knotenzahl nunmehr linear statt exponentiell mit der Anzahl der Zeitstufen wächst.
4 Dies gilt auch, wenn die optimalen Entscheidungen x(i) und x(k) in den Knoten i und k übereinstimmen.
R
R

(i)

(k)

Da im Allgemeinen jedoch die optimalen Entscheidungen xR und xR unterschiedlich
ausfallen werden, ist die Verwendung eines derart rekombinierenden Baumes zur Lösung eines Optimierungsproblems mit zeitkoppelnden Restriktionen nicht ohne Weiteresmöglich. Allerdings kann Eigenschaft (10.2), d.h. die Übereinstimmung von Teilbäumen,
(i)

(k)

von Nutzen sein, da diese die Identität der Kostenfunktionen QR (·, ξ[R] ) und QR (·, ξ[R] )impliziert. Einerseits ermöglicht dies die simultane Approximation5 der Kostenfunktionen, was bereits zu einer verbesserten Konvergenz des Lösungsalgorithmus führt.
Andererseits werden bei einer großen Zahl von Knoten innerhalb der Benders Dekomposition sehr viele Teilprobleme (QR ) gelöst, die sich nur in den auf der rechten Seite auf(i)tretenden Größen xR unterscheiden, vgl. (10.1). Da die Optimalwerte der Teilprobleme
stetig von der rechten Seite abhängen, genügt es, eine Auswahl weniger “repräsentativer”
rechter Seiten zu berücksichtigen, d.h. erfüllen i und k die Eigenschaft (10.2) und hat
(i)

(k)

(i)

(i)

man ähnliche Werte xR und xR , so löst man nur eines der Teilprobleme QR (xR , ξ[R] )
(k)

(k)

und QR (xR , ξ[R] ). Dieses Vorgehen wird im nächsten Abschnitt genauer beschrieben
und kann als dynamisches Rekombinieren von Szenarien interpretiert werden.

10.3

Einfache Dekomposition mit rekombinierenden Bäumen

In diesem Abschnitt wird zunächst die Situation betrachtet, dass der Szenariobaum die
Eigenschaft (10.2) in einem Zeitpunkt R ∈ {1, . . . , T − 1} besitzt. Es wird ein Lösungsalgorithmus vorgestellt, der diese Eigenschaft ausnutzt und das Originalproblem im Zeitpunkt R zerlegt. Im nächsten Abschnitt wird dieses Verfahren zu einer mehrfachen Dekomposition erweitert.
Die Benders Dekomposition betrachtet die dynamische Formulierung (Q0 ) des Problems(i)(P) und approximiert die (konvexen und stückweise linearen) Funktionen QR (·, ξ[R] ),i = 1, . . . , nR , von unten durch lineare Funktionen, vgl. Abbildung 10.3. Diese Unterschätzer werden durch das Lösen der linearen Optimierungsprobleme (QR ) bezüglicheiner adaptiv gewählten Folge von Punkten x̄R ∈ Rr bestimmt. Ist dabei das Problem(QR ) für einen solchen Punkt zulässig, so lässt sich aus dessen Lösung die Ungleichung
(i)

(i)

(i)

q̄ + hπ̄, xR − x̄R i ≤ QR (xR , ξ[R] )(i)für alle xR ∈ Rr

(10.3)(i)

ableiten. Dabei ist q̄ := QR (x̄R , ξ[R] ) der Optimalwert des Problems (QR ) undhi(i)π̄ := E AR+1,1 (ξR+1 )> πR+1 (ξR+1 ) ξ[R] = ξ[R] ,wobei πR+1 (ξR+1 ) der Teil der dualen Lösung von (QR ) ist, welcher zu der zeitkoppelnden Restriktion (10.1) gehört. Die linke Seite der Ungleichung (10.3) stellt eine lineare
5 Das bedeutet, dass die Schnittebenen, die zur Approximation jeweils einer der Kostenfunktionen verwendet werden, ebenfalls zur Approximation der anderen dienen.

(i)Funktion in xR dar und wird als Optimalitätsschnitt bezeichnet, da sie den Optimalwert
(i)

(i)

QR (xR , ξ[R] ) von (QR ) (nach unten) approximiert.
Falls das Problem (QR ) für den gegebenen Referenzpunkt x̄R unzulässig ist, so lässt(i)sich eine Ungleichung hd, xR i ≤ e bestimmen, welche für x̄R verletzt ist, aber für(i)jene xR erfüllt ist, für die (QR ) eine zulässige Lösung besitzt, vgl. [22]. Diese Art vonUngleichung wird als Zulässigkeitsschnitt bezeichnet.

Copt

(i)

QR (·, ξ[R])

Copt

Cfeas

(i)

QLR (·, ξ[R])

xR

(i)

Abbildung 10.3: Approximation einer Funktion QR (·, ξR ) durch Schnittebenen(i)

Die Kostenfunktion QR (·, ξ[R] ) ist durch die bedingte Wahrscheinlichkeitsverteilung des(i)

Zufallsvektors (ξR+1 , . . . , ξT ) unter der Bedingung {ξ[R] = ξ[R] } bestimmt. Diese Verteilung ist wiederum durch den im Knoten i zur Zeit R ansetzenden Teilbaum des Szenariobaums festgelegt. Ist der den Prozess ξ beschreibende Szenariobaum nun derart,
dass einige der in Zeitstufe R ansetzenden Teilbäume übereinstimmen, also die Bedingung (10.2) für zwei Knoten i, k ∈ {1, . . . , nR } erfüllt ist, so stimmen auch die
(i)

(k)

Funktionen QR (·, ξ[R] ) und QR (·, ξ[R] ) überein und können simultan approximiert wer(i)den. Dies bedeutet, dass jede die Funktion QR (·, ξ[R] ) approximierende Schnittebene(k)ebenfalls zur Approximation von QR (·, ξ[R] ) verwendet werden kann und umgekehrt,
vgl. Abbildung 10.4. Um dies zu formalisieren sei eine (nicht-surjektive) Abbildung
λ : {1, . . . , nR } → {1, . . . , nR } definiert, die jedem Knoten i zur Zeit R einen (beliebigen) repräsentativen Knoten k mit dem gleichen Teilbaum zuordnet:
λ(i) := kfür alle i ∈ {1, . . . , nR }, die mit k Bedingung (10.2) erfüllen. (10.4)

Der Wert des Prozesses ξ in dem repräsentativen Knoten λ(i) zur Zeit R wird im Fol(i)(λ(i))
genden mit λ(ξ[R] ) := ξ[R] bezeichnet.
(i)

Für gegebene Mengen von Optimalitätsschnitten Copt (λ(ξ[R] )) und Zulässigkeitsschnit(i)(i)

ten Cfeas (λ(ξ[R] )) für die Funktion QR (·, ξ[R] ) definiert man das folgende Masterpro-

(1)xR (ξ[R] )(QL0 )(π̄, q̄)(d, e)(3)(QR )xR (ξ[R] )

(QR )

[5]

xR (ξ[R] )(QR )Abbildung 10.4: Zeitliche Dekomposition mit rekombinierendem Szenariobaum

blem6 , vgl. [1, 18, 22]:
min E

x,wR

" R
X

#
(QL
0)

hbt (ξt ), xt (ξ[t] )i + wR (ξ[R] )t=1

so dass xt : Rs·t → Xt messbar,wR : Rs·R → R messbar,

At,0 (ξt )xt (ξ[t] ) + At,1 (ξt )xt−1 (ξ[t−1] ) = ht (ξt ),t = 1, . . . , R,t = 2, . . . , R,

hd, xR (ξ[R] )i ≤ e,(d, e) ∈ Cfeas (λ(ξ[R] )),q̄ + hπ̄, xR (ξ[R] ) − x̄R i − wR (ξ[R] ) ≤ 0,(x̄R , q̄, π̄) ∈ Copt (λ(ξ[R] )).

(i)

Da die (unbekannten) Funktionen QR (·, ξ[R] ) hier durch lineare Unterschätzerfunktionen
ersetzt wurden, liefert die Lösung von Problem (QL0 ) eine untere Schranke an den Optimalwert von Problem (P).
Wie in Abschnitt 10.2.3 beschrieben, passt der in Abbildung 10.5 dargestellte Lösungsalgorithmus durch Einfügen von Schnittebenen in Copt (·) und Cfeas (·) das Masterproblem (QL0 ) adaptiv an das Originalproblem (Q0 ) an. Insbesondere werden dabei beim
(k)

(i)

(k)

(k)

Lösen von Teilproblemen QR (xR , ξ[R] ) und QR (xR , ξ[R] ) Schnittebenen in dieselben
(i)

(i)

Mengen Copt (λ(ξ[R] )) und Cfeas (λ(ξ[R] )) eingefügt. Durch dieses simultane Aufdatieren kann bereits eine deutliche Beschleunigung des klassischen Benders Algorithmus
erreicht werden.
6 Die Schreibweise hd, x (ξR [R] )i ≤ e für (d, e) ∈ Cfeas (λ(ξ[R] )) in den Nebenbedingungen vonL(Q0 ) bedeutet, dass diese Bedingung für alle möglichen Realisierungen von ξ[R] gestellt wird, d.h. es wird(i)(i)
hd, xR (ξ[R] )i ≤ e für (d, e) ∈ Cfeas (λ(ξ[R] )) und i = 1, . . . , nR gefordert. Analog gilt dies für die Opti-

malitätsschnitte.

Falls die Anzahl der Knoten nR sehr groß ist, so muss auch eine sehr große Zahl verschie(i) (i)dener Teilprobleme QR (xR , ξ[R] ), i = 1, . . . , nR , gelöst werden. Eigenschaft (10.2)(i)führt nun aber dazu, dass nur noch wenige verschiedene Teilprobleme QR (·, λ(ξ[R] ))(i)

jeweils in einer sehr großen Anzahl von Punkten xR betrachtet werden müssen. Da nun
(i)
(i)
(i)
der Wert QR (xR , λ(ξ[R] )) stetig von xR abhängt, ist es naheliegend, sich auf die Be(i)rechnung von QR (·, λ(ξ[R] )) in einer reduzierten Auswahl von Punkten zu beschränken.

Diese Selektion wird in Schritt 4 des in Abbildung 10.5 dargestellten Algorithmus durchgeführt; die Genauigkeit dieser Operation wird durch den Parameter ρ gesteuert und im
Verlauf der Optimierung sukzessive erhöht. Somit werden die Teilprobleme in Knoten iund k während einer Iteration des Algorithmus’ identifiziert, falls sich die zeitkoppelnden(i)(k)
Variablen xR und xR in diesem Iterationsschritt um weniger als den Wert ρ unterscheiden. Dieses Vorgehen kann als dynamisches Rekombinieren von Szenarien interpretiert
werden.
1. Setze Cfeas (·) := ∅, Copt (·) := {(0, L, 0)} und ρ := ρstart , wobei die Konstante
L eine untere Schranke an QR (·, ·) ist.
2. Löse das Masterproblem (QL
0 ).
3. Falls (QL0 ) unzulässig ist, so stoppe: (Q0 ) ist unzulässig.
Andernfalls sei (x̄, w̄) eine Lösung von (QL
0 ).
(i)

(i)

(i)

(i)

4. Für alle ξ[R] , i = 1, . . . , nR : Falls kein xR mit kxR − x̄R (ξ[R] )k ≤ ρ existiert,
(i)

(i)

für welches QR (xR , λ(ξ[R] )) bereits bestimmt wurde:
(i)

(i)

(a) Bestimme QR (x̄R (ξ[R] ), λ(ξ[R] )) durch Lösung von (QR ).
(b) Falls (QR ) unzulässig ist, so konstruiere einen Zulässigkeitsschnitt und(i)füge diesen zu Cfeas (λ(ξ[R] )) hinzu.
(i)

(i)

(i)

(c) Falls (QR ) zulässig ist und QR (x̄R (ξ[R] ), λ(ξ[R] )) > w̄R (ξ[R] ), so kon(i)struiere einen Optimalitätsschnitt und füge diesen zu Copt (λ(ξ[R] )) hinzu.
5. Falls Optimalitäts- oder Zulässigkeitsschnitte generiert wurden, gehe zurück zu
2.
Sonst: Falls ρ > ρend , verringere ρ und gehe zurück zu 2.
Sonst: Falls ρ ≤ ρend , stoppe. Problem (Q0 ) wurde (ρend -) optimal gelöst.
Abbildung 10.5: Lösung von (Q0 ) mittels einer modifizierten Benders Dekomposition.

Der Lösungsalgorithmus stoppt wenn das Masterproblem (QL0 ) unzulässig ist oder die
Generierung von Schnittebenen mit minimalem ρ = ρend beendet ist. Im ersten Fall

bedeutet dies die Unzulässigkeit des Originalproblems (Q0 ), während im zweiten Falleine (ρend -) optimale Lösung von (Q0 ) bestimmt wurde, vgl. Proposition 10.4.1. Andere
Stoppkriterien, welche es erlauben den Algorithmus mit ρ > ρend zu beenden sobald derApproximationsfehler unter eine vorgegebene Toleranz fällt, werden in Abschnitt 10.4.2
diskutiert.

10.4

Mehrfache Dekomposition mit rekombinierenden Bäumen

Besitzt der den Prozess ξ beschreibende Szenariobaum die Eigenschaft (10.2) in mehreren Zeitpunkten R1 , . . . , Rn mit
0 =: R0 < R1 < . . . < Rn < Rn+1 := T,
so lässt sich der im letzten Abschnitt vorgestellte Dekompositionsalgorithmus rekursiv
erweitern.
Wie in Abschnitt 10.2.2 definiert man dazu Kostenfunktionen durch QRn+1 ≡ 0 und
(i)

(i)

(QRj )QRj (xRj , ξ[Rj ] ) :=


Rj+1
X
(i)
min E hbt (ξt ), xt (ξ[t] )i + QRj+1 (xRj+1 (ξ[Rj+1 ] ), ξ[Rj+1 ] ) ξ[Rj ] = ξ[Rj ] 
t=Rj +1

so dass xt : Rs·t → Xt messbar,t = Rj + 1, . . . , Rj+1 ,(i)

ARj +1,0 (ξRj +1 )xRj +1 (ξ[Rj +1] ) = hRj +1 (ξRj +1 ) − ARj +1,1 (ξRj +1 )xRj
At,0 (ξt )xt (ξ[t] ) + At,1 (ξt )xt−1 (ξ[t−1] ) = ht (ξt ),(i)

t = Rj + 2, . . . , Rj+1 ,(k)für j = 0, . . . , n. Wieder stimmen die Funktionen QRj (·, ξ[Rj ] ) und QRj (·, ξ[Rj ] ) überein, falls die Knoten i und k zur Zeitstufe Rj die Bedingung (10.2) mit R = Rj erfüllen.
Wie in Definition (10.4) seien für j = 1, . . . , n Abbildungen λj definiert, die jedem Knoten i im Zeitpunkt Rj einen repräsentativen Knoten k mit dem gleichen Teilbaum zuordnen. Es sei mj die Anzahl verschiedener Teilbäume, die im Zeitpunkt Rj ansetzen. Für(λ (i))(ī)
ī = 1, . . . , mj bezeichnet λj := ξ[Rjj] den Wert von ξ bis zum ī-ten repräsentativen
Knoten7 zum Zeitpunkt Rj .

7 Der Wert i in der Definition von λ(ī) bezeichnet einen beliebigen Knoten zum Zeitpunkt R , der ebenfallsjj

den ī-ten Teilbaum besitzt, d.h. λj (i) = ī.

10.4.1Lösung mit mehrfach rekombinierenden Bäumen

Analog zu (QL0 ) seien für j = 0, . . . , n und i = 1, . . . , mj die Masterprobleme
(i)

(i)

QL
Rj (xRj , ξ[Rj ] ) :=


Rj+1
X
(i)
min E hbt (ξt ), xt (ξ[t] )i + wRj+1 (ξ[Rj+1 ] ) ξ[Rj ] = ξ[Rj ] x,wRj+1(QL
Rj )t=Rj +1

xt : Rs·t → Xt , wRj+1 : Rs·Rj+1 → R messbar,

so dasst = Rj + 1, . . . , Rj+1 ,(i)

ARj +1,0 (ξRj +1 )xRj +1 (ξ[Rj +1] ) = hRj +1 (ξRj +1 ) − ARj +1,1 (ξRj +1 )xRj ,At,0 (ξt )xt (ξ[t] ) + At,1 (ξt )xt−1 (ξ[t−1] ) = ht (ξt ),hd, xRj+1 (ξ[Rj+1 ] )i ≤ e,t = Rj + 1, . . . , Rj+1 ,(d, e) ∈ Cfeas (λj+1 (ξ[Rj+1 ] )),q̄ + hπ̄, xRj+1 (ξ[Rj+1 ] ) − x̄Rj+1 i − wRj+1 (ξ[Rj+1 ] ) ≤ 0,(x̄Rj+1 , q̄, π̄) ∈ Copt (λj+1 (ξ[Rj+1 ] )),
(ī)

(ī)

definiert, wobei Copt (λj+1 ) und Cfeas (λj+1 ), ī = 1, . . . , mj+1 , Mengen von Optima(ī)

litäts- bzw. Zulässigkeitsschnitten sind, welche die Funktion QL
Rj+1 (·, λj+1 ) approximieren8 , siehe auch [15].
Der Lösungsalgorithmus löst die Masterprobleme (QLRj ), j = 0, . . . , n, abwechselnd ineinem Vorwärts- und einem Rückwärtsmodus, vgl. auch Abbildung 10.6. Der Vorwärtsmodus dient der Weitergabe primaler Information in Form von Kontrollpunkten an zukünftige Zeitstufen. Dabei werden nacheinander die Masterprobleme in den Zeitstufen
R0 , . . . , Rn gelöst. Dabei liefert das Lösen der Probleme zum Zeitpunkt Rj−1 Mengen
(ī)
(ī)
(ī)
(ī)
Zj von Paaren (xRj , wRj ), wobei wRj der Wert der Schnittebenenapproximation von
(ī)

(ī)

QL
Rj (·, λj ) in der Entscheidungsvariablen xRj ist. Im Zeitpunkt Rj werden die Funk(ī)

(ī)

(ī)

tionen QL
Rj (·, λj ), ī = 1, . . . , mj , für die Werte xRj aus Zj berechnet, d.h. die Pro(ī)

(ī)

(ī)

L
bleme (QL
Rj ) gelöst. Falls nun wRj < QRj (xRj , λj ), so ist die Approximation an der(ī)

Stelle xRj noch nicht exakt und ein Optimalitäts- oder Zulässigkeitsschnitt wird aus derLLösung von (QLRj ) generiert, vgl. Abschnitt 10.3 bzw. [22]. Falls (QRj ) zulässig war, so

werden die Lösungen xRj+1 (·) zusammen mit den Werten wRj+1 (·) aus der Lösung von(k̄)(QLRj ) den Mengen Zj+1 , k̄ = 1, . . . , mj+1 , für den nächsten Zeitschritt Rj+1 zugefügt.
Falls ein Teilproblem unzulässig war, d.h. der zugehörige Kontrollpunkt durch einen
Zulässigkeitsschnitt in der vorhergehenden Zeitperiode abgeschnitten werden muss, so
wechselt der Algorithmus kurzfristig in den Rückwärtsmodus.
Der Rückwärtsmodus dient der Rückgabe dualer Information in Form von Optimalitätsund Zulässigkeitsschnitten an zurückliegende Zeitstufen. Dabei werden nacheinander die
Masterprobleme (QLRj ) der Zeitstufen Rn , . . . , R0 in den im Vorwärtsmodus erzeugten
8 Für j = n ist QL
L
Rj+1 (·, ·) ≡ 0, d.h. es sei Cfeas (·) := ∅ und Copt (·) := {(0, 0, 0)} in (QRn ).

(ī)

(ī)

Punkten xRj ∈ Zj gelöst. Die Lösung eines Masterproblems zur Zeit Rj aktualisiert
dabei gegebenenfalls die im Vorwärtsmodus bestimmten Schnittebenen. Die aktualisier(ī)(ī)
ten Schnitte werden zu Copt (λj ) bzw. Cfeas (λj ) hinzugefügt, d.h. die Masterprobleme
zur Zeit Rj−1 werden aufdatiert.
(QL0 )

(1)

(1)

Z1

Z2
(1)

QLR1 (·, λ1 )(1)

QLR2 (·, λ2 )

(2)

Z2

(2)

Z1

(2)

QLR2 (·, λ2 )(2)

QLR1 (·, λ1 )

(3)

(3)

Z1

QLR1 (·, λ1 )

(3)

Z2

(3)

QLR2 (·, λ2 )Abbildung 10.6: Schematische Darstellung des Dekompositionsalgorithmus für mehrfach rekombinierende
Szenariobäume

Der Algorithmus ist beendet wenn entweder das Masterproblem (QL0 ) in der ersten Zeitstufe unzulässig ist, oder alle Masterprobleme zur Optimalität gelöst werden konnten
und keine Schnittebenen mehr generiert werden. Im ersten Fall ist auch das Originalproblem (P) unzulässig, im zweiten Fall wurde es zur Optimalität gelöst. Eine detaillierte
Beschreibung des Nested Benders Algorithmus findet man in den Arbeiten [1, 9, 18].
Im allgemeinen Fall entspricht die Anzahl der Teilprobleme zur Zeit Rj der Anzahl der
Knoten nRj des Baumes zum Zeitpunkt Rj . Im betrachteten Fall übereinstimmender
Teilbäume gibt es nur mj verschiedene Teilprobleme, und damit auch nur mj zu approximierende Kostenfunktionen. Daher kann so eine deutlich schnellere Konvergenz desLösungsalgorithmus erreicht werden. Allerdings muss jedes der mj verschiedenen Teil(ī)(ī)

probleme wiederum für alle Punkte xRj ∈ Zj

gelöst werden. Wie in Abschnitt 10.3

(ī)
(ī)
erlaubt es die Stetigkeit der Funktionen QL
Rj (·, λj ) wieder, die Mengen Zj derart zu

reduzieren, dass die verbleibenden Punkte paarweise jeweils einen Abstand größer als
ein Parameter ρ zueinander haben.
Zu Beginn des Lösungsprozesses sind noch keine Informationen (in Form von Schnittebenen) über die Kostenfunktionen in den Masterproblemen bekannt. Daher können hier(ī)

Punkte xRj generiert werden, die noch weit von einer optimalen Lösung entfernt liegen.
Wählt man zu Beginn nun einen großen Wert für ρ, so werden nur wenige dieser Punkte
ausgewertet und die Kostenfunktionen fernab der optimalen Lösung nicht unnötig genau
approximiert. Hat man den Baum einige Male mit hohem Wert von ρ durchlaufen, so

liefern die Masterprobleme bereits eine grobe Darstellung der Kostenfunktionen und die
im Folgenden generierten Kontrollpunkte befinden sich bereits in der Nähe der optimalenLösung. Reduziert man nun den Parameter ρ sukzessive, so werden die Approximationen
der Kostenfunktionen in der Nähe der optimalen Lösung weiter verfeinert. Die Wahl desletzten, minimalen ρ-Wertes, für welchen Kostenfunktionen ausgewertet und Schnittebenen generiert werden, bestimmt damit auch die Güte der approximierten Lösung:Proposition 10.4.1 ([15]). Es sei v der Optimalwert von (P) und v L der Optimalwert(ī)

9der Approximation (QL0 ). Ist (QRj ) lösbar für alle xRj ∈ XRj , so gilt die folgende

Aussage: Der Nested Benders Dekompositions-Algorithmus mit Reduktion der Mengen(ī)
Zj , ī = 1, . . . , mj , stoppt nach endlich vielen Schritten mit einer Lösung von (QL0 ) dieρ-optimal für (P) ist, d.h. es existiert eine Konstante C ≥ 0, so dass|v − v L | ≤ Cρ.

10.4.2(10.5)

Ein adaptives Abbruchkriterium durch Verwendung oberer Schranken

Die Wahl des Parameters ρ erlaubt eine Steuerung der Approximationsgüte und hat einen
erheblichen Einfluss auf die Laufzeit des Lösungsalgorithmus [15]. Einerseits zeigt Proposition 10.4.1, dass für hinreichend kleine Werte von ρ die Genauigkeit der Approximation beliebig gut wird. Andererseits ist der Wert der Größe C in (10.5) im Allgemeinennicht bekannt und der tatsächliche Approximationsfehler mitunter deutlich kleiner als
Cρ. Daher lässt sich aus der Abschätzung (10.5) im Allgemeinen nicht ableiten, wie ρ
gewählt werden sollte um eine festgelegte Approximationsgüte zu garantieren.
(ī)

Im Folgenden wird eine Möglichkeit studiert den Approximationsfehler QRj (·, λj ) −
(ī)

L
QL
Rj (·, λj ) und damit den Term |v − v | während des Lösens abzuschätzen. Diese(ī)

besteht darin, zusätzlich zu den unteren Approximationen QL
Rj (·, λj ) auch obere Ap(ī)

proximationen QU
Rj (·, λj ) zu bestimmen. Solche oberen Schranken lassen sich rekursiv(ī)

wie folgt bestimmen. Es sei Yj+1 (λj+1 ) eine Menge von Punkten y ∈ XRj+1 und Werten(ī)

q̄yU , von denen bekannt ist, dass sie den Wert QRj+1 (y, λj+1 ) überschätzen. Aufgrund der(ī)
Konvexität von QRj+1 (·, λj+1 ) lässt sich dann mittels Interpolationen der Werte q̄yU eine(ī)
Funktion konstruieren, die QRj+1 (·, λj+1 ) nach oben abschätzt, vgl. Abbildung 10.7. Er(ī)
setzt man im Masterproblem (QRj ) die Kostenfunktion QRj+1 (·, λj+1 ) durch eine solche(ī)
interpolierende Funktion, so erhält man wiederum eine obere Schranke an QRj (·, λj ).
(ī)
(ī)
Formal sei dazu zunächst QU
Rn (·, λn ) := QRn (·, λn ), ī = 1, . . . , mn , definiert. Für
9 Diese Voraussetzung dient hier zur Vereinfachung und ist bei geeigneter Reduzierung der Mengen Z (ī)
jnicht notwendig, vgl. [15].

j = 0, . . . , n − 1 sei weiterhin
(ī)

(ī)

QU
Rj (xRj , λj ) :=
Rj+1
X
min E hbt (ξt ), xt (ξ[t] )i +x,αso dass

t=Rj +1

(QU
Rj )

X

(ī)

αy (ξ[Rj+1 ] ) q̄yU ξ[Rj ] = λj y∈Yj+1 (λj+1 (ξ[Rj+1 ] ))xt : Rs·t → Xt messbar,t = Rj + 1, . . . , Rj+1 ,(i)

ARj +1,0 (ξRj +1 )xRj +1 (ξ[Rj +1] ) = hRj +1 (ξRj +1 ) − ARj +1,1 (ξRj +1 )xRj ,At,0 (ξt )xt (ξ[t] ) + At,1 (ξt )xt−1 (ξ[t−1] ) = ht (ξt ),αy : Rs·Rj+1 → R+ messbar,
X
xRj+1 (ξ[Rj+1 ] ) =t = Rj + 2, . . . , Rj+1 ,y ∈ Yj+1 (λj+1 (ξ[Rj+1 ] )),αy (ξ[Rj+1 ] ) y,y∈Yj+1 (λj+1 (ξ[Rj+1 ] ))Xαy (ξ[Rj+1 ] ) = 1,y∈Yj+1 (λj+1 (ξ[Rj+1 ] ))(ī)

wobei die Mengen Yj+1 (λj+1 ), ī = 1, ..., mj+1 , aus Kontrollpunkten (oder Stützstellen)
(ī)

(ī)

y ∈ Yj+1 (λj ) bestehen, für die der Wert q̄yU := QU
Rj+1 (y, λj+1 ), bereits berechnet wurde. Die Wahl der αy (ξ[Rj+1 ] ) zwingen die Kontrollpunkte xRj+1 (ξ[Rj+1 ] ) einen Punkt inder konvexen Hülle von Yj+1 (λj+1 (ξ[Rj+1 ] )) anzunehmen. Die Kosten dieses Punktes
im Sinne der Kostenfunktion QRj+1 (·, λj+1 (ξ[Rj+1 ] )) werden dann aus der zugehörigen
Konvexkombination der bekannten Werte q̄yU nach oben abgeschätzt.

(ī)

QURj (·, λj )
ŷ

ŷ

ŷ

Abbildung 10.7: Obere Schranken (gestrichelt) an Kostenfunktionen mittels Konvexkombination aus bereit
bestimmten Werten in Punkten ŷ.
(ī)

Die Stützstellen Yj (λj ), ī = 1, . . . , mj , werden analog zu den Optimalitäts- und Zu(ī)

lässigkeitsschnitten gebildet. Zu Beginn des Lösungsprozesses sind die Mengen Yj (λj )leer, wodurch (QURj ) für j < n den Wert +∞ annimmt. Während der Algorithmus den
Baum im Vorwärts- und Rückwärtsmodus durchläuft, werden die Probleme (QURj ) für
(ī)

(ī)

verschiedene xRj ∈ Zj
(ī)

(ī)

gelöst. Liefert (QURj ) eine zulässig Lösung, so wird xRj zu

Yj (λj ) hinzugefügt und damit die Approximationen zur Zeitstufe Rj−1 verbessert.

Die Approximationen (QURj ) liefern sowohl ein “lokales” (d.h. auf einzelne Knoten bezogenes) als auch ein “globales” (d.h. auf den gesamten Baum bezogenes) Abbruchkrite(ī)

(ī)

(ī)

rium: Wurde für einen Kontrollpunkt xRj sowohl die untere Schranke QL
Rj (xRj , λj ) als
(ī)

(ī)

auch die obere Schranke QU
Rj (xRj , λj ) bestimmt und liegen diese nahe genug beieinander, so müssen die Kontrollpunkte xRj+1 (·) aus der Lösung von (QLRj ) nicht mehr an
die nächste Zeitperiode weitergegeben werden, da der zu approximierende (unbekannte)
(ī)
(ī)
Wert QRj (xRj , λj ) zwischen der unteren und oberen Schranke liegen muß. Bezeichnetv U den Optimalwert von (QU
R0 ), so erhält man die Fehlerabschätzungv − vL ≤ vU − vL .
Dies bedeutet, dass der Lösungsalgorithmus beendet werden kann, sobald der Wert v U −v L eine gegebene Toleranz unterschreitet. Abhängig vom Problem kann dies bereits für
relativ große Werte von ρ erreicht werden, vgl. [15].

10.5

Konstruktion rekombinierender Szenariobäume

In diesem Abschnitt wird eine Methode zur Konstruktion von Szenariobäumen vorgestellt, welche die Eigenschaft (10.2) für gewisse Knoten zu verschiedenen Zeitstufen
R1 , . . . , Rn besitzen. Ausgehend von einer Menge von Szenarien, die beispielsweise aus
historischen Daten oder einer Simulation stammen können, wird durch stufenweises Zusammenfassen der Szenarien zu Knoten der Baum konstruiert. In diesem Sinne ist dievorgestellte Methode eng verwandt mit dem in [12] entwickelten Verfahren.
0

Die im Zeitpunkt t = 1 übereinstimmenden N Ausgangs-Szenarien seien mit ζ i =00(ζ1i , . . . , ζTi ), i0 = 1, . . . , N , bezeichnet. In dem in Abbildung 10.8 beschriebenen Algorithmus werden einerseits für jede Zeitperiode [Rj + 1, Rj+1 ] nicht-rekombinierende
Szenariobäume konstruiert (Schritt 2), andererseits werden zum Zeitpunkt Rj+1 Knoten
verschiedener Teilbäume mit ähnlichen Werten identifiziert. Diese Knoten erhalten dannden gleichen Teilbaum für die folgende Zeitperiode [Rj+1 + 1, Rj+2 ] (Schritt 1). Die
verwendete Notation wird in Tabelle 10.1 erklärt. Die Werte mj , nj (h), und st+1 (i),welche die Struktur des Szenariobaumes bestimmen, können entweder a priori festgelegt, oder innerhalb des Algorithmus so bestimmt werden, dass gewisse Fehlerschranken
eingehalten werden.
Die Bestimmung der Mengen A in den Schritten 1a und 2a,b sind sogenannte k−meanclustering Probleme, welche man mit Hilfe der in [4, 11] vorgeschlagenen forward selection Heuristik approximativ lösen kann.

(1,1)

Initialisierung: Setze C1

(1)

(1)

:= C R0 := {1, . . . , N } , ξ1

(h,i)

(h)

:= ζ11 und C Rj := Ct:= ∅für alle h, j, i.
Für j = 0, . . . , n (Konstruktion für den Zeitraum [Rj + 1, Rj+1 ].)1. Falls j > 0: (Rekombinieren.)(a) Finde eine Menge A = {a1 , . . . , amj } ⊂ {1, . . . , N }, so dass der folgende
Abstand aller Knoten (h, i) zur Zeit Rj zum nächsten Element ζ al klein ist:X
h=1,...,mj−1
i=1,...,nRj (h)

0

X

min

al ∈A

0

a

a

i
i
k(ζR
, . . . , ζR
) − (ζRlj −τ , . . . , ζRlj )k.
j −τ
j

(h,i)
j

i0 ∈CR(b) Für alle Knoten (h, i) zur Zeit Rj : (h = 1, . . . , mj−1 , i = 1, . . . , nRj (h))
i. Bestimme ein ah0 ∈ A so dass Szenario ζ ah0 nahe an Knoten (h, i) liegt:
ah0 ∈ arg min
al ∈A

0

X

0

a

a

i
i
k(ζR
, . . . , ζR
) − (ζRlj −τ , . . . , ζRlj )k.
j −τ
j

(h,i)
i0 ∈CR
jii. Knoten (h, i) zur Zeit Rj erhält Teilbaum h0 für den Zeitraum [Rj +1, Rj+1 ]; die Ausgangs-Szenarien die durch Knoten (h, i) gehen, werden
zur Konstruktion des Teilbaums h0 genutzt. Setze dazu
(h0 )

(h0 )

(h,i)

C Rj := C Rj ∪ CRj .
2. Für h = 1, . . . , mj : (Konstruktion der Teilbäume für [Rj + 1, Rj+1 ])(h)(a) Finde eine Menge A = {a1 , . . . , anRj +1 (h) } ⊂ C Rj so dass der folgende TermXai0klein wird:
min kζR− ζRlj +1 k.
j +1
(h)
j

al ∈A

i0 ∈C R

0

(h)

a

i
Für alle i0 ∈ C Rj : Wähle ai ∈ arg minal ∈A kζR− ζRlj +1 k und füge i0 zu
j +1
(h,i)
CRj +1 hinzu.
(h,i)

ai
.
Definiere den Wert von ξRj +1 auf dem Knoten (h, i) durch ξRj +1 := ζR
j +1

Für alle Knoten (h̃, ĩ) zur Zeit Rj (h̃ = 1, . . . , mj−1 , ĩ = 1, . . . , nRj (h̃)) mit(h̃,ĩ)CRj(h)⊂ C Rj : definiere die Übergangswahrscheinlichkeit von (h̃, ĩ) zu (h, i):(h,i)(h̃,ĩ)

(h,i)

(h)

P[ ξRj +1 = ξRj +1 | ξ[Rj ] = ξ[Rj ] ] := |CRj +1 | / |C Rj |.
(b) Für t = Rj + 1, . . . , Rj+1 − 1:Für alle Knoten ĩ des Teilbaumes h zur Zeit t:(h,ĩ)

Finde eine Menge A = {a1 , . . . , ast+1 (ĩ) } ⊂ Ct
so dass der folgende TermXali0klein wird:
min kζt+1 − ζt+1 k.
(h,ĩ)al ∈A

i0 ∈Ct(h,ĩ)

Für alle i0 ∈ Ct

0

a

i
l
: Wähle ai ∈ arg minal ∈A kζt+1− ζt+1
k und füge i0 zu(h,i)
Ct+1 hinzu.
(h,i)
(h,i)
(h̃,ĩ)(h,i)(h,ĩ)
ai
Setze ξt+1 := ζt+1
und P[ ξt+1 = ξt+1 | ξ[t] = ξ[t] ] := |Ct+1 | / |Ct |.

Abbildung 10.8: Konstruktion eines rekombinierenden Szenariobaums.

0

ζti

Wert von Ausgangs-Szenario i0 zur Zeit t(i)ξt

Wert der Zufallsvariablen ξt in Knoten imjAnzahl der verschiedenen Teilbäume die zur Zeit Rj beginnen (m0 = 1)

(h, i)

Knoten i des Teilbaumes hnt (h)

Anzahl der Knoten des Teilbaumes h zur Zeit Rj

st+1 (i)Anzahl der Knoten zur Zeit t + 1, die von Knoten i zur Zeit t ausgehenτLänge des zurückschauenden Zeithorizonts beim Rekombinieren
(h,i)

Ct

0

Teilmenge von {1, . . . , N }, bezeichnet die Ausgangs-Szenarien ζ idie zur Zeit t im Knoten (h, i) liegen

(h)
C Rj

0

Teilmenge von {1, . . . , N }, bezeichnet die Ausgangs-Szenarien ζ i
die in Teilbaum h mit Beginn zur Zeit Rj liegen

Tabelle 10.1: Bezeichnungen innerhalb des Baumkonstruktionsalgorithmus.

10.5.1

Stabilität

Mit dem oben vorgestellten Approximationsverfahren kann man aus einer Menge vonSzenarien ζ eines zugrunde liegenden stochastischen (Original-)Prozesses ζ einen rekombinierenden Szenariobaum konstruieren. Es stellt sich die Frage, unter welchen Bedingungen diese Art der Approximation konsistent ist, d.h. wann beispielsweise der mit
Hilfe eines solchen Szenariobaums ξ berechnete Optimalwert v(ξ) nahe am Optimalwertv(ζ) des Originalproblems liegt. In [14] wurde gezeigt, dass unter gewissen Voraussetzungen an den Ausgangsprozess ζ sowie die Struktur des Optimierungsproblems der
Optimalwert v(ξ) des approximierten Problems bei Verfeinerung der Diskretisierung in
der Tat gegen den Optimalwert v(ζ) des Originalproblems konvergiert. Neben gewissen technischen Voraussetzungen ist dies die folgende Voraussetzung an die Struktur desOriginalprozesses ζ:
Annahme. Der Prozess ζ ist durch die folgende Rekursion definiert:
ζt+1 := gt+1 ((ζt−τ , . . . , ζt ), εt+1 ) ,(10.6)

wobei (εt )t=1,...,T eine Folge unabhängiger Zufallsvariablen ist und gt+1 eine Abbildung von Rs·(τ +1) nach Rs , die die folgende Lipschitz-Bedingung erfüllt: Es gibt eine
Konstante L ≥ 0, so dass
t


X

gt+1 (ζt−τ , . . . , ζt ), εt+1 − gt+1 (ζ̂t−τ , . . . , ζ̂t ), εt+1 ≤ L
kζr − ζ̂r k
r=t−τ

(10.7)
für alle möglichen Werte (ζt−τ , . . . , ζt ), (ζ̂t−τ , . . . , ζ̂t ) des Zufallsvektors (ζt−τ , . . . , ζt )und alle Werte εt+1 von εt+1 . Bedingung (10.7) kann relaxiert werden, indem man einen
Lipschitz Faktor L betrachtet, der von den Werten ζ, ζ̂ und εt+1 abhängen darf.
Diese Bedingungen werden beispielsweise durch die in einer Vielzahl von Anwendungengenutzten autoregressiven Zeitreihenmodelle erfüllt.

10.6

Schlussbemerkungen

In diesem Kapitel wurde eine Modifikation der Nested Benders Dekomposition für mehrstufige stochastische Optimierungsprobleme vorgestellt, die auf dem Konzept rekombinierender Szenariobäume beruht. Der Dekompositionsansatz ermöglicht dabei das “dynamische Rekombinieren” von Teilbäumen in Abhängigkeit von den betrachteten Werten
der Entscheidungsvariablen. Dies erlaubt es, zeitkoppelnde Restriktionen in das Konzept
rekombinierender Szenariobäume zu integrieren. Weiterhin wurde ein auf [12] basierender Algorithmus zur Konstruktion von rekombinierenden Szenariobäumen vorgestellt.
Wie in [7] gezeigt, lässt sich das vorgestellte Dekompositionsverfahren unmittelbar auf
Probleme mit diskreten Entscheidungen in der ersten Stufe verallgemeinern. Diskrete
Entscheidungen nach dem ersten Rekombinierungszeitpunkt führen dagegen im Allgemeinen zum Verlust der Konvexität der zu approximierenden Kostenfunktionen, so dasssich eine wie in Abschnitt 10.2.3 diskutierte Schnittebenenapproximation nicht ohne
Weiteres bestimmen lässt. Eine Möglichkeit zur Berücksichtigung diskreter Entscheidungsvariablen besteht möglicherweise in einer Verallgemeinerung des in [20] vorgestellten Dekompositionsverfahren für zweistufige gemischt-binäre stochastische Optimierungsprobleme auf den mehrstufigen Fall.
Eine weitere mögliche Verallgemeinerung besteht in der Berücksichtigung von Risikomaßen in der Zielfunktion, siehe [5, 6, 17] sowie Kapitel 12. Wie in [6] gezeigt, erlaubt
die Verwendung mehrperiodischer polyhedraler Risikomaße wie dem Average-Value-atRisk weiterhin die Anwendung verschiedener Dekompositionsverfahren, u.a. auch der indiesem Kapitel verwendeten zeitlichen Dekomposition.
Danksagung.

Wir danken dem BMBF für die finanzielle Unterstützung. Weiterhindanken wir Prof. Werner Römisch für seine stetige Unterstützung und Dr. Holger Heitschfür hilfreiche Diskussionen.

10.7 Literatur[1]

B IRGE, J.R.: Decomposition and Partitioning Methods for Multistage Stochastic
Programming. In: Operations Research 33 (1985), Nr. 5, S. 989–1007

[2]C OX, J. ; ROSS, S. ; RUBINSTEIN, M.: Option Pricing: a simplified approach. In:
Journal on Financial Economics 7 (1979), S. 229–263[3]

D UPA ČOV Á, J. ; C ONSIGLI, G. ; WALLACE, S.W.: Scenarios for multistage stochastic programming. In: Annals of Operations Research 100 (2000), S. 25–53[4]

D UPA ČOV Á, J. ; G R ÖWE -K USKA, N. ; R ÖMISCH, W.: Scenarios reduction in stochastic programming: An approach using probability metrics. In: Mathematical
Programming 95(A) (2003), S. 493–511[5]

E ICHHORN, A.: Stochastic Programming Recourse Models: Approximation, Risk
Aversion, Applications in Energy, Humboldt-Universität zu Berlin, Diss., 2007. –Logos Verlag (Berlin)

[6]

[7]

E ICHHORN, A. ; R ÖMISCH, W.: Polyhedral risk measures in stochastic programming. In: SIAM Journal on Optimization 16 (2005), S. 69–95
E PE, A. ; K ÜCHLER, C. ; R ÖMISCH, W. ; V IGERSKE, S. ; WAGNER, H.-J. ; W E Stochastische Optimierung mit rekombinierenden Szena-

BER , C. ; W OLL , O.:

riobäumen – Analyse dezentraler Energieversorgung mit Windenergie und Speichern. In: Optimierung in der Energiewirtschaft. Düsseldorf : VDI-Verlag, 2007
(VDI-Berichte 2018), S. 3–13
[8]

F ELIX, B. ; W EBER, C.: Bewertung von unterschiedlichen Gasspeicherkonfigurationen mittels Verwendung von rekombinierenden Bäumen. In: Optimierung in derEnergiewirtschaft. Düsseldorf : VDI-Verlag, 2007 (VDI-Berichte 2018), S. 107–
120

[9]

G ASSMANN, H.I.: MSLiP: a computer code for the multistage stochastic linear
programming problem. In: Mathematical Programming 47 (1990), S. 407–423[10] G R ÖWE -K USKA, N. ; H EITSCH, H. ; R ÖMISCH, W.: Scenario Reduction and Scenario Tree Construction for Power Management Problems. In: B ORGHETTI, A.
(Hrsg.) ; N UCCI, C.A. (Hrsg.) ; PAOLONE, M. (Hrsg.): IEEE Bologna Power Tech
Proceedings, 2003[11] H EITSCH, H. ; R ÖMISCH, W.: Scenario Reduction Algorithms in Stochastic Programming. In: Computational Optimization and Applications 24 (2003), S. 187–
206[12] H EITSCH, H. ; R ÖMISCH, W.: Scenario tree modeling for multistage stochastic
programs. In: Mathematical Programming (2008). – noch nicht erschienen

[13] I NFANGER, G. ; M ORTON, D.P.: Cut sharing for multistage stochastic linear programs with interstage dependency. In: Mathematical Programming 75 (1996), S.
241–256[14] K ÜCHLER, C.: Stability, Approximation, and Decomposition in Two- and Multistage Stochastic Programming, Humboldt-Universität zu Berlin, Diss., 2008. – inVorbereitung[15] K ÜCHLER, C. ; V IGERSKE, S.: Decomposition of Multistage Stochastic Programswith Recombining Scenario Trees. In: Stochastic Programming E-Print Series 9(2007). – http://www.speps.org[16] P EREIRA, M.V.F. ; P INTO, L.M.V.G.: Multi-stage stochastic optimization applied
to energy planning. In: Mathematical Programming 52 (1991), S. 359–375[17] P FLUG, G.Ch. ; R ÖMISCH, W.: Modeling, Measuring and Managing Risk. Singapore : World Scientific, 2007[18] RUSZCZY ŃSKI, A.: Decomposition Methods. In: [19], Kapitel 3, S. 141–221[19] RUSZCZY ŃSKI, A. (Hrsg.) ; S HAPIRO, A. (Hrsg.): Handbooks in Operations Research and Management Science. Bd. 10: Stochastic Programming. Amsterdam :Elsevier, 2003[20] S EN, S. ; H IGLE, J.L.: The C 3 theorem and a D2 algorithm for large scale stochastic
mixed-integer programming: Set convexification. In: Mathematical Programming
104 (2005), Nr. 1, S. 1–20
[21] S WIDER, D. ; W EBER, C.: The costs of wind’s intermittency in Germany: application of a stochastic electricity market model. In: European Transactions on ElectricalPower 17 (2007), S. 151–172[22] VAN S LYKE, R.M. ; W ETS, R.: L-Shaped Linear Programs with Applications to
Optimal Control and Stochastic Programming. In: SIAM Journal of Applied Mathematics 17 (1969), Nr. 4, S. 638–663