Numerical Evaluation of Approximation Methods in
Stochastic Programming
Christian Küchler and Stefan Vigerske
Humboldt–Universität zu Berlin, Unter den Linden 6,
D–10099 Berlin, Germany,ckuechler@math.hu-berlin.de, stefan@math.hu-berlin.de

Abstract
We study an approach for the evaluation of approximation and solution methods for multistage linear stochastic programs by measuring the performance of theobtained solutions on a set of out-of-sample scenarios. The main point of the approach is to restore the feasibility of solutions to an approximate problem along theout-of-sample scenarios. For this purpose, we consider and compare different feasibility and optimality based projection methods. With this at hand, we study thequality of solutions to different test models based on classical as well as recombining
scenario trees.
Keywords. Multistage Stochastic Programming, Scenario Tree, Out-of-sample Evaluation
AMS subject classification. 90C15, 90C39, 90C59, 49M25

1

Introduction

In general, numerical solution methods for stochastic optimization problems require theunderlying probability measures to have only a finite support. Thus, different techniques
have been developed to approximate random variables or stochastic processes by a limited
number of scenarios or finite scenario trees, respectively. These techniques follow different
principles like random sampling [19], moment matching [8], probability metrics [3, 4, 6, 15],and Quasi Monte-Carlo sampling [14]. Convergence of optimal values and/or solution sets
has been proved for specific techniques and properties of statistical estimates and bounds
have been established (cf., e.g., [19] and the references therein). Stability analysis ofstochastic programs yields further hints how approximations should look like, cf. [5, 11,
13, 16].

1

Numerical Evaluation of Approximation Methods

2

Unfortunately, on the one hand, these theoretical results may require the optimization
problems and underlying random variables to fulfill specific regularity assumptions that
may be hard to verify in some cases of practical interest. On the other hand, quantitative
error bounds and statistical properties are not available for all problem classes. Furthermore, due to the numerical complexity of stochastic programming models, it is sometimes
necessary to use approximations that are too rough to obtain meaningful error bounds orconfidence intervals via asymptotic results.
In such cases, one has to resort to numerical methods to measure the performanceand quality of approximation and solution methods. Since a main task of stochastic
programming is to provide decision strategies that are robust enough to be applicable inreal-world scenarios, it suggests itself to measure the quality of an approximation methodby evaluating the (optimal) solutions obtained from solving the approximate problem.
This can be done, e.g., by evaluating these solutions along out-of-sample scenarios, cf.,
e.g., [10] and [2, 7] for one- and multistage problems, respectively.
In this paper, we study how out-of-sample testing may be used to study the behaviourof approximations to linear multistage stochastic programs (MSP). Thereby, we aim forproblems with many stages, where, due to numerical complexity, the thoroughly construction of out-of-sample strategies as in [7] and the second method of [2] do not apply.
Furthermore, our framework differs since we abstain from a (relatively) complete recourse
assumption. Then, in particular, optimal solutions of an approximate problem are notnecessarily feasible along out-of-sample scenarios. Therefore, the generation of feasible
solutions out of solutions of an approximate problem is an important issue. Furthermore,this question may be of interest whenever one is interested in obtaining practically applicable solutions. For this feasibility restoration we adopt different projection approaches.
Considering MSPs from power scheduling and finance, the proposed feasibility restoration approaches are applied to study the quality of solutions obtained by the decomposition
approach proposed in [12], based on recombining scenario trees, and solutions induced by
non-recombining trees that have been constructed using the Forward Tree Construction
Algorithm of Heitsch and Römisch [4, Algorithm 4.5].
Acknowledgement. We are grateful to Teemu Pennanen for motivating us to this workand to Prof. Werner Römisch for his help and encouragement. This work was supported by
the Bundesministerium für Bildung und Forschung (BMBF) under the grant 03SF0312E,
which is gratefully acknowledged.

Numerical Evaluation of Approximation Methods

2

3

Problem Formulation

On a probability space (Ω, F, P) we consider an Rs -valued discrete time stochastic processξ = (ξ t )t=1,...,T . For t = 1, . . . , T , the vector (ξ 1 , . . . , ξ t ) is denoted by ξ [t] . We study thefollowing linear multistage stochastic program:)(x ∈ Mm , xt ≥ 0,
t = 1, ..., T,(P ) v(ξ) := inf E [ϕ(ξ, x(ξ))] :At,0 xt (ξ [t] ) + At,1 xt−1 (ξ [t−1] ) = ht (ξ t ), t = 2, ..., Twith the cost function
ϕ(ξ, x(ξ)) :=

T
X
­
t=1

®
ct (ξt ), xt (ξ[t] ) .

The set Mm consists of all tuples x = (x1 , . . . , xT ) of Borel-measurable mappings xt :Rs·t → Rm . The costs ct (·) and the right-hand sides ht (·) are affine mappings from Rs to
Rm and Rr , respectively, for t = 1, . . . , T . The assumption of non-random technology andrecourse matrices is due to notational simplicity.
Whenever the optimization problem (P ) is not solvable analytically, a usual approachis to replace the process ξ by a process ξ̄ taking only a finite number of scenarios ξ¯j =(ξ¯tj )t=1,...,T , j ∈ J, with J being some index set. The approximate problem()£¤ x ∈ Mm , xt ≥ 0,
t = 1, ..., T(P̄ ) v(ξ̄) := inf E ϕ(ξ̄, x(ξ̄)) :,At,0 xt (ξ̄ [t] ) + At,1 xt−1 (ξ̄ [t−1] ) = ht (ξ̄ t ), t = 2, ..., T

may then be solved by numerical methods (that may be again approximative). Denotingby x̄ = (x̄t (·))t=1,...,T a solution resulting from the particular approximation and solution
procedure, the optimal value of the approximate problemv(ξ̄) = E[ϕ(ξ̄, x̄(ξ̄))]is often considered as an approximation of v(ξ). However, specific regularity assumptions
on the problem (P ) and the processes ξ and ξ̄ are necessary to ensure certain approximation qualities, cf. [5], [11], and [13]. Indeed, without such conditions ξ̄ may be close to ξin some sense, but passing from (P ) to (P̄ ) may lead to significant changes in the optimalvalue, e.g., by providing arbitrage possibilities, see [11, Example A.4].
Being interested in a good approximation of the unknown value v(ξ), it is thus reasonable rather to evaluate the approximate solution x̄ with regard to the original dataprocess ξ, that is, to consider
E[ϕ(ξ, x̄(ξ))].
Furthermore, whenever the approximate solution x̄ is not feasible (or even not defined)
along the initial process ξ, it may be appropriate to modify x̄ to a feasible strategy x̃.
Then the value(1)

E[ϕ(ξ, x̃(ξ))]

Numerical Evaluation of Approximation Methods

4

provides an upper bound on v(ξ) that can be realized by implementing the strategy x̃.
The value (1) thus appears to be a more reliable approximation of v(ξ) than E[ϕ(ξ, x̄(ξ))].
To evaluate the integral (1), the law of large numbers suggests to draw independent
samples ξ i , i ∈ I, from the distribution of ξ and to consider the out-of-sample value(2)ṽ(ξ̄) :=1 X
ϕ(ξ i , x̃(ξ i )).
|I| i∈I

The value ṽ(ξ̄) can be seen as the real-world performance of the approximative solution′′′x̄(·). Consequently, approximations ξ̄ and ξ̄ of different accuracy or constructed by′different algorithms may be compared by means of their out-of-sample values ṽ(ξ̄ ) and′′ṽ(ξ̄ ). Similarly, approximative solution algorithms can be compared by evaluating theresulting solutions.

3Out-of-sample Evaluation

We consider a solution x̄(·) to the approximate problem (P̄ ) and denote the finite setof scenarios of the approximating process ξ̄ by {ξ¯j : j ∈ J}. Starting from x̄(·), we
aim to construct a strategy x̃(·) that is feasible along a set of out-of-sample scenarios{ξ i , i ∈ I} ⊂ supp Pξ . In order to ensure that x̃(·) is implementable by a non-clairvoyant
decision maker, this feasibility restoration has to be nonanticipative.
To this end, we consider a nonanticipative mapping π : {ξ i : i ∈ I} → {ξ¯j : j ∈ J}
that assigns every out-of-sample scenario ξ i to some scenario of the approximated processthat is close to ξ i , in some sense. We say that π is nonanticipative if it can be writteniis·t
as π(ξ i ) = (π1 (ξ[1]), . . . , πT (ξ[T→] )), where πt are Borel measurable mappings from R
ji¯{ξt : j ∈ J}. Assuming a decision maker who has observed {ξ [t] = ξ[t] } until time t, theiirule π suggests him a scenario (π1 (ξ[1]), . . . , πt (ξ[t])) of the approximate model (and thecorresponding strategy) that is close to his observation. The mapping π may be defined
as a (conditional) projection, see the Appendix for a detailed construction.
The distance between the set of out-of-sample scenarios {ξ i : i ∈ I} and their associated tree scenarios π(ξ i ) can be measured by the term
1 X
dπ (I, J) :=|I| i∈IPTii
t=1 kξt − πt (ξ[t] )k,PTi
t=1 kξt kwhich is the relative euclidean distance between a scenario ξ i and its assigned scenario
π(ξ i ), averaged over all scenarios i ∈ I.
Having related the out-of-sample scenarios ξ i to the approximation scenarios ξ¯j by themapping π, we obtain that x̄(π(·)) ∈ Mm , i.e., x̄(π(ξ)) is nonanticipative w.r.t. the processξ and thus indeed a potential solution to the initial problem (P ). Unfortunately, x̄(π(ξ i ))

Numerical Evaluation of Approximation Methods

5

does not need to be feasible along the scenario ξ i of the initial process ξ, in general. In
order to achieve this feasibility, different projection-based approaches to modify x̄(π(·))
are proposed in the next section.
In the following, we denote the decision x̄(·) along the scenario π(ξ i ) by x̄i and refer
to it as the reference solution. The modification of x̄i along the out-of-sample scenario ξ iis denoted by x̃i .

3.1

Feasibility Restoration

Aiming for a (nonanticipative) solution x̃i that is feasible along the scenario ξ i , we propose
the following straightforward approach. Let x̃i1 := x̄i1 . For t = 2, . . . , T and given x̃it−1 ,we search for a feasible point x̃it that is close to x̄it .
Such a point x̃it may be found by projecting x̄it on the feasible set at timestage t,i.e., on the set {xt ∈ Rm : At,0 xt + At,1 x̃it−1 = ht (ξ i ), xt ≥ 0}. However, in order tocope with possible future infeasibilities in models without relatively complete recourse,
we further restrict the feasible set by incorporating information about future constraints.
More precisely, we consider the value(3)∆it := min kxt − x̄it k∞
xt

s.t. At,0 xt + At,1 x̃it−1 = ht (ξti ),hlow≤ Aτ,0 xτ + Aτ,1 xτ −1 ≤ hupττ ,xτ ≥ 0,τ = t + 1, . . . , T,
τ = t, . . . , T,

being the minimal distance from x̄it onto the (reduced) feasible set at timestage t. The
lowvectors hlow
and hup≤ hτ (ξτi ) ≤ hupττ are chosen such that hττ holds true for all i ∈ I.
The corresponding conditions in (3) are added to avoid decisions x̃it that may lead toup
future infeasibilities, at least to some degree. In particular, we set hlow
τ,j = hτ,j for thosecomponents j of hτ (·) that do not depend on ξ. Observe that this simple approach to
avoid future infeasibilities relies on the assumption of non-random matrices Aτ,0 and Aτ,1 .
However, the approach can be extended, e.g., by demanding the existence of feasible
decisions xτ , τ ≥ t, along all possible future realizations of the process ξ̄.
Basic Restoration
One may think about several techniques for determining a feasible point x̃it based on
previously computed values for x̃i1 , . . ., x̃it−1 . A basic method is to just stay as close asNumerical Evaluation of Approximation Methods

6

possible to the reference solution x̄it and to set x̃it to a solution point of problem (3):x̃it := argmin kxt − x̄it k∞
xt

s.t. At,0 xt + At,1 x̃it−1 = ht (ξti ),hlow≤ Aτ,0 xτ + Aτ,1 xτ −1 ≤ hupττ ,xτ ≥ 0,τ = t + 1, . . . , T,
τ = t, . . . , T.

Myopic Restoration
However, sometimes it may be reasonable to exchange some closeness to x̄it by cost minimality along the out-of-sample scenario ξ i . That is, we allow the decision x̃it to deviate
from the set of closest feasible solutions by a relative fraction εrel ≥ 0 or an absolute value
εabs ≥ 0 in order to minimizing the costs along ξ i . Doing so in a myopic way means tominimize hct (ξ i ), xt i, i.e.,x̃it := argmin hct (ξ i ), xt i + ρt kxt − x̄it k∞
xt

s.t. At,0 xt + At,1 x̃it−1 = ht (ξti ),hlow≤ Aτ,0 xτ + Aτ,1 xτ −1 ≤ hupττ ,xτ ≥ 0,τ = t + 1, . . . , T,
τ = t, . . . , T,

kxt − x̄it k∞ ≤ (1 + εrel )∆it + εabs .
where ρt ≥ 0 is chosen very small, e.g., ρt = 10−4 kct (ξ i )k∞ .
Note that letting ρt → 0 and εrel +εabs → ∞ increases the emphasis on cost minimality
along ξ i by allowing larger deviation from the reference solution x̄it .
Farsighted Restoration
Due to the time-coupling constraints, a decision xt at time t impacts the feasible sets for
future decisions and thus the future costs. These future costs can be taken into accountwithin the feasibility restoration by considering the shadow prices associated to the timecoupling constraints. For this purpose, we consider the dual problem to (P̄ )," T# µ ∈ Mr , µ1 A1,1 ≤ c1 , µT (ξ̄ )AT,0 ≤ cT (ξ̄ ), 

[T ]
T
 X

­
®
¤
£
ht (ξ̄ t ), µt : µt (ξ̄ [t] )At,0 + E µt+1 (ξ̄ [t+1] )At+1,1 |ξ̄ [t] ≤ ct (ξ̄ t ),(D̄) max E
,




t=1
t = 2, .., T − 1.

and denote by µ̄ an optimal solution of (D̄). The shadow price vector corresponding to
the primal decision xt is then equal to µ̄t+1 (ξ̄ [t+1] )At+1,1 . In particular, this value is a
subgradient of the cost-to-go function at time t + 1, that will be defined below. In order7

Numerical Evaluation of Approximation Methods

to maintain the nonanticipativity of the feasibility restoration along the out-of-sample
scenario ξ i , we resort to the expected shadow price vector
¯
¤
£
i
) .
ηt+1 := E µt+1 (ξ̄ [t+1] )At+1,1 ¯ξ̄ [t] = πt (ξ[t]Using this dual information about future costs, we choose x̃it in a farsighted way by settingx̃it := argmin hct (ξ i ) − ηt+1 , xt i + ρt kxt − x̄it k∞
xt

s.t. At,0 xt + At,1 x̃it−1 = ht (ξti ),hlow≤ Aτ,0 xτ + Aτ,1 xτ −1 ≤ hupττ ,xτ ≥ 0,τ = t + 1, . . . , T,
τ = t, . . . , T,

kxt − x̄it k∞ ≤ (1 + εrel )∆it + εabs .
Extensive Restoration
An even more farsighted method is to use not only a single subgradient, but several. Such
subgradients are available, e.g., whenever the approximate problem (P̄ ) was solved by a
Nested-Benders-style algorithm [17]. Thus, for t = 1, . . . , T , let¯£¤
Qt (xt−1 , ξ[t] ) := min hct (ξt ), xt i + E Qt+1 (xt , ξ̄ [t+1] ) ¯ξ̄ [t] = ξ[t]xts.t. At,0 xt + At,1 xt−1 = ht (ξt ),xt ≥ 0,be the cost-to-go function at time t in state ξ[t] (with QT +1 (·, ·) := 0). Having subgradient¯£¤
information πt′ about E Qt+1 (·, ξ̄ [t+1] ) ¯ξ̄ [t] = ξ[t] available in a set of supporting points
x′t , we can set
x̃it := argmin hct (ξ i ), xt i + Θt + ρt kxt − x̄it k∞
xt

s.t. At,0 xt + At,1 x̃it−1 = ht (ξti ),hlow≤ Aτ,0 xτ + Aτ,1 xτ −1 ≤ hupττ ,xτ ≥ 0,¯£¤iΘt ≥ E Qt+1 (x′t , ξ̄ [t+1] ) ¯ξ̄ [t] = πt (ξ[t]) + hπt′ , xt − x′t i,τ = t + 1, . . . , T,
τ = t, . . . , T,
for all x′t ,kxt − x̄it k∞ ≤ (1 + εrel )∆it + εabs .

Remark 1. The latter ‘extensive’ method is related to the first approach in [2], which
has been proposed for MSPs with interstage independence or a weak type of interstage
dependence.
Remark 2. It is also possible to apply a preprocessing step to the introduced restoration
methods similar to the optimal basis prolongation in [1]. That is, having an optimalbasis from the solution of (P̄ ) and a feasible solution x̃it−1 at hand, one can construct aNumerical Evaluation of Approximation Methods

8

corresponding primal solution x̄it that satisfies At,0 x̄it + At,1 x̃it−1 = ht (ξti ), but might violate
the positivity constraint xt ≥ 0. Thus, a feasibility restoration step that uses x̄it insteadof x̃it is applied afterwards. However, within our numerical experiments this preprocessing
affects the out-of-sample evaluation rather adversely.

3.2

Measuring Infeasibility

Without relatively complete recourse, the feasibility restoration might fail if problem (3)is infeasible. However, in some cases it might be possible to relax certain ’soft’ constraints
in order to obtain a feasible solution. Let St be a matrix that indicates the ’soft’ dynamic
constraints at timestage t (of course, the concrete choice of St depends on the considered
model). In order to determine how much these constraints have to be relaxed to make (3)feasible, we solve the auxiliary problem(4)

min i |s̃it |1

xt ,...,xT ,s̃ts.t. At,0 xt + At,1 x̃it−1 + St s̃it = ht (ξti ).
hlow≤ Aτ,0 xτ + Aτ,1 xτ −1 ≤ hupττ ,τ = t + 1, . . . , T,xτ ≥ 0,τ = t, . . . , T.

If problem (4) is feasible, we apply one of the feasibility restoration methods from above
with εrel = εabs = 0 and the right side ht (ξ i ) replaced by ht (ξ i ) − s̃it , to obtain a ‘minimal
infeasible’ solution x̃it−1 . With this solution at hand we can proceed to the next timestage.
If the relaxed problem (4) is feasible for every t = 2, . . . , T , we say that the solution x̃i isweakly infeasible. The minimal weak infeasibility of x̃i along ξ i is then defined as
i

(5)

s̃ :=

T
X

|s̃it |1 ,t=1

where we let s̃it = 0 if (3) was feasible in timestage t.
If even the relaxed problem (4) is infeasible for some t ≥ 2, x̃i is denoted as strongly
infeasible and the feasibility restoration for the out-of-sample scenario ξ i is abandoned.
Piecing this all together, it is reasonable to measure the quality of a solution by thefollowing parameters. By I+ , I− , and I−− we denote the sets of indices i ∈ I such that x̃i
is feasible, weakly infeasible, or strongly infeasible, respectively. Then we consider the
1 X
ϕ(ξ i , x̃(ξ i )),(feasible) out-of-sample value ṽ(ξ̄) :=|I+ | i∈I+

|I− ∪ I−− |
, and
|I|
1 X i
s̃ ,
average (weak) infeasibility s̃ :=|I− | i∈I
infeasibility rate r̃ :=−where s̃i is defined by (5).

Numerical Evaluation of Approximation Methods

4

9

Numerical Examples

We apply the proposed out-of-sample evaluation method for two stochastic programming
models using in both cases classical scenario trees as well as recombining scenario trees.
The classical scenario trees are generated using the Forward Tree Construction Algorithm of Heitsch and Römisch [4, Algorithm 4.5]. The optimal value v(ξ̄) of the scenario
tree based problem is computed by solving the deterministic equivalent with CPLEX [9].
The recombining scenario trees are constructed by a modified version of [4, Alg. 4.5].
Recombination takes place every four time stages and the number of different subtrees
per timeperiod (the time between two recombination points) is varied between two andeight. The tree based problem is solved by an extension of the Nested Benders Decomposition Algorithm as introduced in [12]. Note that the decomposition approach does not
require the solution process to be recombining too. However, the recombining nature of theprocess ξ allows a “dynamic recombination” of solutions due to the re-use of Benders Cuts.
Within this decomposition algorithm the quality of the approximated cost-to-go functions
is controlled by an aggregation parameter γ ∈ [0, 1]. Numerical experiments with various
values of γ have shown that the choice γ = 0.01 yields sufficiently accurate solutions.
For details about the generation of recombining scenario trees and the decomposition
approach, see also [12].

4.1Power Scheduling

Model Description
We consider a power generating system consisting of several thermal units (index set I), a
pumped hydro unit, and a wind power plant. The objective is to find cost-optimal operation levels of the thermal units and hydro units under uncertain production of electricity
from wind. The model makes no claim to depict a real world situation, but is for studying
purposes only.
Let us denote by pi,t the operation level of the thermal unit i ∈ I at time t, by lt thefill level of the water reservoir, and by wt and vt the operation level of the water pumpand turbine, respectively. Deterministic parameters of the problem are operation ranges
for the thermal units pi < p̄i , i ∈ I, the pump w̄ > 0, and the turbines v̄ > 0, the capacityof the water reservoirs ¯l > 0, the fill levels lin and lend of the reservoirs at the beginningand at the end of the considered time horizon, the efficiency of the pump η, the fuel costs
bi , i ∈ I, and the energy demand dt . As stochastic parameter we consider the wind power

10

Numerical Evaluation of Approximation Methods
production κt . The complete model has the form
min E

" T
XX

bi pi,tt=1 i∈I(6)s.t.

l0 = lin ,#lT ≥ lend ,lt = lt−1 − (vt − ηwt ),1|pi,t − pi,t−1 | ≤ (p̄i − pi ),2X
pi,t + (vt − wt ) + κt ≥ dt ,

(7)
(8)
(9)

t = 1, . . . , T,i ∈ I, t = 2, . . . , T,
t = 1, . . . , T,

i∈I

X

(10)

i∈I

pi,t ≤

X
i∈I

p̄i −

1
dt ,
10

t = 1, . . . , T,pi ≤ pi,t ≤ p̄i , 0 ≤ vt ≤ v̄, 0 ≤ wt ≤ w̄, 0 ≤ lt ≤ ¯l,i ∈ I, t = 1, . . . , T.

Constraint (6) models the initial and final fill level of the reservoir, (7) couples the fill levelof the reservoir between successive time stages, (8) bounds the change in the operation of
the thermal units between successive time stages, (9) ensures that the electricity demandis covered, and (10) is a reserve requirement. The model parameters are given by Table 1
in the Appendix.
The time horizon is T = 48 hours. The random nature of future wind power output ismodeled by adding a discrete time Brownian motion to the predicted wind speed curve.
Then 1000 trajectories of wind speed are simulated and transformed into wind power
output using a (piecewise linear) wind park power curve, c.f. Figure 1. The resulting wind
power output scenarios are used for the construction of scenario trees. For the bounds
low
on κt (as required for hlow
and hup≡ 0 and κupt in (3)) we use κtt ≡ maxτ kκτ k∞ fortt = 1, . . . , T .
20
100
17.5

15

80

12.5
60
10

7.5

40

5
20
2.5

10

20

30

40

10

20

30

40

Figure 1: Left side: 1000 simulated trajectories from wind speed in m/s. Right side: The
corresponding wind power output trajectories (percentage of maximal capacity). The
maximal capacity of the wind power park is attained for a wind speed above 15 m/s.
While the capacity of the thermal units is sufficient to cover the maximal load, thisNumerical Evaluation of Approximation Methods

11

model does not possess relatively complete recourse due to the condition (6) on the minimal final fill level of the water storage and the reserve requirements (10). Since a violationof the reserve requirement (10) does not prohibit the power plant to operate inside its operational bounds and such that the demand is satisfied, we have chosen these constraints
as soft constraints (cf. Section 3.1), i.e., weakly infeasible solutions are allowed to violate
the reserve requirement.
Out-of-sample Evaluation
Figure 2 presents the results of the out-of-sample evaluation of classical and two types ofrecombining scenario trees using the farsighted feasibility restoration. As one can observe,
the out-of-sample values ṽ(ξ̄) are higher than the minimal costs v(ξ̄) of the tree based
stochastic programs. While the optimal values v(ξ̄) of the approximate problems do notsignificantly differ for scenario trees having different numbers of nodes, the out-of-sample
values are decreasing with a growing number of nodes in the same manner as the distancedπ between the scenario trees and the set of out-of-sample scenarios. This indicates thatthe quality of the tree based solutions with regard to their usefulness for real-world decision
making improves with increasing accuracy of the scenario tree approximation.
Next, we observe that the out-of-sample values ṽ(ξ̄) computed from solutions of recombining tree based stochastic programs are better than when classical scenario trees on
the same number of nodes are used. This is probably due to the much higher number ofscenarios that can be used within recombining trees. Comparing the out-of-sample values
for recombining scenario trees with two and four different subtrees per timeperiod showsthat a too extensive recombination (i.e., using only two subtrees) worsen the results, eventhough the number of nodes is the same.
Further, we compare the out-of-sample values that are obtained with varying valuesof εrel , i.e., allowing the modified solution to differ more or less from the tree solution infavor of cost minimality. The results for the farsighted feasibility restoration are shownon the left side of Figure 2. The choice εrel = 0 is equivalent to the basic method. As
we can see, the out-of-sample values decrease with increasing values of εrel , i.e., with anincreasing freedom to move solutions towards optimality per timestage. However, we also
observed that the number of out-of-sample scenarios that lead to an infeasible problem(3) in the feasibility restoration phase increases. While for εrel ≤ 0.1 the percentage of
infeasible scenarios r̃ (c.f. Section 3.2) is at most 0.1% for all trees, this value increased to
0.4% for εrel = 0.5. For εrel = 1 the fraction of infeasible scenarios raised to 5%, where theinfeasibilities resulting from recombining scenario tree solutions were always of the weak
type. Both types of infeasibilities could be observed when using classical scenario trees.
We also performed computations with the myopic feasibility restoration and observedthat the out-of-sample values are only slightly worse than with the farsighted approach.

12

Numerical Evaluation of Approximation Methods
Scenario tree constructed with the Forward Tree Construction Algorithm [4]:
0.05
0.04

1.7

0.03

1.67

0.02
0.01

1.64
250

500

1000

2000

5000

250

500

1000

2000

5000

Recombining scenario tree with two subtrees per timeperiod:
0.05
0.04

1.7

0.03
1.67

0.02
0.01

1.64
250

500

1000

2000

5000

250

500

1000

2000

5000

Recombining scenario tree with four subtrees per timeperiod:
0.05
0.04

1.7

0.03
1.67

0.02
0.01

1.64
250

500

1000

2000

5000

250

500

1000

2000

5000

Figure 2: Results for the power scheduling model. Left side: Optimal value v(ξ̄) of thetree based problems (in bold) and corresponding out-of-sample values ṽ(ξ̄) for εabs = 0
and several values of εrel using a series of scenario trees with increasing number of nodes(horizontal axis). The values for εrel are (with decreasing position of the ṽ(ξ̄) curve) 0,0.05, 0.10, 0.25, 0.5, and 1.0, where the farsighted feasibility restoration was used. Rightside: The distance dπ between the scenario trees and the set of out-of-sample scenarios.
Hence, the additional information on (future) costs related to changes in the fill level of
the water reservoir does not change the out-of-sample values considerably.

Numerical Evaluation of Approximation Methods

4.2

13

Swing Option Exercising

Model Description
The following problem of evaluating a swing option basically coincides with Problem 1
considered in [7]. A swing option allows its holder to purchase a total amount of U unitsof energy during the time [1, T ], for a fixed strike price K per unit. The amount of energy
purchased at time t ∈ [1, T ] is denoted by xt and has to lie in some interval [l, u]. Assumingthat the purchased energy is immediately sold on the spot market, the holder of the swing
option is interested in finding a purchase strategy xt , t = 1, . . . , T , that maximizes the
expected accumulated wealth. This problem can be written as follows:
( " T
)
#
X
x ∈ M1 , xt ∈ [l, u], t = 1, . . . , T,
(11)
min E
hK − ξ t , xt i : PT
t=1 xt ≤ U,
t=1

with positive constants l, u, U , and K. The stochastic process ξ = (ξ t )Tt=1 describes
the spot market price per energy unit and is assumed to follow a discrete time geometric
Brownian motion, i.e., ξ 1 = 1 and(12)1ξ t = ξ t−1 exp(εt − σ 2 ),2

t = 2, . . . , T.

Thereby εt , t = 2, . . . , T, are independent, normally distributed random variables withexpectation µ and variance σ 2 . In the following, we assume for the sake of notational
simplicity l = 0, u = 1, and U ∈ N.
In particular, this model allows for an analytic solution whenever the drift µ of the
spot price process ξ is nonnegative. Indeed, using the (negative) payoff function ϕ(ξ t ) :=(K − ξ t )− , problem (11) can be written as
( " T
)
#
X
xt ∈ M1 , xt ∈ [0, 1], t = 1, . . . , T,
(13)
min E
.
xt ϕ(ξ t ) : PT
t=1 xt ≤ U
t=1
If µ ≥ 0, the process ξ is a submartingale. Consequently, due to the negativity, monotonicity, and concavity of ϕ, the process (ϕ(ξ t ))t=1,...,T is a supermartingale. It is thus nosurprise that an early exercise of the the swing option is not optimal. Indeed, one easilyshows that x∗ defined by
0, if t ≤ T − U,∗
(14)
xt (ξ t ) :=
1, if t > T − U and ξ > K,t

is an optimal solution for (11). Observe that the algorithms we have used for the scenario
tree construction do not maintain the martingale property, in general. Hence, the tree
based optimal solutions are not necessarily of type (14).

Numerical Evaluation of Approximation Methods

14

Out-of-sample Evaluation
For the following computations, we used the parameter values T = 52, U = 20, K = 1,
µ = 0, and σ = 0.07. We reformulated the model by adding auxiliary variables ut =Ptτ =1 xt and the (redundant) constraints ut ≤ U , t = 1, . . . , T − 1. Thus, by virtue ofthese modifications, the model has relatively complete recourse. This property and thesmall number of variables and constraints might let this model appear ‘simple’. However,while in the power scheduling model uncertainty appeared only in one component of thepower generating system, the decisions in the swing option model are mainly driven by
the stochastic price process. Thus, one can expect that the quality of the price process
approximations ξ̄ has a major impact on the optimal value of the stochastic program andthe out-of-sample values. Indeed, Figure 3 shows that the relative deviation between thetree based optimal values v(ξ̄) and the out-of-sample values ṽ(ξ̄) is much larger than in
the power scheduling problem, see Figure 2.
For the considered classical scenario trees, the values v(ξ̄) and ṽ(ξ̄) appear to be far
from convergent. This illustrates that the number of scenarios used to approximate theprice process over T = 52 timestages does not lead to a reliable solution yet. In contrast,
scenario trees of similar size have shown to be adequate for the power scheduling model(48 timestages, cf. Section 4.1).
Using recombining scenario trees seems to slightly improve the situation. The optimalvalue v(ξ̄) and out-of-sample values ṽ(ξ̄) are less fluctuating here. Again, we can observethat using only two subtrees per timeperiod is not sufficient and yields a too optimistic
value v(ξ̄), while the use of eight scenario trees shifts both the tree value v(ξ̄) and theout-of-sample value ṽ(ξ̄) closer to the exact value v(ξ).
Next, we compare different feasibility restoration strategies. Since uncertainty appears
only in the objective function coefficients, the scenario-wise solutions x̄i , i ∈ I, are always
feasible for the corresponding scenario ξ i , that is ∆it = 0. Hence, the basic feasibility
restoration method coincides with evaluating the scenario solution x̄i along the costsct (ξti ), t = 1, . . . , T . We compare a basic, a myopic, and a farsighted strategy using
εabs = 1. In difference to the power scheduling model, here the basic approach yields thebest results. For the myopic feasibility restoration, the choice εabs = 1 yields a solutionx̃it that is independent of the scenario tree solution x̄. In particular, the out-of-sample
values based on the myopic approach do not depend on the scenario tree. The farsighted
feasibility restoration with εabs = 1 additionally utilizes in each timestage informationfrom the scenario tree solution in form of shadow prices. As one can observe, this allows
to improve the myopic results, even though the out-of-sample value gets worse with agrowing number of nodes.

15

Numerical Evaluation of Approximation Methods
Scenario tree constructed with the Forward Tree Construction Algorithm [4]:
-2

0.25

-2.5

0.2

-3
0.15
-3.5
0.1
-4
0.05

-4.5
250

500

1000

2000

5000

250

500

1000

2000

5000

Recombining scenario tree with two subtrees per timeperiod:
-2

0.25

-2.5

0.2

-3
0.15
-3.5
0.1
-4
0.05

-4.5
250

500

1000

2000

5000

250

500

1000

2000

5000

Recombining scenario tree with eight subtrees per timeperiod:
-2

0.25

-2.5

0.2

-3
0.15
-3.5
0.1
-4
0.05

-4.5
250

500

1000

2000

5000

250

500

1000

2000

5000

Figure 3: Results for the swing option model. Left side: Optimal value v(ξ̄) of thestochastic program (in bold) and the corresponding out-of-sample values ṽ(ξ̄) for basic(lower solid line) and farsighted (upper solid line) feasibility restoration using a series ofscenario trees with increasing number of nodes (horizontal axis). The myopic value andthe exact value v(ξ) are plotted as the upper and lower dashed lines, respectively. The
parameters εrel = 0 and εabs = 1 were used. Right side: The distance dπ between the
scenario trees and the set of out-of-sample scenarios.

Appendix
For the sake of completeness, we sketch how the mapping π : {ξ i : i ∈ I} → {ξ¯j : j ∈ J},
as introduced in Section 3, may be constructed. Given a scenario tree consisting of theNumerical Evaluation of Approximation Methods

16

scenarios ξ¯j , j ∈ J, a node nt at time t is a subset of J, such that the scenarios ξ¯j , j ∈ nt
are indistinguishable until time t. The set of nodes at time t is denoted by Nt . The setof nodes at time t + 1 succeeding from the node nt is denoted by succ(nt ) ⊂ Nt+1 and wehave ∪nt+1 ∈succ(nt ) nt+1 = nt . The realization of the (tree) process ξ̄ on a node nt is denotedby ξ¯tnt . Recall that there is only one node at time t = 1, i.e. N1 = {n1 } = {J}, and there
are |J| nodes at time horizon T , each of it containing a different singular scenario index.
Now, the mapping π is constructed recursively as follows. Given some out-of-samplei
scenario ξ i , we set π1 (ξ[1]) := ξ¯1n1 . Assume that for some t ∈ {2, . . . , T } we havent−1
i
πt−1 (ξ[t−1]) = ξ¯t−1
for some node nt−1 . Then we defineiπt (ξ[t]) :=argminξ̄ nt :nt ∈succ(nt−1 )kξti − ξ¯tnt k.

parameter

value

b1
b2
b3

fuel cost coalfuel cost gas & steam
fuel cost gas

21
48
154

p̄1
p̄2
p̄3
k̄
v̄

capacity coalcapacity gas & steam
capacity gas
capacity windcapacity hydro turbine

1,000
500
500
1,000
2,000

w̄
η
¯l
lin , lend

capacity hydro pump
pump efficiency
capacity hydro storage
initial/final storage level

2,000
0.75
12,0006,000

Table 1: Parameter of power scheduling model. Maximum demand per hour is 2, 000.

References[1] M.S. Casey and S. Sen. The scenario generation algorithm for multistage stochastic linear
programming. Mathematics of Operations Research, 30:615–631, 2005.
[2] A. Chiralaksanakul and D. Morton. Assessing policy quality in multi-stage, stochastic
programming. Stochastic Programming E-Print Series, 12, 2004.
[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] H. Heitsch and W. Römisch. Scenario tree modeling for multistage stochastic programs.
Mathematical Programming, 118:371–406, 2009.

17

Numerical Evaluation of Approximation Methods[5] H. Heitsch, W. Römisch, and C. Strugarek. Stability of multistage stochastic programs.
SIAM Journal on Optimization, 17:511–525, 2006.
[6] R. Henrion, C. Küchler, and W. Römisch. Scenario reduction in stochastic programmingwith respect to discrepancy distances. Computational Optimization and Applications, 43:67–
93, 2009.
[7] P. Hilli and T. Pennanen. Numerical study of discretizations of multistage stochastic programs. Kybernetika, 44:185–204, 2008.
[8] K. Høyland and S.W. Wallace. Generating scenario trees for multistage decision problems.
Management Science, 47(2):295–307, 2001.
[9] ILOG, Inc. (now IBM). CPLEX 10.0. http://www.ilog.com/products/cplex.
[10] M. Kaut and S.W. Wallace. Evaluation of scenario-generation methods for stochastic programming. Pacific Journal of Optimization, 3(2):257–271, 2007.
[11] C. Küchler. On stability of multistage stochastic programs. SIAM Journal on Optimization,
19(2):952–968, 2008.
[12] C. Küchler and S. Vigerske. Decomposition of multistage stochastic programs with recombining scenario trees. Stochastic Programming E-Print Series, 9, 2007. http://www.speps.
org.
[13] R. Mirkov and G.Ch. Pflug. Tree approximations of dynamic stochastic programs. SIAM
Journal on Optimization, 18(3):1082–1105, 2007.
[14] T. Pennanen. Epi-convergent discretizations of multistage stochastic programs. Mathematics of Operations Research, 30(1):245–256, 2005.
[15] G.Ch. Pflug. Scenario tree generation for multiperiod financial optimization by optimal
discretization. Mathematical Programming, 89:251–271, 2001.
[16] W. Römisch. Stability of Stochastic Programming Problems, chapter 8, pages 483–554.
Volume 10 of Ruszczyński and Shapiro [18], 2003.
[17] A. Ruszczyński. Decomposition Methods, chapter 3, pages 141–221.
Ruszczyński and Shapiro [18], 2003.

Volume 10 of[18] A. Ruszczyński and A. Shapiro, editors. Stochastic Programming, volume 10 of Handbooksin Operations Research and Management Science. Elsevier, Amsterdam, 2003.
[19] A. Shapiro. Monte Carlo Sampling Methods, chapter 6, pages 353–425. Volume 10 ofRuszczyński and Shapiro [18], 2003.