1
Optimization of Dispersed Energy Supply Stochastic Programming with Recombining
Scenario Trees
Alexa Epe1 , Christian Küchler2 , Werner Römisch2 , Stefan Vigerske2 ,Hermann-Josef Wagner1 , Christoph Weber3 , and Oliver Woll3
1

Ruhr-Universität Bochum, Universitätsstraße 150, 44801 Bochum, Germany,{epe,lee}@lee.rub.de2
Humboldt–Universität zu Berlin, Unter den Linden 6, 10099 Berlin, Germany,{ckuechler,romisch,stefan}@math.hu-berlin.de
3
Universität Duisburg-Essen, Universitätsstraße 2, 45141 Essen, Germany,{Christoph_Weber,Oliver.Woll}@uni-duisburg-essen.de
Summary. The steadily increasing share of wind energy within many power generating systems leads to strong and unpredictable fluctuations of the electricity supply
and is thus a challenge with regard to power generation and transmission. We investigate the potential of energy storages to contribute to a cost optimal electricity
supply by decoupling the supply and the demand. For this purpose we study astochastic programming model of a regional power generating system consisting ofthermal power units, wind energy, different energy storage systems, and the possibility for energy import. The identification of a cost optimal operation plan allows
to evaluate the economical possibilities of the considered storage technologies.
On the one hand the optimization of energy storages requires the consideration of long-term planning horizons. On the other hand the highly fluctuating windenergy input requires a detailed temporal resolution. Consequently, the resulting
optimization problem can, due to its dimension, not be tackled by standard solution approaches. We thus reduce the complexity by employing recombining scenario
trees and apply a decomposition technique that exploits the special structure ofthose trees.

1.1 Introduction
Electric power, one of the most important fields within energy supply, has
two main characteristics: on the one hand supply and demand have to bebalanced at every time, on the other hand it is storable at only small rates. For
these reasons, power plants have to regulate any imbalances between supply
and demand, and, in particular, need to cope with unpredictable changes in
the customer load. For that purpose regulating power plants are used, whichmostly run in part load and with reduced efficiency. Alternatively fast power

2

Alexa Epe et. al.

plants such as open cycle gas turbines may be used, which can start up withinshort time. Beyond the cover of the fluctuating load of the customer side, thesepower plants must also adjust to the increasing share of time-varying power
production on the supply side, mostly caused from fluctuating renewables,
notably wind.
Germany is the country with the highest installed wind power capacities
worldwide. In the year 2006, there was approximately 20 GW installed (about
16.6% of the total installed power in Germany) and with the planned offshoredevelopment it could be up to 50 GW in 2030. Thereby, the sometimes strong
and rapid fluctuations of the wind energy fed into the electrical network as
well as the regional concentration in the north of the country increasingly pose
problems to the network operators and power suppliers [7, 17]. Conventional
fuel consumption may be saved by down-regulating conventional (back-up)
power plants, but investments in the power plant park can hardly be saved.
In this context, electrical energy storages offer a possibility to decouple
supply and demand and to achieve a better capacity utilization as well as ahigher efficiency of existing power plants. The changing context has led to anincreased interest in such possibilities over the last few years. Yet with theliberalization of the electricity markets, the economics of storages have to bevalued against market prices as established at the energy exchanges. Also theoperation of storages will mostly not follow local imbalances of demand andsupply, but rather try to benefit from market price variations. In particular,
the (partial) unpredictability of market prices as well as of wind energy supply
have to be taken into account. Things are complicated further through daily,
weekly, seasonal, and other cyclic patterns in demand, supply, and prices. Thisrequires a valuation of storages (and other options) over periods as long as
one year.
Cost optimal operation planning under uncertainty for such long time periods poses a huge challenge to conventional stochastic programming methods.
In this paper we investigate a novel approach, reducing complexity by applying recombining scenario trees. The latter are used to analyze a regionalenergy system model that is described in Section 1.2. Section 1.3 presents thedecomposition approach based on recombining trees, whereas Sections 1.4 and1.5 are devoted to the results obtained so far.

1.2 Model Description
To study the economics of storages, a fundamental model is used. Combining
technical and economical aspects, the model describes the energy supply of alarge city, the available technologies for electricity generation, and the demand.
An optimal load dispatch has to consider the marginal generation costs as wellas the impact of other system restrictions such as start up costs, etc. Most
important restriction of the model is the covering of the demand according
to a given profile. For this purpose, energy can be produced by conventional

1 Optimization of Dispersed Energy Supply

3

power plants, procured as wind energy, and purchased on the spot market, cf.
Figure 1.1.
External power supply

Wind power

Demandtime

time

Network

Customer
Offshorewind park

Conventional power
plant

StorageFig. 1.1. Scheme of the fundamental model

Uncertainty in the amount of available wind energy and electricity prices
is modelled by a multivariate stochastic process that can be represented bya recombining scenario tree. Thus, the proposed model combines many features of generation scheduling models (unit commitment and load dispatch) as
found typically in energy system models [14, 15]. In the following, the model
is discussed in detail. Table 1.1 gives an overview of the notation used.
Under the assumption of power markets with efficient information treatment and without market power, the market results correspond to the outcomes of an optimization carried out by a fully informed central planner. If
electricity demand is assumed to be price inelastic, welfare maximization isequivalent to cost minimization within the considered power network. Thereby,
the total costs T C are given as the sum of import costs ICt , operating costs
OCt,i , and startup costs SCt,i over all time steps t and unit types i:4Alexa Epe et. al.
Variables
Q
H
L

Production
Storage level
Capacity

IndicestTime step
T
Final timei
Unit typestu
Start-up
Parameters
D
DemandW
Wind powerℓLoad factor
η0 , ηm
Efficiency

IC
SC
OC
TC

Import costs
Start-up costs
Operating costs
Total costs

com
pum
imp
wind

Compressing power
Pumping power
Import power
Wind power

cstu
cimp
coth
cfuel

specific start-up costs
specific import costs
other variable costs
fuel price

Table 1.1. Notation used by the model.

TC =

T
X
t=1

ICt +

X

OCt,i + SCt,i

i

!

.

(1.1)

The costs for power import at time t are given by
ICt = cimp
Qimp
.
t
t

(1.2)

For the operating costs OCt,i , an affine function of the plant output Qt,i isassumed. An exact description of the plant operation costs requires a mixedbinary nonlinear formulation due to the dependency of the plant efficiency on
the power output and the startup behaviour. This is hardly feasible due tothe high level of time detail. An appropiate linearization can be done by defining an additional decision variable for each plant type, the capacity currently
online Lonlt,i [18]. The capacity online forms an upper bound on the actual
output. Multiplied with the minimum load factor, it is also a lower bound on
the output for each power plant. Hence, operating costs can be decomposedin fuel costs for operation at minimum load, fuel costs for incremental output,
and other variable costs:
OCt,i =

cfuel
cfuel
i,t
i,t
oth
ℓi Lonl(Qt,i − ℓi Lonlt,i +t,i ) + ci Qt,i .
0
ηi
ηim

(1.3)

Here, ηim denotes the marginal efficiency for an operating plant and ηi0 the
efficiency at the minimum load factor ℓi . With ηim > ηi0 , the operators havean incentive to reduce the capacity online (for details see [18]).
Besides operating costs, start-up costs may influence the power schedulingdecisions considerably. The start-up costs of unit i at time t are given bystu
SCt,i = cstu
i Lt,i ,(1.4)

1 Optimization of Dispersed Energy Supply

5

where Lstu
t,i is the start-up capacity given by
onl
onl
Lstu
t,i = max(0, Lt,i − Lt−1,i ).

(1.5)

Covering the demand at time step t is ensured by
X pum X
X
Qcom
Qt,i +
Qt,i + Qwind+ Qimp≥ Dt +
t,i ,
t
t
i

i

(1.6)

i

i.e., the supply at time t is given by the sum of the power production Qt,i , theimported energy Qimp
, and the wind energy supply Qwind. The total demandtt
equals the sum of the exogenously given domestic demand Dt and the pumping and compressing energies Qpum
and Qcomt,i used to fill the pumped hydrot,istorage and compressed-air storage, respectively.
The operation levels of the units, pumps, and air compressors are constrained by the available capacity,
Qt,i ≤ Lt,i ,

pum
Qpum
t,i ≤ Lt,i ,

com
Qcom
t,i ≤ Lt,i ,(1.7)

whereas the wind energy supply is bounded by the available wind energy attime t,Qwind≤ Wt .
(1.8)
t
For the storage plants, storage constraints need to be considered and thefilling and discharging has to be described. This leads to the following storagelevel equation, linking the storage level Ht,i at time t with the level Ht−1,i attime t − 1, both expressed in energy units. For the pumped hydro units, thisreads as11
1
Ht,i = Ht−1,i − m Qt,i −− m ℓi Lonl
(1.9)
t,i
ηi
ηi0
ηi
0,pum
onl,pum+ηim,pumQpum− ηim,pum )ℓi Lt,it,i + (ηifor t = 1, . . . , T , where H0,i denotes the initial fill level. Additionally, as anadequate terminal condition we require the initial and terminal fill levels of
the reservoirs to be fixed at the minimum fill level Himin . Further, the storagelevel at time step t is also limited by the minimum and maximum storagelevels,Himin ≤ Ht,i ≤ Himax .

(1.10)

Similar capacity constraints are formulated for the compressed-air units.
Additionally, all variables have to fulfill non-negativity conditions.
The objective of the optimization is to find a decision process satisfying theconstraints (1.5)–(1.10), being nonanticipative with respect to the stochastic
process (Wt , cimp
)t , and minimizing the expected total costs E[T C].
t

6Alexa Epe et. al.

1.3 Decomposition using Recombining Scenario Trees
In this section, we present the solution method based on recombining scenario
trees that has been developed in [11], and sketch a method for generating
recombining scenario trees.
1.3.1 Problem Formulation
The optimization problem presented in Section 1.2 can be written as a linear
multistage stochastic program:
" T
#
X
min E
hbt (ξ t ), xt i
(1.11)
t=1

s.t. xt ∈ Xt ,xt ∈ σ(ξ t ),At,0 xt + At,1 xt−1 = ht (ξ t ),t = 1, . . . , T,t = 2, . . . , T.

(1.12)

Thereby, the vector xt contains all decision variables at time stage t. The sets
Xt are closed and polyhedral and model deterministic, static linear constraints
at time t, i.e., the conditions (1.6), (1.7), and (1.10). The identities (1.12)
describe the random and time-coupling constraints (1.5), (1.8), and (1.9).
The uncertainty concerning the future wind energy input and spot prices ismodeled by the bivariate discrete time stochastic process ξ = (ξ t )t=1,...,T , thatenters into the optimization model through the costs bt (·) and the right-handsides ht (·), which are assumed to depend affinely linear on ξ t for t = 1, . . . , T .
Furthermore, ξ defines the nonanticipativity constraints, i.e., a decision xt attime t must depend exclusively on observations made until t. This is formalized
by the condition xt ∈ σ(ξ t ), where ξt denotes the vector (ξ 1 , . . . , ξt ).
To render possible a numerical solution of (1.11), every ξt is assumedtt}. Consequently, the pro, . . . , ξ(nto take values in a finite set Ξ t = {ξ(1)t)
cess ξ can be represented by a scenario tree, cf., e.g., [2]. Then, (1.11)
can be formulated as a (large scale) deterministic linear optimization problem that can be solved, in principle, by means of available solvers. However, with growing time horizon T , problem (1.11) becomes too large tobe solved as a whole and one has to resort to decomposition techniques,e.g., temporal decomposition. To this end, one considers certain time stages
0 = R0 < R1 < . . . < Rn < Rn+1 = T, and defines the cost-to-go function atRtime Rj and state (xRj , ξ(i)j ) ∈ XRj × Ξ Rj recursively by QRn+1 (·, ·) := 0 and
the Bellman Equation

1 Optimization of Dispersed Energy SupplyR

QRj (xRj , ξ(i)j ) :=

min E (QRj )


Rj+1

X

7

R

hbt (ξ t ), xt i + QRj+1 (xRj+1 , ξ Rj+1 ) ξRj = ξ(i)j 

t=Rj +1

s.t. xt ∈ Xt , xt ∈ σ(ξ t ),At,0 xt + At,1 xt−1 = ht (ξ t ),t = Rj + 1, . . . , Rj+1 ,t = Rj + 1, . . . , Rj+1 ,for j = 1, . . . , n. Using this notation, problem (1.11) can be reformulated interms of Dynamic Programming:
"R
#
1
X
R1
min E
(Q0 )
hbt (ξ t ), xt i + QR1 (xR1 , ξ )t=1

s.t. xt ∈ Xt ,xt ∈ σ(ξ t ),At,0 xt + At,1 xt−1 = ht (ξ t ),t = 1, . . . , R1 ,t = 2, . . . , R1 ,and solved by, e.g., the Nested Benders Decomposition method [1, 12, 16].
Furthermore, a modification of this algorithm as proposed in [11] allows to
exploit the structure of recombining scenario trees for simultaneous cutting
plane approximations. This approach indeed enables to solve problem (Q0 )for longer time horizons T and large number of scenarios.
1.3.2 Recombining Scenario Trees
At time t, the scenario tree representing ξ has nt = |Ξ t | nodes, that aretdenoted by u = 1, . . . , nt . The node u corresponds to the event {ξt = ξ(u)}.
A special situation is given whenever the subtrees associated at some nodes u
and k at time Rj coincide, i.e., the corresponding conditional distributions of(ξ t )t=Rj +1,...,T are equal:
h
i
h
i
R
R
P (ξ t )t=Rj +1,...,T ∈ · |ξ Rj = ξ(u)j = P (ξ t )t=Rj +1,...,T ∈ · |ξ Rj = ξ(k)j .

(1.13)

As far as it concerns the tree representation of the process ξ, property (1.13)
would allow to recombine the nodes u and k, and recombining at several
time stages Rj may prevent the node number to grow exponentially with thenumber of time stages. Unfortunately, recombining is not allowed under time
coupling constraints (1.12) since the scenario-dependent control xRj (ξ Rj ) willRR
not be equal on {ξRj = ξ(u)j } and {ξRj = ξ(k)j }, in general. However, (1.13)
R

can be useful since it entails equality of the cost-to-go functions QRj (·, ξ(u)j )R

and QRj (·, ξ(k)j ). This is exploited by the solution algorithm presented in Section 1.3.3.
In the remaining part of this section we sketch a method for generating
scenario trees with property (1.13) for some nodes and several time stages

8Alexa Epe et. al.

t=1

2

3

4

5

Fig. 1.2. Scenario tree with property (1.13), R1 = 3, and mR1 = 2, i.e., two different
subtrees are associated at time stage 3. (The black and the gray subtrees coincide,respectively.)

Rj , j = 1, . . . , n. It is a modification of the forward tree construction [9],also based on successive stagewise clustering of a set of sampled trajectories
ζ i = (ζ1i , . . . , ζTi ), i = 1, . . . , N , that coincide in t = 1. Basically, it consistsof constructing non-recombining subtrees for every time period [Rj + 1, Rj+1 ](Step 2), and assigning to several nodes at time Rj+1 the same subtree forthe subsequent time period (Step 1). Thereby, two nodes at time Rj+1 obtainthe same subtree h, whenever the values of ξ in these nodes are close forsome time t before Rj+1 . Table 1.2 explains the notation used. The valuesof mRj , nRj (h), and st+1 (u), determining the structure of the scenario tree,
may be predefined or, as proposed in [9], determined within the algorithm tonot exceed certain local error levels. Whenever the sampled trajectories comefrom a time series model, the parameter t may be chosen according to thelatter.
Algorithm 1 (Generation of a recombining scenario tree).
(1,1)
(1)
(1)
(h)
Initialization: Set C1:= C R0 := {1, . . . , N } , ξ1 := ζ11 and C Rj := ∅ forall h and j ≥ 1.
For j = 0, . . . , n (j-th recombination time stage):
1. If j > 0: Short-term history clustering for subtree assignment.
Find an index set A = {a1 , . . . , amRj } ⊂ {1, . . . , N } with minimalXX

h=1,...,mRj−1 u=1,...,nRj (h)

min

al ∈A

X

al
al
i
i
)−(ζR
, . . . , ζR
, . . . , ζR
)k.
k(ζR
j
j −t
j −t
j

(h,u)

i∈CR

j

Pass through all nodes (h, u) at time Rj :a) Consider an ah′ ∈ A that is close to node (h, u) in the sense that

1 Optimization of Dispersed Energy Supply

9

ζtivalue of trajectory i at time tξt,(u) value of the random variable ξ t in node unumber of subtrees with root node at time Rj (mR0 := 1)m Rj(h, u) node u of some subtree hnRj (h) number of nodes at time Rj of some subtree hst+1 (u) number of nodes at time t + 1 descending from some node u at time tttime parameter for short-term history clustering(h,u)
subset of {1, . . . , N }, indicating trajectories ζ i
Ctgoing at time t through node (h, u)(h)subset of {1, . . . , N }, indicating trajectories ζ i
C Rj
lying in subtree h with root node at time Rj
Table 1.2. Notation used by Algorithm 1.

ah′ ∈ arg min
al ∈A

X

al
al
i
i
) − (ζR, . . . , ζR
, . . . , ζR
)k.
k(ζR
j
j −t
j −t
j

(h,u)

i∈CR

j

b) Node (h, u) obtains subtree h′ , i.e., one has to update
(h′ )

(h′ )

(h,u)

C Rj := C Rj ∪ CRj.

In particular, all trajectories ζ i belonging to node (h, u) will be used
for the construction of subtree h′ of the subsequent timeperiod.
2. Subtree generation.
For every subtree h = 1, . . . , mRj of period [Rj + 1, Rj+1 ]:(h)a) Find an index set A = {a1 , . . . , anRj (h) } ⊂ C Rj with minimalXali
min kζR
− ζR
k.
j +1
j +1

(h)

i∈C R

al ∈A

j

(h,u)

(h)

Find a partition CRj +1 , u = 1, . . . , nRj +1 (h), of C Rj with
(h,u)

(h)

al
i
− ζR
k}.
CRj +1 ⊂ {i ∈ C Rj : au ∈ arg min kζR
j +1
j +1
al ∈A

au
Define the value of ξRj +1 on node (h, u) by ξRj +1,(h,u) := ζt+1.
For every subtree h̃ of the preceding period [Rj−1 + 1, Rj ] and every
node (h̃, ũ) at time Rj : Define the transition probability from node(h̃, ũ) to node (h, u) by (h,u) |CRj +1 | if C (h̃,u) ⊂ C (h) ,
(h)
Rj
Rj
|CR |
PRj +1|Rj [(h, u)|(h̃, ũ)] :=
j

0
else.

10

Alexa Epe et. al.

b) For t = Rj + 1, . . . , Rj − 1:
For every node ũ of subtree h at time t:(h,ũ)
Find an index set A = {a1 , . . . , ast+1 (ũ) } ⊂ Ct
with minimalX(h,ũ)ali
min kζt+1− ζt+1
k.

al ∈A

i∈Ct
(h,u)

(h,ũ)

Find a partition Ct+1 , u = 1, . . . , st+1 (ũ), of Ct(h,u)(h,ũ)

Ct+1 ⊂ {i ∈ Ct

with

al
i
− ζt+1: au ∈ arg min kζt+1
k}.
al ∈A

au
Set ξt+1,(h,u) := ζt+1
and Pt+1|t [(h, u)|(h, ũ)] :=(h,u)|Ct+1 |(h,ũ)|Ct|.

The determination of the index sets A is a k−mean problem and, thus, an
NP-hard combinatorial optimization problem [5]. While it is possible for smallvalues mRj , nRj (h), and st+1 (u) to find optimal sets A by enumeration, larger
values demand for heuristics, e.g., the forward selection proposed in [3, 8].
In Step 1 of Algorithm 1, several nodes (h̃, ũ) at time Rj obtain the same
subtree h. For notational convenience, we pick out a representative amongst
Rthem and denote the associated value of ξ Rj by λh j . The following function
will be used in Section 1.3.3 and maps a node (h̃, ũ) with subtree h to thecorresponding representative node:RRjj
λRj : Ξ Rj → {λ1 j , . . . , λRmR } =: Λ ,
j

λ

Rj

R
Rj
(h)
(h̃,ũ)) := λh j whenever CRj ⊂ C Rj .
(ξ(h̃,ũ)1.3.3 Solution Algorithm
In [11], it was shown how to modify a Nested Benders Decomposition [1, 12, 16]of problem (1.11) to exploit the recombining property (1.13) of the process ξ.
In the following, we sketch this modified algorithm.
Let us consider the formulation (Q0 ) of problem (1.11). A Nested Benders
Decomposition successively approximates the piecewise-linear convex funcRtions xRj 7→ QRj (xRj , ξ(u)j ) by a set of supporting hyperplanes and evaluatesthem in an adaptively chosen sequence of points xRj . Whenever two nodes
RRu and k at time Rj fulfill (1.13), the functions QRj (·, ξ(u)j ) and QRj (·, ξ(k)j )coincide, and, thus, they may be approximated simultaneously.
To this end, we define the following underestimating functions: We set
Rj
Rj
let
QLC
Rn+1 (·, ·) := 0 and for j = n, . . . , 0, x̄Rj ∈ XRj , and λi ∈ Λ1 Optimization of Dispersed Energy Supply

11

R

j
(QL
QL
Rj (x̄Rj , λi ) :=
Rj )


Rj+1
X
R
Rj+1 Rj+1
min E )) ξ Rj = λi j (ξ
hbt (ξt ), xt i + QLC
Rj+1 (xRj+1 , λt=Rj +1

s.t. xt ∈ Xt ,xt ∈ σ(ξ t ),t = Rj + 1, . . . , Rj+1 ,At,0 xt + At,1 xt−1 = ht (ξ t ),xRj = x̄Rj .

t = Rj + 1, . . . , Rj+1 ,(1.14)

R

R

j+1
j+1
Thereby, QLC) is an approximation of QL) by supporting
Rj+1 (·, λi
Rj+1 (·, λi
hyperplanes that is easy to evaluate and that will be properly defined in equation (1.15) below. Problem (QLRj ) is often referred to as the master problem.
Note, that in contrast to the classical Nested Benders Decomposition, the same
Rj+1
) can be used in the objective function of (QL
approximation QLC
Rj )Rj+1 (·, λiRfor all nodes with the same subtree i, i.e., whenever λRj+1 (ξ Rj+1 ) = λi j+1 .
Thus, this decomposition into subproblems for each timeperiod allows one to
exploit the recombining nature of the process ξ. This had not been possible
with other decomposition algorithms like scenario decomposition [12].
Rj
The function QLC
Rj (·, λi ) is used to induce a feasible solution at stage RjRj
and to approximate the value of QL
Rj (·, λi ) on its domain. For the latterRjpurpose, given a point x̄ ∈ XRj with QL
Rj (x̄, λi ) < ∞, an optimality cut

R

R

j
j
L
supporting QL
Rj (·, λi ) is given by QRj (x̄, λi ) + hπ, xRj − x̄i ≤ 0, where
π denotes the dual variables corresponding to the constraint (1.14) in an
optimal solution of problem (QLRj ). To induce feasibility at time stage Rj , apoint x̄ ∈ XRj that is infeasible for (QLRj ) is cut off using a feasibility cut
hd, xi + e ≤ 0. This cut is computed by solving an auxiliary problem, cf. [11],and has the property hd, x̄i + e > 0 and hd, xi + e ≤ 0 for all x ∈ XRj withRj
QL
Rj (x, λi ) < ∞.

RjHence, an approximation of QL
Rj (·, λi ) by means of optimality cuts

R

R

Copt (λi j ) and feasibility cuts Cfeas (λi j ) is given byRj
QLC
Rj (xRj , λi ) :=

R

max

R
(x̄,π̄)∈Copt (λi j )

s.t.

j
QL
Rj (x̄, λi ) + π̄, xRj − x̄d, xRj + e ≤ 0,

(1.15)R(d, e) ∈ Cfeas (λi j ).

The solution algorithm processes the master problems (QLRj ), j = 0, . . . , n,of the decomposed scenario tree in a forward or backward manner. At eachRj
Rj
∈ ΛRj , is evaluatedtime stage Rj , each master problem QL
Rj (·, λi ), λi
R

R

R

j
j
L
for a set Zj (λi j ) of controls xRj . If QLC
Rj (xRj , λi ) < QRj (xRj , λi ), the

R

R

j
j
LC
approximation QLC
Rj (·, λi ) (and all master problems that use QRj (·, λi ))is updated by generating new optimality or feasibility cuts. Further, in the12Alexa Epe et. al.

forward mode, new control points xRj+1 are generated from the solution ofR
Rj+1
the master problem (QL) for λi j+1 ∈ ΛRj+1 .
Rj ) to form the sets Zj+1 (λiSince each such evaluation contributes several new controls xRj+1 to the sets
R
Zj+1 (λi j+1 ), the latter can grow exponentially with increasing j. It was shownin [11] how the problem structure allows to deal with this difficulty.
The algorithm stops when either the first timeperiod master problem (QL0)
is infeasible, or all master problems could be solved to optimality and the generation of cuts has stopped. In the former case, also problem (Q0 ) is infeasible,
in the latter, the problem has been solved to optimality. Another stopping criteria which allows to stop the algorithm when the error falls below a given
tolerance, is also discussed in [11]. A more detailed description of the Nested
Benders Decomposition Algorithm can be found in [1, 6, 12].

1.4 Case study
We study a power generating system, consisting of a hard coal power plant
to cover the minimum and medium load, and two fast gas turbines on different power levels to cover the peaks. The operating parameters of these units
rely on real data. Furthermore, the model contains an offshore wind park, a
pump-storage power plant (PSW) with the basic data of the PSW Geesthacht,
Germany, and a compressed-air energy storage (CAES) with the operating parameters of the CAES Huntorf, Germany. Further source of power supply is
the EEX spot market. The time horizon considered for the optimization is
one year and a hourly discretization is used, i.e., the model contains T = 8760
time stages.
The stochastic wind power process is represented by a time series model
fitted to historical data and scaled to the size of the offshore wind park regarded. To take into account the interdependency between wind power andspot price behaviour, the expected spot market prices are calculated from afundamental model that is based on the existing power plants in Germany
and their reliability, prices for fuels and CO2 , the German load, and the wind
power process above. Fluctuation of the spot prices around their expectedvalue are modeled by a further time series model. This hybrid approach wasused to generate 1 000 trajectories, containing hourly values of wind power andspot prices in the course of one year. These trajectories were used to generatea recombining tree by Algorithm 1 of Section 1.3.2. The resulting scenario tree
branches three times per day in a binary way. Recombination into mRj = 3
different subtrees took place once a day, i.e., Rj = j · 24, j = 1, . . . , 364.

1.5 Numerical Results
The optimization problem was solved with varying model parameters. To thisend, a base setting was defined, with wind power of approximately 50% of the1 Optimization of Dispersed Energy Supply

13

100%
1879264

Load / Maximum load

75%1409448

Wind
Import
Gas
Coal
Load

50%
93963225%4698160
1

Monday25

Tuesday 49 Wednesday 73 Thursday

97

Friday121 Saturday

145 Sunday

169

Fig. 1.3. Optimal power scheduling in a winter week.

80

Power price
CAES discharge
CAES level70

100,00%60€ / MW50

75,00%

40
50,00%
30

20
25,00%
10

0
1 Monday

25 Tuesday

49Wednesday

73 Thursday

97

Friday121 Saturday

145 Sunday

0,00%
169

Fig. 1.4. Spot market price and CAES output in a winter week.

totally installed plant power and storage sizes corresponding to the aforementioned CAES and PSW units. Coming from this setting, variations with higher
and lower levels of installed wind power and different storage dimensions were
calculated. In the following some results are presented.
The optimal operation levels along a randomly chosen scenario from thebase setting during a winter week are depicted in Figure 1.3. Whenever the14Alexa Epe et. al.

power production exceeds the demand curve, energy is put into the storages,
whereas the white spaces under the demand curve represent the output of
the storage plants. The operation levels of the thermal units show the usual
characteristics and availability of wind power obviously reduces imports fromthe spot market. The storage units are mainly used to cover the peaks and are
only marginally used during the weekend. In this model, the contribution of
the operating costs to the power supply costs amount to 2.08 Eurocents/kWhwith using storage plants and 2.10 Eurocents/kWh without using storage
plants. Figure 1.4 shows the optimal output and fill level of the CAES (as afraction of maximum discharge power and maximum fill level, respectively)in comparison to the actual power price. The minimum fill level of the CAES
is 60%. Obviously, the storage plant discharges in times of high spot prices
on weekdays. The aforementioned marginal usage of storage plants during theweekend coincides with lower power prices over this period.

100%
1879264

Load / Maximum load

75%1409448

Wind
Import
Gas
Coal
Load

50%
93963225%4698160
1

Monday25

Tuesday 49 Wednesday73 Thursday 97

Friday121 Saturday

145 Sunday

169

Fig. 1.5. Optimal power scheduling in a winter week with doubled wind power
capacity.

To study the impact of the share of wind power on the system, the optimization problem was solved again with doubled wind power capacity. The
results along the same scenario and for the same winter week are depicted in
Figure 1.5. While this extension does not lead to significant changes of thethermal units, it enables to largely reduce the amount of energy bought at
the spot market. Figure 1.6 shows the operation of the CAES in the course of
the week. Again, the CAES is mainly used at peak times to avoid expensive

1 Optimization of Dispersed Energy Supply

1580

Power price
CAES discharge
CAES level70

100,00%60€ / MW50

75,00%

40
50,00%
30

20
25,00%
10

0
1

Monday25

Tuesday

49 Wednesday

73

Thursday97

Friday121

Saturday145

Sunday0,00%
169

Fig. 1.6. Spot market price and CAES output in a winter week with doubled wind
power capacity.

imports from the spot market. It can be seen that the availability of morewind power in the system can lead both to more and to less extraction of
stored energy. This is due to the fact that, on the one hand, with more wind
power more energy may be stored and therefore extracted (Sunday). On theother hand, less power has to be generated in times with high wind power(from Wednesday to Friday).
The optimization problem was solved further times with varying quantities
of installed wind power and storage capacities. Figure 1.7 shows the minimal
expected costs depending on the wind power capacity for different storage
capacities. Thereby, a wind factor of y stands for an amount of wind power
being y times the wind power of the base setting. In relation to the wind power
capacity, the impact of an extension of the storage system on the costs appears
to be rather marginal and the individual curves are almost superposed. Thus,
to analyze the latter, Figure 1.8 shows the relative reduction of costs that canbe achieved by the use of storage systems of different dimensions, where amodel without storages generates operating costs of 100%. Again, a storage
system dimension of y corresponds to y times the dimension of the base setting.
The results clearly show, that the relative cost reduction due to storage use is
the highest in the twice-wind-setting, and in all settings the most prevailing
gradient is between no use and the use of the half dimension of storage sizes.
Hence, it seems promising to study expansion models for cost-optimal storage
sizes, taking into account operational as well as investment costs.
The optimization algorithm was implemented in C++ and the masterproblems were solved with CPLEX 10.0 [10]. Running time on a PC with 2.4
GHz CPU and 2 GB RAM was 10 minutes, approximatively.

16

Alexa Epe et. al.
130,00%relative costs

120,00%110,00%100,00%90,00%
80,00%
70,00%
60,00%

Wind x 0.5

Wind x 1.0

Wind x 1.5

Wind x 2.0

Fig. 1.7. Minimal expected costs depending on the wind power capacity installed
for different storage capacities installed.

2,50%

relative cost reduction

wind x 2.0
wind x 1.5

2,00%

wind x 1.0
wind x 0.75

1,50%wind x 0.5

1,00%

0,50%

0,00%
0

0.5

1

1.5

2

Storage system dimension
Fig. 1.8. Reduction of minimal expected costs depending on the storage capacity
installed for different wind power capacities installed.

1.6 Conclusions and Outlook
We applied a decomposition method for linear multistage stochastic optimization problems proposed by [11] to optimal scheduling within a regional energy

1 Optimization of Dispersed Energy Supply

17

system including wind power and energy storages. It has been shown that
this approach relying on recombining scenario trees allows to handle multistage problems with large numbers of scenarios and including time-coupling
constraints, and, therefore, it is suitable for optimizing and analyzing energy
systems.
In principle, the recombining tree decomposition approach allows for discrete decision variables in the first time stage and, hence, this method seems
to be also appropriate to find optimal first-stage investment decisions within
expansion models. This could be one aspect of future studies. However, further research is needed to extend the decomposition approach to more general
optimization models, in particular those including discrete variables in latertime stages. The latter would allow to adapt numerous aspects of the energy system model to achieve a more detailed picture of the system and its
constraints.
Another aspect of future research could be the extension of the decomposition approach to optimization problems including multiperiod risk measures.
For this purpose, especially the class of polyhedral risk measures [4] seems tobe suitable, since the linear structure of the optimization problem is maintained.
Acknowledgement
This work was supported by the German Ministry of Education and Research (BMBF) within the topic “Dezentrale regenerative Energieversorgung:Innovative Modellierung und Optimierung” and the “Wiener Wissenschafts-,
Forschungs- und Technologiefonds” in Vienna (http://www.univie.ac.at/
crm/simopt).

References
1. J.R. Birge. Decomposition and partitioning methods for multistage stochastic
programming. Operations Research, 33(5):989–1007, 1985.
2. J. Dupačová, G. Consigli, and S.W. Wallace. Scenarios for multistage stochastic
programming. Annals of Operations Research, 100:25–53, 2000.
3. J. Dupačová, N. Gröwe-Kuska, and W. Römisch. Scenarios reduction in stochastic programming: An approach using probability metrics. Mathematical Programming, 95(A):493–511, 2003.
4. A. Eichhorn and W. Römisch. Polyhedral risk measures in stochastic programming. SIAM Journal on Optimization, 16:69–95, 2005.
5. M.R. Garey and D.S. Johnson. Computers and Intractability - A Guide to the
Theory of NP-Completeness. W.H. Freeman, 1979.
6. H.I. Gassmann. MSLiP: a computer code for the multistage stochastic linear
programming problem. Mathematical Programming, 47:407–423, 1990.
7. E. Handschin, F. Neise, H. Neumann, and R. Schultz. Optimal operation ofdispersed generation under uncertainty using mathematical programming. International Journal of Electrical Power & Energy Systems, 28:618–626, 2006.

18

Alexa Epe et. al.

8. H. Heitsch and W. Römisch. Scenario reduction algorithms in stochastic programming. Computational Optimization and Applications, 24:187–206, 2003.
9. H. Heitsch and W. Römisch. Scenario tree modeling for multistage stochastic
programs. Mathematical Programming, to appear, 2008.
10. ILOG, Inc. CPLEX 10.0. http://www.ilog.com/products/cplex.
11. C. Küchler and S. Vigerske. Decomposition of multistage stochastic programswith recombining scenario trees. Stochastic Programming E-Print Series, 9,2007. http://www.speps.org.
12. A. Ruszczyński. Decomposition Methods, chapter 3, pages 141–221. In
Ruszczyński and Shapiro [13], 2003.
13. A. Ruszczyński and A. Shapiro, editors. Stochastic Programming. Handbooksin Operations Research and Management Science. Elsevier, Amsterdam, 2003.
14. D. Swider, P. Vogel, and C. Weber. Stochastic model for the european electricity
market and the integration costs for wind power. Technical report, GreenNet
Report on WP 6, 2004.
15. D. Swider and C. Weber. The costs of wind’s intermittency in Germany: Application of a stochastic electricity market model. European Transactions on
Electrical Power, 17(2): 151–172.
16. R.M. Van Slyke and R. Wets. L-shaped linear programs with applications to
optimal control and stochastic programming. SIAM Journal of Applied Mathematics, 17(4):638–663, 1969.
17. H.-J. Wagner. Wind Energy and Present Status in Germany. In: N. Bansal and
J. Mathur (Editors), Wind Energy Utilization. Anamaya Publishers, New Delhi,
2002.
18. C. Weber. Uncertainty in the Electric Power Industry: Methods and Models for
Decision Support. New York et al., 2005.