Scenario Reduction Algorithms in Stochastic Programming.
We consider convex stochastic programs with an (approximate) initial probability distribution P having finite support supp P, i.e., finitely many scenarios. The behaviour of such stochastic programs is stable with respect to perturbations of P measured in terms of a Fortet-Mourier probability metric. The problem of optimal scenario reduction consists in determining a probability measure that is supported by a subset of supp P of prescribed cardinality and is closest to P in terms of such a probability metric. Two new versions of forward and backward type algorithms are presented for computing such optimally reduced probability measures approximately. Compared to earlier versions, the computational performance (accuracy, running time) of the new algorithms has been improved considerably. Numerical experience is reported for different instances of scenario trees with computable optimal lower bounds. The test examples also include a ternary scenario tree representing the weekly electrical load process in a power management model.
Introduction
Many stochastic decision problems may be formulated as convex stochastic
programs of the form
minf
Z
where X  R m is a given nonempty closed convex
set,
a closed subset
of R s , the function f 0
from
R m to R is continuous with respect to !
and convex with respect to x, and P is a xed Borel probability measure on
P(
). For instance, this formulation covers (convex) two- and
multi-stage stochastic programs with recourse.
Typical integrands f 0 (; x), x 2 X , in convex stochastic programming problems
are nondierentiable, but locally Lipschitz continuous
on
. In the fol-
lowing, we assume that there exist a continuous and nondecreasing function
nondecreasing function
and some xed element s such that
Heitsch, Romisch
for each x 2 X , where the function c
R is given by
c(!; ~
This means that the function h(k  ! 0 describes the growth of the local
Lipschitz constants of f 0 (; x) in balls around ! 0 with respect to some norm
k on R s . An important particular case is that the function h grows poly-
nomially, i.e.,  1. For instance, it is
shown in [11] that the choice is appropriate for two-stage models with
stochasticity entering prices and right-hand sides.
It is shown in [4,11] that the model (1) behaves stable with respect to small
perturbations in terms of the probability metric
c (P; Q) := sup
Z
Z
where F c is the class of continuous functions dened by
!) for all !; ~
and probability measures P and Q in the set
P(
Z
(see also the earlier work in [13]). The distance  c is a probability metric on
with -structure and also called a Fortet-Mourier (type) metric. In this
generality, it is introduced in [14] and further studied in [10,12]. In particular,
the metric  c has dual representations in terms of the Kantorovich-Rubinstein
functional (cf. Section 5.3 in [10] and [12]).
An important instance is that the initial probability measure P is itself discrete
with nitely many atoms (or scenarios) or that a good discrete approximation
of P is available. Its support may be very large so that, for reasons
of computational complexity and time limitation, this probability measure
is further approximated by a probability measure Q carried by a (much)
smaller subset of scenarios. In this case the distance  c (P; Q) represents the
optimal value of a nite-dimensional linear program. More precisely, from (4)
we obtain for
c
where J
P(
denotes the Dirac measure placing unit
mass at !. In particular, the metric  c can be used to evaluate distances of
specic probability measures obtained during a scenario-reduction process.
Various reduction rules appear in the literature in the context of recent large-scale
real-life applications. We refer to the corresponding discussion in [4], to
Scenario Reduction Algorithms in Stochastic Programming 3
the recent work [3] on scenario generation and reduction, and to the paper
[9], in which an approach to scenario generation based on Fortet-Mourier
distances is given.
In the present paper, we follow the approach for reducing scenarios of a given
discrete probability measure
developed in [4]. It consists in
determining an index set J   of given cardinality #J
and a probability measure Q
such that
Problem (7) may be reformulated as a mixed-integer program. In Section 2
we derive bounds for (7) and develop two new algorithms (fast forward selection
and simultaneous backward reduction), which constitute heuristics for
solving (7). We study their complexity and their relations to the algorithms
in [4]. Indeed, the fast forward selection algorithm turns out to be an ecient
implementation of the forward selection procedure of [4], producing the
same reduced probability measures. In order to compare the performance of
the algorithms we provide, in Section 3, explicit formulas for the minimal
distances (7) in case that P is a regular (binary or ternary) scenario tree
(i.e., a tree having a specic structure) and Q  is a reduced tree with xed
cardinality n. In Section 4 we report on numerical experience for the reduction
of regular binary and ternary scenario trees. The test trees also include
the ternary scenario tree representing the weekly electrical load process in a
power management model which was considered in [4]. It turns out that the
new implementation of the fast forward selection algorithm is about 10-100
times faster than the earlier version. Furthermore, fast forward selection is
the best algorithm when comparing accuracy. The results of the simultaneous
backward reduction algorithm are more accurate than the backward reduction
variant of [4] in most cases, but at the expense of higher running times.
When comparing running times, fast forward selection (simultaneous backward
reduction) is preferable in case that n < N
Reduction
We consider the stochastic program (1) and select the function c of form (3)
such that the Lipschitz condition (2) is satised. Let the initial probability
distribution P be discrete and carried by nitely many scenarios ! iwith
weights
and consider the probability
measure
4 Heitsch, Romisch
i.e., compared to P , the measure
reduced by deleting all
by assigning new probabilistic weights q j to each scenario
. The optimal reduction concept described above recommends
to consider the probability distance
depending on the index set J and q. The optimal reduction concept (7)
says that the index set J  and the optimal weight q  are selected such
that D(J  ; q
Jg. First we recall the following bounds for min q D(J ;
when the index set J
is xed ([4], Theorem 3.1).
Theorem 2.1 (redistribution)
For any index set J
where C := max
Furthermore, we have equality in (9) and, hence, optimality of  q if h  1.
For convenience of the reader the proof of Theorem 2.1 is displayed in the


Appendix

. The interpretation of the formula (10) is that the new probability
of a preserved scenario is equal to the sum of its former probability and of
all probabilities of deleted scenarios that are closest to it with respect to c.
If h  1, we call (10) the optimal redistribution rule.
Next we discuss the optimal choice of an index set J for scenario reduction
with xed cardinality #J . Theorem 2.1 motivates to consider the following
formulation of the optimal reduction problem for given n 2 N, n <
minfD J :=
ng (11)
Problem (11) means that the set has to be covered by two sets
such that J has xed cardinality N n and
the cover has minimal cost D J . Thus, (11) represents a set-covering problem.
It may be formulated as a 0-1 integer program (cf. [7]) and is NP-hard. Since
e-cient solution algorithms are hardly available in general, we are looking
Scenario Reduction Algorithms in Stochastic Programming 5
for (fast) heuristic algorithms exploiting the structure of the costs D J . In the
specic cases of solving (11) becomes quite easy.
In case that #J = 1, the problem (11) takes the form
min
l2f1;:::;Ng
l min
If the minimum is attained at l  i.e., the scenario ! l  is deleted,
the redistribution rule (10) yields the probability distribution of the reduced
measure
Q. If j  2 arg min j 6=l  c(! l  holds that  q
and  q g. Of course, the optimal deletion of a single
scenario may be repeated recursively until a prescribed number N n of
scenarios is deleted. This strategy recommends a conceptual algorithm called
backward reduction.
In case that #J = N 1, the problem (11) is of the form
min
u2f1;:::;Ng
If the minimum is attained at u  only the scenario ! u is kept
and the redistribution rule (10) provides
strategy provides the basic concept of a second conceptual algorithm called
forward selection.
First, we take a closer look at the backward reduction strategy. A backward
type algorithm was already suggested in [4,6]. It determines a reduced scenario
set by reducing N n scenarios from the original set of scenarios as
follows. Let the indices l i be selected such that
l min
It can be shown that
is a lower bound of the optimal value of (11). Furthermore, it holds that the
index set is a solution of (11) if for all n, the set
is nonempty ([4,6]). This
property motivates the following algorithm. In the rst step, an index n 1 with
determined using formula (14) such that J
is a solution of (11) for Next, the redistribution rule of Theorem 2.1 is
used. This leads to the reduced probability measure P 1 containing all scenarios
indexed by reduced
by deleting all scenarios belonging to some index set J 2 with #J
and n  which is obtained in the same way using formula (14).
This procedure is continued until, in step r, we have n
6 Heitsch, Romisch
Finally, the redistribution rule (10) is used again for the index set J . This
algorithm is called backward reduction of scenario sets. There are still many
degrees of freedom to choose the next scenario in each step. Often there exist
several candidates for deletion. In Section 4 we use one particular implementation
of backward reduction of scenario sets.
Another particular variant consists in the case that #J
This variant (without the nal redistribution) was already announced
in [2,5]. However, numerical tests have shown that the backward
reduction of scenario sets provides slightly more accurate results compared
to backward reduction of single scenarios.
Next we are going to present a new modication of the backward reduction
principle. The major dierence is to include all deleted scenarios into each
backward step simultaneously. Namely, the next index l i is determined as a
solution of the optimization problem
A more detailed description of the whole algorithm, which will be called
simultaneous backward reduction, is given in
Algorithm 2.2 (simultaneous backward reduction)
Sorting of fc
c [1]
ll := min
z [1]
l := p l c [1]
l2f1;:::;Ng
z [1]
Step i: c [i]
kl := min
z [i]
l :=
z [i]
Step N-n+1: Redistribution by (10):
Algorithm 2.2 allows the following interpretation. Its rst step corresponds
to the optimal deletion of only one scenario. For i > 1, l i is chosen such that
where D J [i 1] [flg is dened in (11). Hence, the index l i is dened recursively
such that the index set is optimal subject to the constraint
Scenario Reduction Algorithms in Stochastic Programming 7
that the previous indices are xed.
Since running times are important characteristics of scenario reduction al-
gorithms, we study the computational complexity, i.e., the number of elementary
arithmetic operations, of Algorithm 2.2. It is shown in [6] that a
proper implementation (without sorting) of backward reduction of scenario
sets requires a complexity of O(N 2 ) operations which holds uniformly with
respect to n. When comparing formulas (14) and (16), one notices an increase
of complexity in the cost structure of (16) for determining l i . More precisely,
step i requires the computation of N sums each consisting of i summands
and N i+1 comparisons. Each summand represents a product of two
numbers. One of these factors requires about 2 operations for determining
the minimum. The sorting process in step 1 requires O(N 2 log N) operations
([1], Chapter 1). When excluding the complexity of evaluating the function
c and of the redistribution rule, altogether we obtain
where a(N) := N 3+O(N 2 log N)
operations for selecting a subset of n scenarios. Hence, we have
Proposition 2.3 The computational complexity for reducing a set of N 2 N
scenarios to a subset containing n 2 consists of b N (n)
(see (18)) operations when using simultaneous backward reduction.
Hence, the complexity of simultaneous backward reduction is increasing with
decreasing n and is, of course, minimal at This result corresponds to
the running time observations of our numerical tests reported in Section 4.
Next, we describe a strategy that is just the opposite of backward reduction.
Its conceptual idea is based on formula (13) and consists in the recursive
selection of scenarios that will not be deleted. The basic concept of such
an algorithm is given in [4] and called forward selection. Forward selection
determines an index set fu un g such that
g. The rst step of
this procedure coincides with solving problem (13). After the last step, the
optimal redistribution rule has to be used to determine the reduced probability
measure. Formula (19) allows the same interpretation as in case of
simultaneous backward reduction. It is again closely related to the structure
of D J in (11). Now, let us consider the following algorithm, which is easy to
implement and is called fast forward selection.
8 Heitsch, Romisch
Algorithm 2.4 (fast forward selection)
z [1]
u2f1;:::;Ng
z [1]
Step i: c [i]
z [i]
z [i]
Step n+1: Redistribution by (10):
Theorem 2.5 The index set fu un g determined by Algorithm 2.4 is a
solution of the forward selection principle, i.e., u i satises condition (19) for
each
holds for each
D J [i] is dened in (11).
Proof: For the result is immediate. For holds that
z [i]
Hence, the index u i satises condition (19) and it holds that
z [i]
The conditions (17) and (20) show that both algorithms are based on the
same basic idea for selecting the next (scenario) index. The only dierence
Scenario Reduction Algorithms in Stochastic Programming 9
consists in the use of backward and forward strategies, respectively, i.e., in
determining the sets of deleted and remaining scenarios, respectively.
As in the case of backward reduction, the computational complexity of Algorithm
2.4 is of interest. Step i requires operations for computing
c [i]
operations for z [i]
operations for determining u i . Altogether, we obtain
operations for selecting a subset of n scenarios. Hence, we have
Proposition 2.6 The computational complexity of fast forward selection for
reducing a set of N 2 N scenarios to a subset containing n 2
scenarios consists of fN (n) (see (21)) operations.
Hence, the complexity of fast forward selection is increasing with increasing
n and is maximal if . Thus, the use of fast forward selection is recommendable
if the number n of remaining scenarios satises the condition
(n). The number n such that fN (n
is a zero of a polynomial of degree 3 which depends nonlinearly on N . It turns
out that n   Nfor large N .
3 Minimal distances of scenario trees
All algorithms discussed in the previous section provide only approximate
solutions of (11) in general. Since error estimates for these algorithms are not
available, we need test examples of discrete original and reduced measures
of dierent scale with known (optimal)  c -distances. Because of their practical
importance, we consider probability measures with scenarios exhibiting a
tree structure. In particular, we derive optimal distances of certain regularly
structured original scenario trees and of their reduced trees containing different
numbers of scenarios.
We consider a scenario tree that represents a discrete parameter stochastic
process with a parameter set f0; and with scenarios
(or paths) branching at each parameter k 2 f0; branching
degree d (i.e., each node of the tree has d successors). In case of
the tree is called binary and ternary, respectively. Hence, the tree
consists of N := d K scenarios
, and has
d K as its root node. Furthermore, we let all scenarios have equal
probabilities
. Such a scenario tree is called regular if,
for each k 2 there exist symmetric sets V k := f- k
d g  R
such that
Heitsch, Romisch
where a (K + 1)-tuple of indices (i corresponds
to each index We say that
In case of means that the sets V k are
of the form V respectively, for some
Clearly, it holds that - 0
trees. Figure 1
c
c
c c
c
c c
c c
c
c
@
@ @
@
@ @
@
@ @
@
@ @

Figure

1: Binary scenario tree
shows an example of a regular binary scenario tree with
scenarios. We specify the function c in (3) by setting h  1 and by choosing
the maximum norm k  k1 on R K+1 , i.e.,
c(!; ~
Our rst result provides an explicit formula for the minimal distance between
a regular binary tree and reduced subtrees with at least
Proposition 3.1 (3/4-solution)
Let a regular binary scenario tree with
it holds for each n 2 N with N
Proof: We use the representation (22) of each scenario ! i for
Let the
corresponding (K + 1)-tuples of indices. Let l 2 Kg be such that
l 1. Then we obtain
(- r
jr )j
l
(- r
Scenario Reduction Algorithms in Stochastic Programming 11
Hence, it holds for each J  that
It remains to show that there exists an index set J  such that #J
and that the lower bound is attained, i.e., D
. To this end, we
consider the index set
I  := fi
I  . Let T r  denote the tree consisting of all
I  . Figure 2 illustrates a detail of T r  starting at a node
r
r
r
r
r r
@ @
@ @
@ @

Figure

2: Detail of the subtree T r
at level k 0 1 and ending at level k 2. Hence, for the cardinality of I  and
J  we obtain that
4 and #J
Now we want to show that for each j 2 J  , there exists an index i 2 I  such
that be the related scenario. Let
us consider the behaviour of ! j on the branching levels k 2.
we have to distinguish three cases each for - k0
Case
Case
Case (3): - k0+1
Now, we consider the following (K
all k 62 fk
12 Heitsch, Romisch
denote the corresponding index. Clearly, i 2 I  and,
consequently, it holds for the distance between ! i and ! j that
(- r
jr )j
(- r
jr )j
case (1)
case (2)
maxfj2- k0 j; j2- k0 2- k0+1 jg , in case (3)
The latter equation holds due to the assumption that - k0  maxf-
2- k0 . Hence, D
. By considering subsets of J  having
cardinality in [1; 3
the result follows for the general case, too. 2
The second result provides a similar formula for the minimal distance between
a regular ternary tree and reduced subtrees containing n  2
9 N scenarios.
Proposition 3.2 (7/9-solution)
Let a regular ternary scenario tree with
Then it holds for each n 2 N with 2
Proof: Similarly as in Proposition 3.1 we obtain
for all
for each subset J of Again we have to show
that there exists an index set J  such that #J and that the lower
bound N n
attained with D J . We consider the index set
I  := fi
_ (- k0
denote the tree consisting of all
I  . Figure 3 illustrates a detail of T r  starting at a
Scenario Reduction Algorithms in Stochastic Programming 13
r
r
r
r
r r
@ @
@ @
@ @
@ @

Figure

3: Detail of the subtree T r
node at level k 0 1 and ending at level k 2. We obtain for the cardinality
of I  and J  that
9 3
9 N and #J
Similarly as in Proposition 3.1 it can be shown that for each j 2 J  , there
exists an index i 2 I  such that it holds that k!
9 - k0 . By considering subsets of J  having cardinality in
9 N ], the result follows for the general case, too. 2
Similar results are available under additional assumptions in case of the Euclidean
norm instead of the maximum norm (see also [6]).
The aim of this section is to report on numerical experience of testing and
comparing the algorithms described in Section 2, namely, backward reduction
of scenario sets, simultaneous backward reduction, fast forward selection. All
algorithms were implemented in C. The test runs were performed on an HP
9000 (780/J280) Compute-Server with 180 MHz frequency and 768 MByte
main memory under HP-UX 10.20, i.e., the same conguration as for the numerical
tests in [4]. We consider the situation where the function c is dened
by c(!; ~
!) and the original discrete probability
measure P is given in scenario tree form. More precisely, we use a test battery
of three binary and ternary scenario trees, respectively. All test trees
are regular and, thus, the results of Section 3 apply. They provide minimal
distances of P to reduced measures supported by n scenarios
when n is not too small.
Example 4.1 (binary scenario tree)
(0:5; 0:6; 0:7; 0:9; 1:1; 1:3; 1:6; 1:9; 2:3; 2:7).

Figure

4 illustrates the original scenario
tree. Proposition 3.1 applies with k
N holds for
each N= 256  n < N .
Example 4.2 (ternary scenario tree)
(0:7; 0:9; 1:2; 1:5; 2:6; 3:3). The tree is shown in Figure 5. Proposition 3.2 applies
with
N holds for each 2N
14 Heitsch, Romisch
Example 4.3 (ternary load scenario tree)
We consider the scenario tree construction in Section 4 of [4] for the weekly
electrical load process of a German power utility (see also [5,8] for a description
of a stochastic power management model and its solution by Lagrangian
relaxation). The original construction is based on an hourly discretization of
the weekly time horizon with branching points at t
and on a piecewise linear interpolation between the t k . The corresponding
mean shifted tree is illustrated in Figure 6. For a moment, we disregard all
non-branching points of the time discretization and consider the corresponding
mean shifted tree. The latter tree is a regular ternary scenario tree with
denotes the standard deviation of the stochastic
load process at time t. Since in this case  t increases with increasing t,
Proposition 3.2 applies with k it holds that D min
N for
Finally, it remains to remark that, due to the piecewise
linear structure of the scenarios and the choice of the maximum norm for
dening c, the minimal distance D min
n does not change when including all
non-branching points.
The original scenario trees of the Examples 4.1{4.3 were reduced to trees containing
n scenarios by using all 3 reduction algorithms. The corresponding
tables contain the relative accuracy and the running time of each algorithm
needed to produce a reduced tree with n scenarios. In addition, the tables
provide the (relative) lower bound (15) and the (relative) minimal distance
n in percent if available. Here, \relative" always means that the corresponding
quantity is divided by the minimal  c -distance of P and one of its
scenarios endowed with unit mass. In particular, the relative accuracy is dened
as the quotient of the  c -distance of the original measure P and the
reduced measure Qn (having n scenarios) and of the  c -distance of P and the
rel
c
denotes the set of scenarios of P and ! i  is dened by
c
c
Our numerical experience shows that all algorithms work reasonably well.
All algorithms reduce 50% of the scenarios of P in an optimal way. As ex-
pected, simultaneous backward reduction and fast forward selection produce
more accurate trees than backward reduction of scenario sets at the expense
of higher running times. Our results also indicate that fast forward selection
is slightly more accurate than simultaneous backward reduction, although
both backward reduction variants are sometimes competitive. Fast forward
selection works much faster than the implementation of forward selection in
Scenario Reduction Algorithms in Stochastic Programming 15
Figure

4: Original binary scenario tree
Number Backward of Simultaneous Fast Lower Minimal
n of Scenario Sets Backward Forward Bound Distance
Scenarios  rel
c
Time  rel
c
Time  rel
c
Time

Table

1: Results of binary scenario tree reduction
Heitsch, Romisch
Figure

5: Original ternary scenario tree
Number Backward of Simultaneous Fast Lower Minimal
n of Scenario Sets Backward Forward Bound Distance
Scenarios  rel
c
Time  rel
c
Time  rel
c
Time

Table

2: Results of ternary scenario tree reduction
Scenario Reduction Algorithms in Stochastic Programming 17
-50050024 48 72 96 120 144 168

Figure

Original load scenario tree
Number Backward of Simultaneous Fast Lower Minimal
n of Scenario Sets Backward Forward Bound Distance
Scenarios  rel
c
Time  rel
c
Time  rel
c
Time

Table

3: Results of load scenario tree reduction
Heitsch, Romisch
[4]. For instance, fast forward selection required 35 seconds to determine a
load scenario subtree (Example 4.3) containing 600 scenarios instead of 8149
seconds reported in [4]. Especially, in the case of deeply reduced trees, fast
forward selection works very fast and accurately.
Furthermore, it turned out that the lower bound is very good (even optimal)
for large n, but extremely pessimistic for small n. Another observation is
that the reduction of half of the scenarios implies only a loss of about 10%
of the relative accuracy. For instance, in case of Example 4.2 it is possible to
determine a subtree containing just 6 out of the originally 729 scenarios that
still carries about 50% of the relative accuracy.
Finally, we take a closer look at the numerical results of the load scenario
tree reduction. In particular, we compare the running times of simultaneous
backward reduction and fast forward selection in this case. Figure 7 displays10300 100 200 300 400 500 600
Time
in
seconds
Number of scenarios
fast forward
simultaneous backward

Figure

7: Running time for reducing the load scenario tree
the running times of both algorithms and shows clearly their opposite algorithmic
strategies. It re
ects the corresponding theoretical complexity results
(Propositions 2.3 and 2.6) and shows that the running time of fast forward
selection is smaller if n  N(approximately). This conrms again that the
forward selection concept is favourable if n is small. Figures 8, 9 and 10
show the reduced load trees with 15 scenarios obtained by all algorithms.
The gures display the scenarios with line width proportional to scenario
probabilities.
Scenario Reduction Algorithms in Stochastic Programming 19
-50050024 48 72 96 120 144 168

Figure

8: Backward reduction / load tree
-50050024 48 72 96 120 144 168

Figure

9: Simultaneous backward reduction / load tree
-50050024 48 72 96 120 144 168

Figure

10: Fast forward selection / load tree
Heitsch, Romisch

Acknowledgement

This research was partially supported by the BMBF project 03-ROM5B3.


Appendix


Proof (Theorem 2.1): Let J  I := be an arbitrary index set.
We set c ij := c(!
programming duality implies for any feasible q
In particular, we consider
It holds
ki and
I . Hence, we obtain
Next, we set u i := min
for each i 2 I . Noting that u
all
for all Hence, we obtain for any feasible q that
and the proof is complete. 2



--R

Introduction to Algorithms




Reduktion von Szenariob
in: Handbooks in Operations Research and Management Science


Probability Metrics and the Stability of Stochastic Models



Theory of Probability and its Applications 28
--TR
Introduction to algorithms
Integer programming
Quantitative Stability in Stochastic Programming
scenario reduction;electrical load;probability metric;scenario tree;stochastic programming