Recent progress in two-stage mixed-integerstochastic programming with applications
to power production planning
Werner Römisch and Stefan Vigerske

Abstract We present recent developments in two-stage mixed-integer stochastic
programming with regard to application in power production planning. In particular,
we review structural properties, stability issues, scenario reduction and decomposition algorithms for two-stage models. Furthermore, we describe an application tostochastic thermal unit commitment.
Key words: stochastic programming, two-stage, mixed-integer, stability, scenario
reduction, decomposition algorithms, unit commitment, discrepancy

1 Introduction
Since its beginnings in the late 1980s mixed-integer stochastic programming has
undergone a considerable development both in theory and computations. We referto the excellent overviews in [LoS03, Sch03, Se05] and to the very comprehensive
bibliography [vdV07].
The aim of this paper is to look at some of the more recent developmentsthat bear further potential for applications to power systems modeling and optimization. First, we mention new results on structures and convex approximations[KSV06, vdV04], on estimating and approximating the underlying probability distribution [ER07, RV08], and as a consequence of the latter on scenario reduction
in two-stage mixed-integer stochastic programs [HKR08, HKR09]. Another line ofwork deals with the consequences of replacing the (traditional) expectation functional in the objective by risk functionals on structural properties and algorithms(see [ER05, ST06] for example). Much work was directed to algorithmic issues and,Werner Römisch
Humboldt University, 10099 Berlin, Germany, e-mail: romisch@math.hu-berlin.de
Stefan Vigerske
Humboldt University, 10099 Berlin, Germany, e-mail: vigerske@math.hu-berlin.de1

2

Werner Römisch and Stefan Vigerske

in particular, to decomposition schemes [AEO03, ATS00, CS99, DR04, EG+07,
LuS04, NS08b, SSV98, SeH05, SeS06], where much is due to the pioneering workof S. Sen and his co-workers.
In the following, we review some of the recent work. We start with a review ofstructural properties, discuss stability issues, methods for scenario reduction, anddecomposition algorithms. As an illustration, we finally discuss an application to
the stochastic unit commitment problem in power production planning.

2 Models and structural properties
Stochastic programs with mixed-integer recourse arise as deterministic equivalentsof linear programs containing a random parameter vector ξ (varying in Ξ ) and beingof the form
min{hc, xi | x ∈ X, T (ξ )x ≥ h(ξ )},where X is a closed subset of Rm , c ∈ Rm , the (technology) matrix T (ξ ) andthe vector h(ξ ) may depend on ξ . Given a realization of ξ , a possible violation of h(ξ ) − T (ξ )x ≤ 0 is compensated by the recourse cost hq1 (ξ ), y1 (ξ )i +hq2 (ξ ), y2 (ξ )i, where the pair (y1 (ξ ), y2 (ξ )) with integral y2 satisfies the constraint W1 y1 + W2 y2 ≤ h(ξ ) − T (ξ )x. Here, the cost coefficients q1 (ξ ) and q2 (ξ )may depend on ξ . The modeling idea consists in adding the expected recourse costE(hq1 (ξ ), y2 (ξ )i + hq2 (ξ ), y2 (ξ )i) to the original cost hc, xi and in minimizing thetotal cost with respect to (y1 , y2 ). This leads to the stochastic program with mixedinteger recourse
¯
¾
½Z
¯
¯
(1)
min
f0 (x, ξ )dP(ξ ) ¯ x ∈ X ,Ξ

where the function f0 is given by
f0 (ξ , x) := hc, xi + Φ(q(ξ ), h(ξ ) − T (ξ )x) ((x, ξ ) ∈ Rm × Ξ ),(2)

Φ is the infimum function of a mixed-integer linear programΦ(u,t) := inf{h u1 , y1 i + hu2 , y2 i | y1 ∈ Rm1 , y2 ∈ Zm2 ,W1 y1 +W2 y2 ≤ t}(3)

for all pairs (u,t) ∈ Rm1 +m2 × Rr , Ξ is a polyhedron in Rs , W1 and W2 are (r, m1 )and (r, m2 )-matrices, respectively, q(ξ ) ∈ Rm1 +m2 , h(ξ ) ∈ Rr , and the (r, m)-matrix
T (ξ ) are affine functions of ξ ∈ Rs , and P is a probability distribution on the set Ξ(shortly P ∈ P(Ξ )). Since the decisions x and y(ξ ) are made before and after therealization of ξ , they are called first and second stage decisions, respectively.
The following conditions are imposed to have the model (1) well-defined:(C1) The matrices W1 and W2 have only rational elements.
(C2) For each pair (x, ξ ) ∈ X × Ξ it holds that h(ξ ) − T (ξ )x ∈ T , where
T := {t ∈ Rr | ∃y = (y1 , y2 ) ∈ Rm1 × Zm2 such that W1 y1 +W2 y2 ≤ t} .

Progress in two-stage mixed-integer stochastic programming3(C3) For each ξ ∈ Ξ the recourse cost q(ξ ) belongs to the dual feasible set¯on¯U := u = (u1 , u2 ) ∈ Rm1 +m2 ¯ ∃z ∈ Rr− such that W1> z = u1 , W2> z = u2 .
(C4) P ∈ P2 (Ξ ), i.e., P ∈ P(Ξ ) and Ξ kξ k2 P(dξ ) < +∞.
Condition (C2) means that a feasible second stage decision always exists (relatively complete recourse). Both (C2) and (C3) imply Φ(u,t) to be finite for all(u,t) ∈ U × T . Clearly, it holds (0, 0) ∈ U × T and Φ(0,t) = 0 for every t ∈ T .
With the convex polyhedral coneR

K := {t ∈ Rr | ∃y1 ∈ Rm1 such that t ≥ W1 y1 } = W1 (Rm1 ) + Rr+
one obtains the representation
T =[(W2 z + K ).

(4)

z∈Zm2

The two extremal cases are (i) W1 has rank r implying K = Rr = T (complete
recourse) and (ii) W1 = 0 (pure integer recourse) leading to K = Rr+ .
In general, the set T is connected (i.e., there exists a polygon connecting two
arbitrary points of T ) and condition (C1) implies that T is closed. If, for eacht ∈ T , Z(t) denotes the set
Z(t) := {z ∈ Zm2 | ∃y1 ∈ Rm1 such that W1 y1 +W2 z ≤ t} ,the representation (4) implies that T may be decomposed into subsets of the form
T (t0 ) := {t ∈ T | Z(t) = Z(t0 )} =\(W2 z + K ) \ (z∈Z(t0 )[(W2 z + K )) (5)z∈Zm2 \Z(t0 )for every t0 ∈ T . In general, the set Z(t0 ) is finite or countable, but condition (C1)
implies that Z(t0 ) in the intersection in (5) may be replaced by a single element of
T and Zm2 \ Z(t0 ) in the union by a finite subset of Zm2 , respectively (see [BG+82,
Lemmas 5.6.1 and 5.6.2]). Hence, if (C1) is satisfied, there exist countably many
elements ti ∈ T and zi j ∈ Zm2 for j belonging to a finite subset Ni of N, i ∈ N, such
that
[
[
T =
T (ti ) with T (ti ) = (ti + K ) \(W2 zi j + K ).
(6)
i∈N

j∈Ni

The sets T (ti ), i ∈ N, are nonempty and connected (even star-shaped cf. [BG+82,
Theorem 5.6.3]), but nonconvex in general (see the illustration in Figure 1). If forsome i ∈ N the set T (ti ) is nonconvex, it can be decomposed into a finite numberof subsets of T (ti ) whose closures are convex polyhedra with facets parallel tosuitable facets of W1 (Rm1 ) or of Rr+ (see Figure 1). By renumbering all such subsets(for every i ∈ N) one obtains countably many subsets B j , j ∈ N, of T which forma partition of T . Since the sets Z(t) of feasible integer decisions do not change if

4Werner Römisch and Stefan Vigerske

W2zi,1W2zi,2W2zi,3

B1

B2

B3

B4

ti

Fig. 1 Illustration of T (ti ) (see (6)) for W1 = 0 and r = 2, i.e., K = R2+ , with Ni = {1, 2, 3} and its
decomposition into the sets B j , j = 1, 2, 3, 4, whose closures are convex polyhedral (rectangular).

t varies in some B j , the function Φ(u, ·) is Lipschitz continuous (with modulus notdepending on j) on B j for every j ∈ N and every fixed u ∈ U .
Now, let (C1)–(C3) be satisfied. Then the function Φ is lower semicontinuousand the function (u,t) 7→ Φ(u,t) from U × T to R has the (convex) polyhedral
continuity regions U × B j , j ∈ N. More precisely, the estimate|Φ(u,t) − Φ(ũ, t˜)| ≤ L(max{1, ktk, kt˜k}ku − ũk + max{1, kuk, kũk}kt − t˜k) (7)holds for all pairs (u,t), (ũ, t˜) ∈ U × B j and some constant L > 0. For proofs andfurther details the interested reader is referred to [BG+82, Chapter 5.6].
Next, we consider the integrandf0 (x, ξ ) = hc, xi + Φ(q(ξ ), h(ξ ) − T (ξ )x)
for all pairs (x, ξ ) ∈ X × Ξ and study the continuity properties and growth behaviorof f0 (x, ·) on Ξ for fixed x ∈ X. The properties of Φ imply that, for every x ∈ X,
there exists a partition {Ξx, j } j∈N of Ξ given by
Ξx, j = {ξ ∈ Ξ | h(ξ ) − T (ξ )x ∈ B j }( j ∈ N).

(8)

Furthermore, the function f0 (x, ·) (on Ξ ) satisfies the properties| f0 (x, ξ ) − f0 (x, ξ˜ )| ≤ L̂ max{1, kξ k, kξ˜ k}kξ − ξ˜ k
| f0 (x, ξ )| ≤ C max{1, kxk} max{1, kξ k }2(x ∈ X, ξ , ξ˜ ∈ Ξx, j ), (9)(x ∈ X, ξ ∈ Ξ ),(10)with some positive constants L̂ and C. Due to (10), condition (C4) implies the existence of the integral in (1). We note that f0 (x, ·) is globally Lipschitz continuouson Ξx, j if the recourse cost q(ξ ) does not depend on ξ . It is even globally Lipschitz continuous on Ξ if only q(ξ ) depends on ξ . In both cases | f0 (x, ·)| grows onlyProgress in two-stage mixed-integer stochastic programming

5

linearly with kξ k and a finite first order moment of P, i.e., P ∈ P1 (Ξ ) (instead of(C4)), implies the existence of the integral.
Since the objective function of (1) is lower semicontinuous if the conditions(C1)–(C4) are satisfied, solutions of (1) exist if X is closed and bounded. If theprobability distribution P has a density, the objective function of (1) is continuous,
but nonconvex in general. If the support of P is finite, the objective function is piecewise continuous with a finite number of polyhedral continuity regions. The latter isillustrated by Fig. 2, which shows the expected recourse function
x 7→Z

Φ(q, h(ξ ) − T x)dP(ξ )(x ∈ [0, 5]2 )Ξwith r = s = 2, h(ξ ) = ξ , m1 = 0, W1 = 0, m2 = 4, q = (−16, −19, −23, −28)> , thematrices
µ
¶
µ
¶
2345
2/3 1/3
W2 =and T =,
6132
1/3 2/3and binary restrictions on the second stage variables as in [SSV98], but with a uniform probability distribution P having a smaller finite support than in [SSV98],
namely, supp (P) = {5, 10, 15}2 .

4
2
0
-20

-30

-40

-50

4

2

0

Fig. 2 Illustration of an expected recourse function with pure 0 − 1 recourse, random right-handside and discrete uniform probability distribution.

3 Stability
In this section, we review stability results for mixed-integer two-stage stochastic
programs (1), i.e., results on the dependence of their solutions and optimal values
on the underlying probability distribution P. Such results also provide information

6

Werner Römisch and Stefan Vigerske

how an underlying probability distribution should be approximated such that approximate solutions and optimal values get close to the original ones.
In this context it is well known that the behavior of the (first stage) solution set
¯Z
¾
½
¯
¯
f0 (x, ξ )P(dξ ) = v(P)S(P) := x ∈ X ¯
Ξwith respect to changes of P requires knowledge on the growth of the objectivefunctionZ
x 7→ FP (x) := E( f0 (x, ξ )) =f0 (x, ξ )P(dξ )Ξnear S(P). Here, v(P) denotes the infimum of the objective function or optimalvalue, i.e.,
¯
¾
½Z
¯
¯
v(P) := inf
f0 (x, ξ )P(dξ ) ¯ x ∈ X .
Ξ
However, the growth behavior of FP depends essentially on properties of the underlying probability distribution P. The situation is different for optimal values v(P).
Their behavior with respect to changes of P depends essentially on structural properties of the function f0 , which are well studied (cf. Section 2).
It is shown in [RV08] that the following distances of probability distributions areimportant for mixed-integer two-stage stochastic programs:
¯¯
¾
½¯Z
Z
¯¯
¯
ζ`,B (P, Q) := sup ¯¯ f (ξ )P(dξ ) − f (ξ )Q(dξ )¯¯ ¯¯ f ∈ F` (Ξ ), B ∈ B , (11)B

B

where ` ∈ {1, 2} and B is a set of convex polyhedra, which contains the closures of
Ξx, j , j ∈ N, x ∈ X (see (8)), and F` (Ξ ) contains all functions f : Ξ → R such that| f (ξ )| ≤ max{1, kξ k` } and | f (ξ ) − f (ξ˜ )| ≤ max{1, kξ k`−1 , kξ˜ k`−1 }kξ − ξ˜ kholds for all ξ , ξ˜ ∈ Ξ . While the set F` (Ξ ) of functions has its origin in property(9) of the integrand f0 , but depends on the specific structure of the second stage
program only with respect to ` ∈ {1, 2}, the class B of convex polyhedra strongly
depends on that structure.
If the conditions (C1)–(C4) are satisfied and X is closed and bounded, there exists
a constant L > 0 such that the estimate|v(P) − v(Q)| ≤ LϕP (ζ`,B (P, Q))(12)holds for every Q ∈ P` (Ξ ) with ` ∈ {1, 2} and ` = 2 if ξ enters q(ξ ) and, in addition, h(ξ ) or T (ξ ). Here, the function ϕP is defined by ϕP (0) = 0 and
½
¾
Z
r+1
`
ϕP (t) := inf R t +kξ k P(dξ )(t > 0).
R≥1{ξ ∈Ξ | kξ k>R}

Progress in two-stage mixed-integer stochastic programming

7

The function characterizes thetail behavior of P and is continuous at t = 0. If P has
Ra finite pth moment, i.e., if Ξ kξ k p P(dξ ) < +∞, for some p > `, the estimatep−`(t ≥ 0)ϕP (t) ≤ Ct p+r−1

is valid for some constant C > 0 and if Ξ is bounded, the estimate (12) simplifies to|v(P) − v(Q)| ≤ Lζ`,B (P, Q).
If the set Ξ ⊂ Rs belongs to B, we obtain from (11) by choosing B := Ξ and f ≡ 1,
respectively,max{ζ` (P, Q), αB (P, Q)} ≤ ζ`,B (P, Q)(13)
for all P, Q ∈ P` (Ξ ). Here, ζ` and αB denote the `th order Fortet-Mourier metric(see [Ra91, Section 5.1]) and the polyhedral discrepancy
¯¯
½ ¯Z
¾
Z
¯
¯¯
ζ` (P, Q) := sup ¯¯ f (ξ )P(dξ ) −f (ξ )Q(dξ )¯¯ ¯¯ f ∈ F` (Ξ ) , (14)ΞΞαB (P, Q) := sup |P(B) − Q(B)|,

(15)

B∈B

respectively. Hence, convergence of probability distributions with respect to ζ`,B
implies their weak convergence, convergence of `th order absolute moments, andconvergence with respect to the polyhedral discrepancy αB . For bounded Ξ thetechnique in [Sch96, Proposition 3.1] can be employed to obtain1

ζ`,B (P, Q) ≤ Cs αB (P, Q) s+1(P, Q ∈ P(Ξ ))(16)

for some constant Cs > 0. In view of (13) and (16), the metric ζ`,B is stronger thanαB in general, but in case of bounded Ξ both distances metrize the same convergence on P(Ξ ).
For more specific models (1), improvements of the stability estimate (12) may
be obtained by exploiting specific recourse structures, i.e., by using additional information on the shape of the sets B j , j ∈ N, and on the behavior of the function Φon these sets. This may lead to stability estimates with respect to distances that are(much) weaker than ζ`,B . For example, if W1 = 0, Ξ is rectangular, T is fixed andsome components of h(·) coincide with some of the components of ξ , the closuresof Ξx, j , x ∈ X, j ∈ N, are rectangular subsets of Ξ , i.e., belong to©
ªBrect := I1 × I2 × · · · × Is | 0/ 6= I j is a closed interval in R, j = 1, . . . , s(17)
and the stability estimate (12) is valid with respect to ζ`,Brect . As shown in [HKR09]
convergence of a sequence of probability distributions with respect to ζ`,Brect isequivalent to convergence with respect to both ζ` and αBrect . If, in addition to theprevious assumptions, q is fixed and Ξ is bounded, the estimate (12) is valid withrespect to the rectangular discrepancy αBrect (see also [Sch96, Section 3]).

8

Werner Römisch and Stefan Vigerske

4 Scenario reduction
A well known approach for solving two-stage stochastic programs computationallyconsists in replacing the original probability distribution by a discrete distribution
based on a finite number of scenarios. Let P be such a discrete distribution withscenarios ξ i and probabilities pi , i = 1, . . . , N. The corresponding stochastic programming model is of the form¯¯ W1 yi +W2 yi ≤ h(ξ i ) − T (ξ i )x 

N
1
2
¯
min hc, xi + ∑ pi (hq1 (ξ i ), yi1 i + hq2 (ξ i ), yi2 i) ¯¯ yi1 ∈ Rm1 , yi2 ∈ Zm2 , i = 1, . . . , N, .
¯x ∈ X
i=1
It may turn out that the computing times for solving the resulting mixed-integer linear programs are not acceptable. In such a case one might wish to reduce the numberof scenarios entering the stochastic program. In [DGR03, HR07] a stability-based
approach for scenario reduction in two-stage models without integrality requirements is developed. This approach suggests to look at stability results for optimalvalues and to use the corresponding distance of probability distributions for determining discrete distributions based on a smaller and prescribed number of scenarios
as best approximations of P. According to the stability estimate (12) in Section 3,
the distances ζ`,B or αB (if Ξ is bounded) appear as the right choice, where B is aset of convex polyhedra that depend on the structure of the stochastic program (1).
In [HKR08] the scenario reduction approach is elaborated for αB and a relevant
set B of convex polyhedra. The numerical results show that the complexity of scenario reduction algorithms increases if B gets more involved. To avoid this effect,
the distance ζ`,Brect or, equivalently,
dλ (P, Q) := λ αBrect (P, Q) + (1 − λ )ζ` (P, Q)(18)

for some λ ∈ (0, 1) is considered in this section (and in [HKR09b]).
Let QJ denote a probability distribution whose support supp(QJ ) contains thefollowing subset of {ξ 1 , . . . , ξ N }:©
ªsupp(QJ ) = ξ i | i ∈ {1, . . . , N} \ J
and J ⊂ {1, . . . , N}.
Let qi (i 6∈ J) denote the probability of scenario ξ i of QJ . Now, the aim is to determine QJ such that the distance dλ (P, QJ ) is minimal, i.e., for arbitrary subsets J of
{1, . . . , N}, we are interested in
¯
(
)
¯
¯
DJ := min dλ (P, QJ ) ¯ qi ≥ 0, i ∈/ J, ∑ qi = 1 .
(19)
¯
i∈J
/
In the following, we show that DJ can be computed as optimal value of a linear
program. To this end, we assume without loss of generality that J = {n + 1, . . . , N},i.e., supp(QJ ) = {ξ 1 , . . . , ξ n } for some 1 ≤ n < N. We consider the system of indexsets

Progress in two-stage mixed-integer stochastic programming

9

IBrect := {I(B) := {i ∈ {1, . . . , N} | ξ i ∈ B} | B ∈ Brect }and obtain the following representation of the rectangular discrepancy
¯
¯
¯
¯
¯
¯
q
αBrect (P, QJ ) = sup |P(B) − QJ (B)| = max ¯∑ pi −¯ (20)j∑¯I∈IBrect ¯i∈I
B∈Brectj∈I∩{1,...,n}¾½ ¯¯ − ∑ j∈I∩{1,...,n} q j ≤ tα − ∑i∈I pi , I ∈ IBrect(21)= min tα ¯¯∑ j∈I∩{1,...,n} q j ≤ tα + ∑i∈I pi , I ∈ IBrect
Since the set IBrect may be too large to solve the linear program (21) numerically,
we consider the system of reduced index sets
IB∗ rect := {I(B) ∩ {1, . . . , n} | B ∈ Brect }and the quantities¯)¯¯γ := max ∑ pi ¯ I ∈ IBrect , I ∩ {1, . . . , n} = I ∗
¯
i∈I
¯
(
)
¯
¯
∗
γI ∗ := min ∑ pi ¯ I ∈ IBrect , I ∩ {1, . . . , n} = I
¯
i∈I
I∗

(

for every I ∗ ∈ IB∗ rect . Since any such index set I ∗ corresponds to some left-hand side∗of the inequalities in (21), γ I and γI ∗ correspond to the smallest right-hand sides in(21). Hence, the rectangular discrepancy may be rewritten as½ ¯¾¯ − ∑ j∈I ∗ q j ≤ tα − γ I ∗ , I ∗ ∈ I ∗B¯rect
αBrect (P, QJ ) = min tα ¯.
(22)∑ j∈I ∗ q j ≤ tα + γI ∗ , I ∗ ∈ IB∗ rect
Since the number of elements of IB∗ rect is at most 2n (compared to 2N in IBrect ),passing from (21) to (22) indeed drastically reduces the maximum number of inequalities and may make the linear program (22) numerically tractable.
Due to duality arguments, the Fortet-Mourier distance ζell (P, QJ ) (see (14)) allows the representation as linear program (cf. [HR07])
¯
(
)
¯
N n
N
η
=
q
,
j
=
1,
.
.
.
,
n
∑
¯
i,
j
j
ζ` (P, QJ ) = inf ∑ ∑ ηi j ĉ` (ξ i , ξ j ) ¯ ηi j ≥ 0, ni=1
∑ j=1 ηi, j = pi , i = 1, . . . , N¯
i=1 j=1
where c` (ξ , ξ˜ ) := max{1, kξ k`−1 , kξ˜ k`−1 }kξ − ξ˜ k for all ξ , ξ˜ ∈ Ξ = {ξ 1 , . . . , ξ N }and ĉ` denotes the reduced costs
¯
)
(
¯
K
¯
i
i
i
i
ĉ` (ξ , ξ˜ ) := inf ∑ c` (ξ k−1 , ξ k ) ¯ K ∈ N, ik ∈ {1, . . . , N}, ξ 0 = ξ , ξ K = ξ˜ .
¯
k=1

10

Werner Römisch and Stefan Vigerske

Hence, extending the representation (22) of αBrect we obtain the following linear
program for determining DJ and the probabilities q j , j = 1, . . . , n, of the discrete
reduced distribution QJ ,¯¯ tα ,tζ ≥ 0, q j ≥ 0, ∑nj=1 q j = 1, 



¯



¯ ηi j ≥ 0, i = 1, . . . , N, j = 1, . . . , n, 



¯





¯
N
n
i
j


t
≥
ĉ
(ξ
,
ξ
)η
,


¯
∑
∑
ij
ζ
i=1 j=1 `¯n¯DJ = min λ tα + (1 − λ )tζ ¯ ∑ j=1 ηi j = pi , i = 1, . . . , N,(23)¯ ∑N η = q , j = 1, . . . , n,




j
¯ i=1 i j




¯




I∗ , I ∗ ∈ I ∗
¯


−
q
≤
t
−
γ
∗
∑


j
α
j∈I
Brect ¯¯ ∑ j∈I ∗ q j ≤ tα + γI ∗ , I ∗ ∈ I ∗ 
BrectWhile the linear program (23) can be solved efficiently by available software, the∗determination of the index set IB∗ rect and the coefficients γ I , γI ∗ is more intricate.
∗
It is shown in [HKR08, HKR09b, Section 3] that the parameters IB∗ rect and γ I ,γI ∗ can be determined by studying the set R of supporting rectangles. A rectangle B
in Brect is called supporting if each of its facets contains an element of {ξ1 , . . . , ξn }in its relative interior (see also Fig. 3). Based on R the following representations are

Fig. 3 Non-supporting rectangle (left) and supporting rectangle (right). The dots represent theremaining scenarios ξ 1 , . . . , ξ n for s = 2.

valid according to [HKR08, Prop. 1 and 2]:[ ©
¯ªIB∗ rect =I ∗ ⊆ {1, . . . , n} ¯ ∪ j∈I ∗ {ξ j } = {ξ 1 , . . . , ξ n } ∩ int B
B∈R
I∗

¯
©
ª
= max P(int B) ¯ B ∈ R, ∪ j∈I ∗ {ξ j } = {ξ 1 , . . . , ξ n } ∩ int B¯½¾¯
γI ∗ = ∑ pi where I := i ∈ {1, . . . , N} ¯¯ min∗ ξl j ≤ ξli ≤ max∗ ξl j , l = 1, . . . , s
j∈I
j∈I

γ

i∈I

for every I ∗ ∈ IB∗ rect . Here, int B denotes the interior of the set B.

Progress in two-stage mixed-integer stochastic programming

11

An algorithm is developed in [HKR09b] that constructs recursively l-dimensional
supporting rectangles for l = 1, . . . , s. Computational experiments show that its¡ run¢sning time grows linearly with N, but depends on n and s via the expression n+1.
2Hence, while N may be large, only moderately sized values of n given s are realistic.
Since an algorithm for computing DJ is now available, we finally look at determining a scenario index set J ⊂ {1, . . . , N} with cardinality # J = n such that DJ is
minimal, i.e., at solving the combinatorial optimization problem
min{DJ | J ⊂ {1, . . . , N}, # J = n}(24)

which is known as n-median problem and as NP-hard. One possibility is to reformulate (24) as mixed-integer linear program and to solve it by standard software.
Since, however, approximate solutions of (24) are sufficient, heuristic algorithms
like forward selection are of interest, where uk is determined in its kth step such thatit solves the minimization problem¯no¯min DJ [k−1] \{u} ¯ u ∈ J [k−1] ,where J [0] = {1, . . . , N}, J [k] := J [k−1] \ {uk } (k = 1, . . . , n) and J [n] := {1, . . . , N} \{u1 , . . . , un } serves as approximate solution to (24). Recalling that the complexity of
1

0.8

0.6

0.4

0.2

0
0

0.2

0.4

0.6

0.8

1

Fig. 4 N = 1 000 samples ξ i from the uniform distribution on [0, 1]2 and n = 25 points ξ uk , k =1, . . . , n, obtained via the first 25 elements zk , k = 1, . . . , n, of the Halton sequence (in the bases 2
and 3) (see [Ni92, p. 29]). The probabilities qk of ξ uk , k = 1, . . . , n, are computed for the distancedλ with λ = 1 (gray balls) and λ = 0.9 (black circles) by solving (23). The diameters of the circles
are proportional to the probabilities qk , k = 1, . . . , 25.

12Werner Römisch and Stefan Vigerske¡ ¢sevaluating DJ [k−1] \{u} for some u ∈ J [k−1] is proportional to k+1
shows that even2
the forward selection algorithm is expensive.
Hence, heuristics for solving (24) become important that require only a low number of DJ evaluations. For example, if P is a probability distribution on [0, 1]s withindependent marginal distributions Pj , j = 1, . . . , s, such a heuristic can be based on
Quasi-Monte Carlo methods (cf. [Ni92]). The latter provide sequences of equidistributed points in [0, 1]s that approximate the uniform distribution on the unit cube
[0, 1]s . Now, let n Quasi-Monte Carlo points zk = (zk1 , . . . , zks ) ∈ [0, 1]s , k = 1, . . . , n,be given. Then we determine³´
yk := F1−1 (zk1 ), . . . , Fs−1 (zks )(k = 1, . . . , n),where Fj is the (one-dimensional) distribution function of Pj , i.e.,
Fj (z) = Pj ((−∞, z]) =

N

∑i

pi

(z ∈ R)i=1, ξ j ≤z

and Fj−1 (t) := inf{z ∈ R | Fj (z) ≥ t} (t ∈ [0, 1]) its inverse ( j = 1, . . . , s). Finally, we
determine uk as solution ofmin kξ u − yk ku∈J [k−1]

and set again J [k] := J [k−1] \ {uk } for k = 1, . . . , n, where J [0] = {1, . . . , N}. Figure
4 illustrates the results of such a Quasi-Monte Carlo based heuristic and Figure 5
shows the discrepancy αBrect for different n and the running times of the QuasiMonte Carlo based heuristic.
discrepancy

time
600

0.8
400

0.6
0.4

200
0.2
0

10

20

30

40

50

Fig. 5 Distance αBrect between P (with N = 1000) and equidistributed QMC-points (dashed),
QMC-points, whose probabilities are adjusted according to (23) (bold), and running times of the
QMC-based heuristic (in seconds).

Progress in two-stage mixed-integer stochastic programming

13

5 Decomposition algorithms
When the size of an optimization problem becomes intractable for standard solution
approaches, a decomposition into small tractable subproblems by relaxing certain
coupling constraints is often a possible resort. The task of the decomposition algorithm is then to coordinate the search in the subproblems in a way that their solutions
can be combined into one that is feasible for the overall problem and has a “good”
objective function value. Often, the algorithm also provides a certified lower boundon the optimal value which allows to evaluate the quality of a found solution. Since,on the one hand, mixed-integer stochastic programs easily reach a size that is intractable for standard solution approaches, but, on the other hand, are also very structured, many decomposition algorithms have been developed [LoS03, Sch03, Se05].
In the following, we discuss some of them in more detail.
Let us assume that the set of first stage feasible solutions X is given in the form
X = {x ∈ Zm0 × Rm−m0 | Ax ≤ b},
where m0 denotes the number of first stage variables with integrality restrictions, A
is a (r0 , m)-matrix, and b ∈ Rr0 . Further, we denote byX̄ := {x ∈ Rm | Ax ≤ b}
the linear relaxation of X. We recall the value function (3),Φ(u,t) = inf{hu1 , y1 i + hu2 , y2 i | y1 ∈ Rm1 , y2 ∈ Zm2 ,W1 y1 +W2 y2 ≤ t},
and define the expected recourse function of model (1) byΨ (x) :=Z

Φ(q(ξ ), h(ξ ) − T (ξ )x)P(dξ )(x ∈ X̄).

Ξ

For continuous (m2 = 0) stochastic programs, the Benders decomposition is anestablished method [VSW69, Bi85]. It decomposes the decision on the first stage
from the recourse decisions on the second stage by replacing the value function
Φ(u,t) in (1) by an explicit approximation based on supporting hyperplanes. Unfortunately, Benders Decomposition relies heavily on the convexity of the valuefunction t 7→ Φ(u,t). Thus, in the view of Section 2, it cannot be directly applied tothe case where discrete variables are present.
However, there are several approaches to overcome this difficulty. One of the firstis the Integer L-shaped method [LaL93], which assumes that the first stage problem
involves only binary variables. This property is exploited to derive linear inequalities
that approximate the value function Φ(u,t) pointwise. While the algorithm makes
only moderate assumptions on the second stage problem, its main drawback is the
weak approximation of the value function due to lacking first order information
about the value function. Thus, the algorithm might enumerate all feasible first stage
solutions in order to find an optimal solution.

14

Werner Römisch and Stefan Vigerske

A cutting-plane algorithm is proposed in [CT97]. Here, the deterministic equivalent of (1) is solved by improving its linear relaxation with lift-and-project cuts.
Decomposition arises here in two ways. First, the linear relaxation (including additional cuts) is solved by Benders Decomposition. Second, lift-and-project cuts arederived scenariowise. Further, in case of a fixed technology matrix T (·) ≡ T , cutcoefficients that have been computed for one scenario can also be reused to derive
cuts for other scenarios. This algorithm can be seen as a predecessor of the dual
decomposition approach presented in [SeH05]. While the cuts in [CT97] include
variables from both stages, [SeH05] extends the Benders Decomposition approachto the mixed-integer case by sequentially convexifying the value function Φ(u,t). It
is discussed in detail in Section 5.2.
In [KSV06] it is observed, that even though the value function Φ(u,t) might
be nonconvex and difficult to handle, under some assumptions on the distribution of ξ , the expected recourse function Ψ (x) can be convex. Starting with simple integer recourse models and then extending to more general classes of problems, techniques to compute tight convex approximations of the expected recourse
function by perturbing the distribution of ξ are developed in a series of papers[KSV06, vdV04, vdV05]. We sketch this approach in more detail in Section 5.1.
In the case that the second stage problem is purely integer (m1 = 0), the valuefunction Φ(u,t) has the nice property to be constant on polyhedral subsets of U ×
T . Therefore, in case of a finite distribution, also the expected recourse function
Ψ (x) is constant on polyhedral subsets of X̄. This property allows to reduce theset X to a finite set of solution candidates that can be enumerated [SSV98]. Since
the expected recourse function Ψ (x) has to be evaluated for each candidate, many
similar integer programs have to be solved. In [SSV98] a Gröbner basis for thesecond stage problem is computed once in advance (which is expensive) and thenused for evaluation of Ψ (x) for every candidate x (which is then cheap).
Another approach based on enumerating the sets where Ψ (x) is constant is presented in [ATS00]. Instead of a complete enumeration, here a branch-and-boundalgorithm is applied to the first stage problem to enumerate the regions of constantΨ (x) implicitely. Branching is thereby performed along lines of discontinuity ofΨ (x), thereby reducing its discontinuity in generated subproblems.
While all approaches discussed so far explore the structure of the value or expected recourse function in some way, Lagrange decomposition is a class of algorithms where decomposition is achieved by relaxation of problem constraints. By
moving certain coupling restrictions from the set of constraints into the objectivefunction as penalty term, the problem decomposes into a set of subproblems, eachof them often much easier to handle than the original problem. This relaxed problem
then yields a lower bound onto the original optimal value, which is further improvedby optimization of the penalty parameters. Since, in general, a solution of the relaxed problem violates the coupling constraints, heuristics and branch-and-boundapproaches are applied to obtain good feasible solutions of the original problem.
While there are several alternatives to choose a set of coupling constraints for relaxation, each one providing a lower bound of different quality [DR04], in general, scenario and geographical decomposition are the preferred strategies [CS99, NR00]. In

Progress in two-stage mixed-integer stochastic programming

15

scenario decomposition, nonanticipativity constraints are relaxed, so that the problem decomposes into one deterministic subproblem for each scenario. We discuss
this approach in more detail in Section 5.3. In geographical decomposition, modelspecific constraints are relaxed, which leads to one subproblem for each component of the model. Even though each subproblem then corresponds to a stochastic
program itself, its structure often allows to develop specialized algorithms to solve
them very efficiently. Similarly, the modelers knowledge can be explored to makesolutions from the relaxed problem feasible for the original problem. Geographicaldecomposition is demonstrated for a unit commitment problem in Section 6.

5.1 Convexification of the expected recourse function
In a simple integer recourse model, the second stage variables are purely integer(m1 = 0) and are partitioned into two sets y+ , y− ∈ Zs+ with 2s = m2 . The cost-vectorq(ξ ) ≡ (q+ , q− ) and technology matrix T (ξ ) ≡ T are fixed, r = 2s, h(ξ ) ≡ ξ , andthe value function takes the form¯¯y+ ≥ ξ − T x,¯y− ≥ −(ξ − T x) .
Φ(q(ξ ), h(ξ ) − T (ξ )x) = inf hq+ , y+ i + hq− , y− i ¯¯¯ y+ , y− ∈ Zs+
The simple structure of the value function allows to write the expected recourse
function in a separable form,s+−−Ψ (x) = ∑ q+i E[dξi − (T x)i e ] + qi E[bξi − (T x)i c ],i=1

where dαe denotes the smallest integer that is at least α, bαc the largest integer thatis at most α, α + = max(0, α), and α − = min(0, α). Thus, it is sufficient to consider
one-dimensional functions of the form Q(z) := q+ Eζ [dζ − ze+ ] + q− Eζ [bζ − zc− ](with ζ a random variable).
In [KSV06], convex approximations of Q(z) are derived from a piecewise linear
function in the points (z, Q(z)), z ∈ α + Z, where α ∈ [0, 1) is a parameter. Further,if ζ has a continuous distribution, then the approximation of Q(z) can be realized as
expected recourse function of a continuous simple recourse model,
Qα (z) = q+ Eζα [(ζα − z)+ ] + q− Eζα [(ζα − z)− ] +q+ q−
q+ + q−,where ζα is a discrete random variable with support in α + Z [KSV06].
The results in [KSV06] are extended to derive convex approximations of the
expected recourse function for models of the form (1), where m1 = 0, h(ξ ) ≡ ξ , andq(ξ ) ≡ q and T (ξ ) ≡ T are fixed [vdV04]. Further, the parameter α can be chosensuch that the derived convex approximation underestimates the original expected16

Werner Römisch and Stefan Vigerske

recourse function. Since this convex underestimator is at least as good as an LPbased underestimator (obtained by relaxing the integrality condition on y) and evenyield the convex hull of Ψ (x) in the case that T is unimodular, it can be utilized toderive lower bounds in a Branch-and-Bound search for a solution of (1).
Another extension of the methodology from [KSV06] considers mixed-integer
recourse models where r = 1 and the value function is semi-periodic, c.f. [vdV05].

5.2 Convexification of the value function
From now on we assume that the random vector ξ has only finitely many outcomesξ i with probability pi > 0, i = 1, . . . , N. Thus, we can write the expected recourse
function asNΨ (x) = ∑ pi Φ(q(ξ i ), h(ξ i ) − T (ξ i )x)(x ∈ X̄).

i=1

As discussed in Section 2, the nonconvexity of the function Φ(u,t) forbids a representation by supporting hyperplanes as used in a Benders decomposition. However, while in the continuous case (m2 = 0) the hyperplanes are derived from dual
feasible solutions of the second stage problem, it is conceptually possible to carryover these ideas to the mixed-integer case by introducing (possibly nonlinear) dual
price functions [TiW81]. Indeed, Chvátal and Gomory functions are sufficiently
large classes of dual price functions that allow to approximate the value function
Φ(u,t) [BJ82]. These dual functions can be obtained from a solution of (3) with abranch-and-bound or Gomory cutting plane algorithm [Wo81]. In [CT98] this approach is used to carry over the Benders decomposition algorithm for two-stage
linear stochastic programs to the mixed-integer linear case by replacing the hyperplane approximation of the expected recourse function by an approximation basedon dual price functions. While [CT98] does not discuss how the master problem
with its (nonsmooth and nonconvex) dual price functions can be solved, the seriesof papers [SeH05, NS08b, SeS06, NS05, NS08a] show that a careful constructionof dual price functions combined with a convexification step based on disjunctive
programming allows to implement an efficient Benders decomposition for mixedinteger two-stage stochastic programs.
We consider the following master problem obtained from (1) by replacing thevalue functions x 7→ Φ(q(ξ i ), h(ξ i ) − T (ξ i )x) by approximations Θi : Rm → R:
¯
(
)
¯
N
¯
min hc, xi + ∑ piΘi (x) ¯ x ∈ X ,
(25)
¯
i=1
where each function Θi (·), i = 1, . . . , N, is given in the form
Θi (x) := max{min{η1 (x), . . . , ηk (x)} | (η1 (·), . . . , ηk (·)) ∈ Ci },x ∈ X̄,Progress in two-stage mixed-integer stochastic programming

17

and a tuple η := (η1 (·), . . . , ηk (·)) ∈ Ci consists of k (where k is allowed to vary withη) affine linear functions η j (·) : Rm → R, j = 1, . . . , k. The tuple η takes here therole of an optimality cut in Benders decomposition for the continuous case. That is,each η ∈ Ci is constructed in a way such that for all x ∈ X
Φ(q(ξ i ), h(ξ i ) − T (ξ i )x) ≥ η j (x) for at least one j ∈ {1, . . . , k}.

(26)

Hence, we have Φ(q(ξ i ), h(ξ i ) − T (ξ i )x) ≥ Θi (x) and the optimal value of problem(25) is a lower bound to the optimal value of (1). Before discussing the constructionof the tuples η, we shortly discuss an algorithm to solve problem (25).

5.2.1 Solving the master problem
Note that problem (25) can be written as a disjunctive mixed-integer linear problem:
¯
(
)
¯
N
¯ x ∈ X,min hc, xi + ∑ pi θi ¯
. (27)¯ θi ≥ η1 (x) ∨ . . . ∨ θi ≥ ηk (x), η ∈ Ci , i = 1, . . . , N
i=1
Problem (27) can be solved by a branch-and-bound algorithm [NS08b]. To thisend, assume that for each tuple η ∈ Ci an affine linear function η̄(·) : Rm → R
is known which underestimates each η j (·), j = 1, . . . , k, on X, i.e., η j (x) ≥ η̄(x) for
j = 1, . . . , k and x ∈ X. η̄(·) allows to derive a linear relaxation of problem (27):
¯
(
)
¯
N
¯
min hc, xi + ∑ pi θi ¯ x ∈ X̄, θi ≥ η̄(x), η ∈ Ci , i = 1, . . . , N .
(28)
¯
i=1
Let (x̂, θ̂ ) be a solution of (28). If x̂ is feasible for (27), then an optimal solutionfor (27) has been found. Otherwise, x̂ either violates an integrality restriction on avariable x j , j = 1, . . . , m0 , or a disjunction θi ≥ min{η1 (x), . . . , ηk (x)} for some tuple η ∈ Ci (with k > 1) and some scenario i. In the former case, that is, x̂ j 6∈ Z,
two subproblems of (28) are created with additional constraints x j ≤ bx̂ j c andx j ≥ dx̂ j e, respectively. In the latter case, the tuple η is partitioned into two tuples
η 0 = (η1 (·), . . . , ηk0 (·)) and η 00 = (ηk0 +1 (·), . . . , ηk (·)), 1 ≤ k0 < k, corresponding
linear underestimators η̄ 0 (·) and η̄ 00 (·) are computed (where η̄ 0 = η1 if k0 = 1 andη̄ 00 = ηk if k0 = k − 1), and two subproblems where the tuple η ∈ Ci is replaced byη 0 and η 00 , respectively, are constructed. Next, the same method is applied to each
subproblem recursively. The first feasible solution for problem (27) is stored as “incumbent solution”. In the following, new feasible solutions replace the incumbent
solution if they have a better objective value. If a subproblem is infeasible or thevalue of its linear relaxation exceeds the current incumbent solution, then it can bediscarded from the list of open subproblems. Since in each subproblem the number of feasible discrete values for a variable x j or the length of a tuple η ∈ Ci is
reduced with respect to the ascending problem, the algorithm can generate only afinite number of subproblems and thus terminates with a solution of (27).

18

Werner Römisch and Stefan Vigerske

5.2.2 Convexification of disjunctive cuts
A linear function η̄(·) in (28) that underestimates min{η1 (·), . . . , ηk (·)} can be constructed by means of disjunctive programming [Ba98, SeH05]: For a fixed scenario
index i and a tuple ηS∈ Ci , an inequality θ ≥ η̄(x) is valid for the feasible set of(27), if it is valid for kj=1 {(x, θ ) ∈ Rm+1 | x ∈ X̄, θ ≥ η j (x)}. That is, we requireη̄(x) ≤ η j (x)for all x ∈ X̄,j = 1, . . . , k.

(29)

We write η̄(x) = η̄0 + hη̄x , xi and η j (x) = η j,0 + hη j,x , xi for some η̄0 , η j,0 ∈ R andη̄x , η j,x ∈ Rm , j = 1, . . . , k. Then (29) is equivalent toη̄0 − η j,0 ≤ min{hη j,x − η̄x , xi | x ∈ Rm , Ax ≤ b}r

= max{hλ j , bi | λ j ∈ R−0 , A> λ j = η j,x − η̄x }.
r

Therefore, choosing λ j ∈ R−0 and η̄x ∈ Rm such that A> λ j + η̄x = η j,x , and settingη̄0 := η j,0 + min{hλ j , bi| j = 1, . . . , k} yields a function η̄(x) that satisfies (29).
[SeH05] note, that given an extreme point x̂ of X̄, the linear underestimator η̄(·)
can be chosen such that η̄(x̂) = min{η1 (x̂), . . . , ηk (x̂)}. Thus, if only extreme pointsof X̄ are feasible for (1), then it is not necessary to branch on disjunctions η to solve(27). This is the case, e.g., if all first stage variables are restricted to be binary.

5.2.3 Approximation of Φ(u,t) by linear optimality cuts
The simplest way to construct a tuple η with property (26) is to derive a supporting
hyperplane for the linear relaxation of Φ(u,t), which we denote byΦ̄(u,t) := min{hu1 , y1 i + hu2 , y2 i | y1 ∈ Rm1 , y2 ∈ Rm2 , W1 y1 +W2 y2 ≤ t}. (30)
It is well known, that Φ̄(u,t) is piecewise linear and convex in t. Thus, if, for fixed(û, tˆ) ∈ U × T , π̂ is a dual solution of (30), we obtain the inequality Φ̄(û,t) ≥Φ̄(û, tˆ) + hπ̂,t − tˆi = hπ̂,ti (t ∈ T ). Letting û = q(ξ i ) and tˆ = h(ξ i ) − T (ξ i )x̂ for a
fixed scenario ξ i and first stage decision x̂ ∈ X̄, we obtainΦ̄(q(ξ i ), h(ξ i ) − T (ξ i )x) ≥ hπ̂, h(ξ i ) − T (ξ i )xi =: η1 (x).

(31)

Since Φ(u,t) ≥ Φ̄(u,t), (31) yields the optimality cut η := (η1 (·)) (i.e., k = 1).
Due to the polyhedrality of Φ̄(u,t), a finite number of such cuts for each scenario issufficient to obtain an exact representation of Φ̄(u,t) in the master problem (27).

5.2.4 Approximation of Φ(u,t) by lift-and-project
However, in order to capture the nonconvexity of the original value function Φ(u,t),nonconvex optimality cuts are necessary, i.e., tuples η of length k > 1. For the caseProgress in two-stage mixed-integer stochastic programming

19

that the discrete variables in the second stage are all of binary type, the following method is proposed in [SeH05]: Let x̄ ∈ X be a feasible solution of the masterproblem (25), let (ȳi1 , ȳi2 ) be a solution of the relaxed second stage problem (30) foru = q(ξ i ) and t = h(ξ i ) − T (ξ i )x̂, i = 1, . . . , N. If ȳi2 ∈ Zm2 for all i = 1, . . . , N, thena linear optimality cut (31) is derived, c.f. (31). Otherwise, let j ∈ {1, . . . , m2 } be anindex such that 0 < ȳi2, j < 1. We now seek for inequalities hπ1i , y1 i+hπ2i , y2 i ≥ π0i (x),i = 1, . . . , N, which are valid for (3) for all x ∈ X̄, but cut off the solution ȳi from(30) for at least one scenario i with fractional ȳi2, j . That is, we search for inequalities
that are valid for the disjunctive sets
¯
¯
½
¾ ½
¾
¯ W y +W2 y2 ≤ t,
¯ W y +W2 y2 ≤ t,y ∈ Rm1 +m2 ¯¯ 1 1∪ y ∈ Rm1 +m2 ¯¯ 1 1, (32)
y2, j ≤ 0−y2, j ≤ −1
where t = h(ξ i ) − T (ξ i )x̄, i = 1, . . . , N. Observe, that points with fractional y2, j are
not contained in (32). With an argumentation similar to the derivation of η̄(·) before,
it follows that, for fixed x, valid inequalities for (32) are described by the systemi
W1> λ1,1= π1i ,ii
W2> λ1,1+ e j λ1,2= π2i ,ii ≥ π0i (x),hh(ξ i ) − T (ξ i )x, λ1,1ii
λ1,1∈ Rr− , λ1,2
∈ R− ,i
W1> λ2,1= π1i ,ii
W2> λ2,1− e j λ2,2= π2i ,iihh(ξ i ) − T (ξ i )x, λ2,1
− λ2,2i ≥ π0i (x),iiλ2,1∈ Rr− , λ2,2∈ R− ,

(33a)
(33b)
(33c)
(33d)

where e j ∈ Rm2 is the j-th unit vector. Observe further, that the coefficients in (33a)and (33b) (i.e., W1 , W2 , e j ) are scenario independent. Thus, it is possible to use
common cut coefficients (π1 , π2 ) ≡ (π1i , π2i ) for all scenarios, thereby reducing thecomputational effort to the solution of a single linear program [SeH05]:¯¯ λ1,1 , λ2,1 ∈ Rr− , λ1,2 , λ2,2 ∈ R− ,


¯


¯ π1 ∈ Rm1 , π2 ∈ Rm2 , π0i (x̄) ∈ R,¯¯ W > λ = π ,W > λ + e λ = π , 
j 1,2
1
2 

¯ 1 1,1
2 1,1
¯ W > λ = π ,W > λ − e λ = π , 
j 2,2
1
2
2 2,1
i
i
i ¯ 1 2,1
max ∑ pi (π0 (x̄) − hπ1 , ȳ1 i − hπ2 , ȳ2 i) ¯¯ hh(ξ i ) − T (ξ i )x̄, λ1,1 i ≥ π0i (x̄),



i=1
¯



i ) − T (ξ i )x̄, λ i − λi (x̄), 
¯


hh(ξ
≥
π
2,1
2,2


0

¯



¯ kπ k ≤ 1, kπ k ≤ 1, |π i (x̄)| ≤ 1, 



∞
∞
1
2


¯
0


¯ i = 1, . . . , N












N

The objective function of this simple recourse problem maximizes the average violation of the computed cuts by (ȳi1 , ȳi2 ). The functions π0i (·), i = 1, . . . , N, with
hπ1 , y1 i + hπ2 , y2 i ≥ π0i (x) for all x ∈ X̄ is derived from a solution of this LP as
π0i (x) := min{hh(ξ i ) − T (ξ i )x, λ1,1 i, hh(ξ i ) − T (ξ i )x, λ2,1 i − λ2,2 }.

(34)

Adding these new cuts to (30) for u = q(ξ i ) and t = h(ξ i ) − T (ξ i )x̄, i = 1, . . . , N,
yields the updated second stage linear relaxations

20

Werner Römisch and Stefan Vigerske¯¯
W1 y1 +W2 y2 ≤ h(ξ i ) − T (ξ i )x̄ ¯. (35)min hq1 (ξ i ), y1 i + hq2 (ξ i ), y2 i ¯¯ −hπ1 , y1 i − hπ2 , y2 i ≤ −π0i (x̄)
¯ y1 ∈ Rm1 , y2 ∈ Rm2
A dual solution (µ, µ0 ) of (35) can then be used to derive the inequality
Φ(q(ξ i ), h(ξ i ) − T (ξ i )x) ≥ hµ, h(ξ i ) − T (ξ i )xi − µ0 π0i (x).
However, the nonconvexity of the right-hand side π0i (x) yields a nonconvex optimality cut η := (η1 (·), η2 (·)), whereη1 (x) :=hµ − µ0 λ1,1 , h(ξ i ) − T (ξ i )xi,
η2 (x) :=hµ − µ0 λ2,1 , h(ξ i ) − T (ξ i )xi + µ0 λ2,2 .
In a next iteration, when the second stage problems are revisited with a different firststage solution x̄, the updated relaxation (35) takes the place of the original relaxation(30). Since the functions π0i (·) are known, the right-hand side of the added cut in (35)is updated when x̄ changes.

5.2.5 Approximation of Φ(u,t) by branch-and-boundFor the general case where the discrete second stage variables can also be of integer
type, the second stage problem (3) can be solved by a (partial) branch-and-boundalgorithm and a (nonlinear) optimality cut η can be derived from the dual solutionsof the linear programs in each leaf of the branch-and-bound tree [SeS06]: Let x̄ ∈X be again a feasible point to problem (25) and fix a scenario i. Assume that (3)with ū = q(ξ i ) and t¯ = h(ξ i ) − T (ξ i )x̄ is (partially) solved by a branch-and-boundalgorithm. Denote by Q the set of terminal nodes of the generated branch-andbound tree. For any node q ∈ Q let yq2,l and yq2,u denote the vectors that define lower
and upper bounds on the y2 variables in the subproblem at node q. Then the LP
relaxation of (3) for node q ∈ Q is given as¯¾½¯
y2 ≤ yq2,u ,. (36)min hū1 , y1 i + hū2 , y2 i ¯¯ y ∈ Rm1 +m2 ,W1 y1 +W2 y2 ≤ t¯,q−y2 ≤ −y2,lWe assume, that subproblems have been pruned if they are infeasible or their lower
bound exceeds a known upper bound. Thus, all terminal nodes are associated witha feasible LP relaxation. The dual problem to (36) is¯½¾¯µ ∈ Rr− ,W1> µ = ū1qq
max hµ, t¯i + hπu , y2,u i − hπl , y2,l i ¯¯, (37)m2
πl , πu ∈ R− , W2> µ + πu − πl = ū2where we assume that a dual variable πl, j , πu, j is fixed to 0 if the corresponding
bound yq2,l, j , yq2,u, j is −∞ or +∞, respectively, j = 1, . . . , m2 . Based on a dual solution(µ q , πlq , πuq ) of (37), a supporting hyperplane of each nodes LP value function can

Progress in two-stage mixed-integer stochastic programming

21

be derived, c.f. (31). Since the branch-and-bound tree represents a partition of thefeasible set of (3), it allows to state a disjunctive description of the function t 7→Φ(ū,t) by combining the LP value function approximations in all nodes q ∈ Q:
Φ(ū,t) ≥ hµ q ,ti + hπuq , yq2,u i − hπlq , yq2,l i for at least one q ∈ Q.

(38)

This result directly translates into a nonlinear optimality cut η := (η1 (·), . . . , ηq (·))
by letting ηq (x) := hµ q , h(ξ i ) − T (ξ i )xi + hπuq , yq2,u i − hπlq , yq2,l i.
5.2.6 Full Algorithm
We can now state a full algorithm for the solution of (1):
1. Solve the master problem (27) by branch-and-bound. If it is infeasible, then (1)is infeasible. Otherwise, let (x̄, θ̄ ) be a solution of (27).
2. Solve (3) for each scenario i = 1, . . . , N. Let φi := Φ(q(ξ i ), h(ξ i ) − T (ξ i )x̄) be
the optimal value of (3) for the first stage decision x̄ in scenario i.
3. For scenarios i where φi > θ̄i , derive an optimality cut η of the value function x 7→ Φ(q(ξ i ), h(ξ i ) − T (ξ i )x) either via linearization of Φ̄(u,t) (see (31)),
via lift-and-project (Section 5.2.4), or from a (partial) branch-and-bound search(Section 5.2.5). Add η to Ci in (27).
4. If no new tuples η have been constructed, i.e., the master problem has not been
updated, then finish: x̄ is an optimal solution to (3). Otherwise, go back to 1.
Some remarks are in order:• At the beginning, the sets Ci are empty, i.e., no information about the value function Φ(u,t) is available in (27). Thus, (27) should be solved either with the variables θi removed or bounded from below by a known lower bound on Φ(u,t).
• In the first iterations, when almost no information about Φ(u,t) is available, it isunnecessary to solve the master problem (27) and the second stage problems (3)to optimality. Instead, at first it is more efficient to ignore the integrality conditions and to construct a representation of the LP value function Φ̄(u,t) by a usual
Benders decomposition. Later, partial solves of (27) and the introduction of nonlinear optimality cuts η into (27) based on lift-and-project or partial branch-andbound searches should be performed to capture the nonconvexity of Φ(u,t) in themaster problem. Finally, to ensure convergence, first and second stage problems
need to be solved to optimality, see also [SeH05, NS08b].

5.2.7 Extension to multistage problems
While the algorithms discussed so far allow an efficient extension of the Benders
decomposition to two-stage mixed-integer stochastic programs, a further extension
to the multistage case seems possible. While in the two-stage case we have a nonconvex value function only in the first stage, in the multistage setting we are faced

22

Werner Römisch and Stefan Vigerske

with such a function in each node of the scenario tree other than the leaves. Thatis, the master problems in each node before the last stage are of the form (27). Approximation of the value function of such a master problem then requires to takethe nonlinear optimality cuts which approximate the value functions of successor
nodes into account. For that matter, we have seen how such a master problem canbe solved by branch-and-bound (Section 5.2.1, [NS08b]) and how an optimality cut
can be derived from a (partial) branch-and-bound search (Section 5.2.5, [SeS06]).
However, the efficiency of such an approach might suffer under the large numberof disjunctions that are induced from optimality cuts on late stages into the masterproblems on early stages. That is, while in the two-stage case the disjunctions in(27) are caused only by integrality constraints on the second stage, in the multistage
setting we have to deal with disjunctions that are induced by disjunctions on succeeding stages. Therefore, solving a fairly large mixed-integer multistage stochastic
program to optimality with this approach seems questionable.
Nevertheless, an interesting application are multistage problem that can only besolved efficiently by a temporal decomposition, e.g., stochastic programs with recombining scenario trees [KV07]. For the latter, the recombining nature of the scenario tree leads to coinciding value functions, a property that can be explored bya nested Benders decomposition. Therefore, an extension to the mixed-integer caseby application of the ideas discussed in this section seems promising.

5.3 Scenario Decomposition
Consider the following reformulation of (1) where the first stage variable x is replaced by one variable xi for each scenario i = 1, . . . , N and an explicit nonanticipativity constraint is added:Nmin∑ pi (hc, xi i + hq1 (ξ i ), y1 (ξ i )i + hq2 (ξ i ), y2 (ξ i )i)

i=1i

yi1 ∈ Rm1 ,(39a)yi2 ∈ Zm2 ,i = 1, . . . , N,

(39b)

T (ξ )x +W1 yi1 +W2 yi2 ≤ h(ξ i ),1
2
N

i = 1, . . . , N,

(39c)such that x ∈ X,ii

x = x = ... = x .

(39d)

Problem (39) decomposes into scenariowise subproblems by relaxing the coupling
constraint (39d) [CS99]. The violation of the relaxed constraints is then added as apenalty into the objective function. That is, each subproblem has the form¯()¯
xi ∈ X, yi1 ∈ Rm1 , yi2 ∈ Zm2 ,¯i iDi (λ ) := min Li (x , y ; λ ) ¯,(40)¯ T (ξ i )xi +W1 yi1 +W2 yi2 ≤ h(ξ i ),where λ := (λ 1 , . . . , λ N ) ∈ RmN is the Lagrange multiplier andProgress in two-stage mixed-integer stochastic programming23

Li (xi , yi ; λ ) := pi (hc, xi i + hq1 (ξ i ), y1 (ξ i )i + hq2 (ξ i ), y2 (ξ i )i + hλ i , xi − x1 i),
i = 1, . . . , N. For every choice of λ , a lower bound on (39) is obtained by computingN

D(λ ) := ∑ Di (λ ).

(41)

i=1

That is, to compute (41), the deterministic problem (40) is solved for each scenario.
To find the best possible lower bound, one now searches for an optimal solution to
the dual problem
max{D(λ ) | λ ∈ RmN }.
(42)
The function D(λ ) is a piecewise linear concave function for which subgradientscan be computed from a solution of (40). Thus, solution methods for the nonsmoothconvex optimization problem (42) use a bundle of subgradients of D(λ ) to findpromising values of λ [Ki90].
The primal solutions (xi , yi ), i = 1, . . . , N, of (40), associated with a solution of(42), yield in general not a feasible solution to the original problem. To regain therelaxed nonanticipativity constraint, heuristics are employed that, e.g., select for x
an average or a frequently occurring value among the xi and then possibly resolve
each second stage problem to ensure feasibility.
To find an optimal solution to (39), a branch-and-bound algorithm can be employed. Here, nonanticipativity constraints are insured by branching on the firststage variables. Since the additional bound constraints on xi become part of theconstraints in (40), the lower bound (42) improves by a branching operation.
An alternative to solving the dual problem (42) by a bundle method is proposed in[LuS04]: As shown in [CS99], the problem (42) is equivalent to the primal problemNmin∑ pi (hc, xi i + hq1 (ξ i ), y1 (ξ i )i + hq2 (ξ i ), y2 (ξ i )i)

(43a)

i=1

such that¯½¯(xi , yi ) ∈ conv (x, y1 , y2 ) ¯¯

x ∈ X, y1 ∈ Rm1 , y2 ∈ Zm2 ,T (ξ i )x +W1 y1 +W2 y2 ≤ h(ξ i )¾, (43b)i = 1, . . . , N,x = x2 = . . . = xN .
1

(43c)

This problem is solved by a column generation approach, which constructs an inner approximation of the convex hull in (43b). Feasible solutions for the original
problem are obtained by application of branch-and-bound.
For problems where all first stage variables are restricted to be binary, [AEO03]
propose to relax both nonanticipativity and integrality constraints. Thereby, eachscenario is associated with a branch-and-bound tree that enumerates the integer feasible solutions to the scenario’s subproblem (i.e., the feasible set of (40)). Since eachbranch-and-bound fixes first stage variables to be either 0 or 1, a coordinated searchacross all n branch-and-bound trees allows to select feasible solutions from each
subproblem that satisfy the nonanticipativity constraints. If also continuous vari-24

Werner Römisch and Stefan Vigerske

ables are present in the first stage, [EG+07] propose to “cross over” to a Benders
decomposition whenever the coordinated branch-and-bound search yields solutions
which binary first stage variables satisfy the nonanticipativity constraints and second
stage integer variables are fixed.

6 Application to stochastic thermal unit commitment
We consider a power generation system comprising thermal units and contracts fordelivery and purchase, and describe a model for its cost-minimal operation under
uncertainty in electrical load and in prices for fuel and electricity. Contracts for delivery and purchase of electricity are regarded as special thermal units. It is assumed
that the time horizon is discretized into uniform (e.g., hourly) intervals. Let T andI denote the numbers of time periods and thermal units, respectively. For thermal
unit i in period t, uit ∈ {0, 1} is its commitment decision (1 if on, 0 if off), and xit its
production, with
uit xitmin ≤ xit ≤ xitmax uit(i = 1, . . . , I, t = 1, . . . , T ),(44)

where xitmin and xitmax are the minimum and maximum capacities. Additionally, there
are minimum up/down-time requirements: when unit i is switched on (off), it mustremain on (off) for at least τ̄i (τ i , resp.) periods, i.e.,uiτ − ui,τ−1 ≤ uit(τ = t − τ̄i + 1, . . . ,t − 1),(45)ui,τ−1 − uiτ ≤ 1 − uit(τ = t − τ i + 1, . . . ,t − 1).

(46)

for all t = 1, . . . , T and i = 1, . . . , I. Let Ui denote the set of all pairs (xi , ui ) satisfying
the constraints (44), (45), and (46) for all t = 1, . . . , T . The basic system requirement
is to meet the electrical load dt during all time periods t = 1, . . . , T , i.e.,I∑ xit ≥ dt

(t = 1, . . . , T ).

(47)

i=1

The expected total system cost is given by the sum of expected startup and operating
costs of all thermal units over the whole scheduling horizon, i.e.,
!
Ã
T

E

I

∑ ∑ (Cit (xit , uit ) + Sit (ui )) .

(48)

t=1 i=1

The fuel cost Cit for operating thermal unit i is assumed to be piecewise linear convex (concave for purchase contracts) during period t, i.e.,
Cit (xit , uit ) := max { ailt xit + bilt uit }
l=1,...,l¯Progress in two-stage mixed-integer stochastic programming25

with cost coefficients ailt , bilt . The startup cost of unit i depends on its downtime;
it may vary between a maximum cold-start value and a much smaller value when
the unit is still relatively close to its operating temperature. This is modeled by thestartup costÃ!
Sit (ui ) := max c ciτ
τ=0,...,τiτ

uit − ∑ ui,t−κ,κ=1where 0 = ci0 < . . . < ciτic are cost coefficients, τic is the cool-down time of unit i, ciτic
T , and u ∈ {0, 1} for τ = 1 − τ c , . . . , 0
is its maximum cold-start cost, ui := (uit )t=1iτi
are given initial values.
T is given by
It is assumed that the stochastic input process ξ = {ξ }t=1
ξt := (at , bt , ct , dt ) (t = 1, . . . , T ).
or by some of its components. Furthermore, it is assumed that ξ1 , . . . , ξt1 (i.e., theinput data for the first time period for which reliable forecasts are available), and,thus, the (first stage) decisions {(xit , uit ) | t = 1, . . . ,t1 , i = 1, . . . , I} are deterministic.
Minimizing the expected total cost (48) such that the operational constraints(44), (45), (46), and (47) are satisfied, represents a two-stage (linear) mixedinteger stochastic program with (random) second stage decision {(xit , uit ) | t =t1 + 1, . . . , T, i = 1, . . . , I}.
In many cases it is possible to derive a model for the probability distribution Pof ξ via time series analysis based on historical data (see, e.g., [ERW05, SYG06]).
Sampling from P together with applying scenario reduction (see Section 4) thenleads to a finite number of scenarios ξ j = (atj , btj , ctj , dtj ) with probabilities p j , j =1, . . . , N, for the stochastic process ξ and to the corresponding decision scenarios(xij , uij ) (for unit i). The scenario-based unit commitment problem then reads¯ j j¯ (x , u ) ∈ Ui , i = 1, . . . , I, i i¯min ∑ ∑ ∑ p j (Citj (xitj , uitj ) + Sitj (uij )) ¯¯ ∑Ii=1 pitj ≥ dtj , t = 1, . . . , T,
 j=1 t=1 i=1¯ j = 1, . . . , N

N

T

I

(49)

where Citj and Sitj denote the cost functions for scenario j.
Since the optimization problem (49) only contains NT (unit) coupling constraints
while the number 2NT I of decision variables is typically (much) larger, geographical decomposition based on Lagrangian relaxation of the coupling constraints (47)
seems to be promising. The Lagrangian function is of the form
Ã
!
N

T

L(x, u; λ ) = ∑ ∑ p j
j=1 t=1
NT= ∑ ∑ pj
j=1 t=1

I

I

i=1

i=1

∑ (Citj (xitj , uitj ) + Sitj (uij )) + λt j (dtj − ∑ xitj )ÃI(Citj (xitj , uitj ) + Sitj (uij ) − λt j xitj )) + λt j dtj
i=1∑

which leads to the dual function

!
,

26

Werner Römisch and Stefan VigerskeÃND(λ ) = inf L(x, u; λ ) = ∑ p j
(x,u)

I

T

∑ Di j (λ ) + ∑
j

j=1

i=1

!
λt j dtj

t=1

T

Di j (λ ) = inf ∑ (inf(Citj (xit , uit ) − λt xit ) + Sitj (ui ))
ui

t=1

xit

decomposing into unit subproblems for every scenario ξ j and, hence, λ j . While theinner minimization (with respect to xit ) can be solved explicitly, the outer minimization (with respect to ui ) can efficiently be done by dynamic programming. The dual
concave (nondifferentiable) maximization problem©
ªmax D(λ ) | λ ∈ RNT(50)+
can be solved by bundle subgradient methods (e.g., [Ki90]). If (x̄, ū, λ̄ ) is an (approximate) solution of (50), D(λ̄ ) is a lower bound of the infimum of (49), but, ingeneral, the (maximal) load constraints
I∑ xitmax ūitj ≥ dtj(t = 1, . . . , T, j = 1, . . . , N)

(51)

i=1

are violated for some scenarios j and some time intervals t, respectively. However,as shown in [Be82, Section 5.6.1], the relative duality gap gets small if the number Iof units is large. In many practical situations this allows to apply simple Lagrangianheuristics (like [ZG88]) to modify ū scenariowise such that (51) is satisfied for every pair (t, j). After having the commitment decision ū fixed, a final scenariowise
economic dispatch [vBL87] leads to good primal solutions (x̄, ū).
The approach can be extended to multistage models by requiring in addition thatthe decisions (xt , ut ) in (49) only depend on (ξ1 , . . . , ξt ) (for t > t1 ). We refer to therelevant work [CC+96, GK+02, GR05, NR00, PCW00, SYG06, TBL96, TKW00].
Furthermore, instead of the expected total system cost, a mean-risk objective of
the formtIγρ(Yt1 , . . . ,YT ) − (1 − γ)E(YT ), Yt := − ∑ ∑ (Cit (xit , uit ) + Sit (ui )) (t = t1 , . . . , T )τ=1 i=1

may be considered, where γ ∈ (0, 1) and ρ is a multi-period risk functional (see[EHR09]). In this way, risk management is integrated into unit commitment. If therisk functional is polyhedral [ER05, EHR09], the scenario-based unit commitment
model may be reformulated as mixed-integer linear program.
Extensions of the two-stage stochastic unit commitment model are discussed
in [NR02] and [NSW05], respectively. In [NR02], a planning model is described
whose (deterministic) first stage and (stochastic) second stage decisions are givenon the whole time horizon {1, . . . , T }. The first stage decisions are determined suchthat a transition from the first to the second stage and vice versa is always feasible

Progress in two-stage mixed-integer stochastic programming

27

and compatible. In [NSW05], day-ahead trading at a power exchange is incorporated into unit commitment.

7 Conclusions
We reviewed recent progress in two-stage mixed-integer stochastic programming.
First we reviewed structural properties of optimal value functions of mixed-integer
linear programs from the literature and discussed conclusions for continuity properties of integrands in two-stage mixed-integer stochastic programs. If the probability distribution has finite support, the expected recourse function is piecewise
continuous with a finite number of polyhedral continuity regions. When perturbing
or approximating the underlying probability distribution, the optimal value function behaves continuous with respect to a discrepancy distance of the original andperturbed probability measures. This result allowed to extend the stability based
scenario reduction algorithm from [DGR03, HR07] to the mixed-integer two-stage
situation.
For solving a two-stage mixed-integer stochastic program, several decomposition
algorithms are reviewed. First, methods to convexify the expected recourse functionof simple and more complex integer-recourse models by perturbing the probability
measure are discussed. This allows to obtain tight bounds on the original optimalvalue. Secondly, algorithms that decompose the stochastic program in a Benders
decomposition style are detailed. Here, the nonconvexity in the second-stage valuefunctions is captured by nonlinear optimality cuts, which might make a solutionof the master problem by branch-and-bound necessary. Further, scenario decomposition based algorithms based on relaxation of nonanticipativity constraints arereviewed. In a Lagrangian decomposition, the subproblems are coupled via a dual
problem which comprises the maximization of a piecewise linear concave function.
Finally, a geographical Lagrangian decomposition method is illustrated on astochastic thermal unit commitment problem. Here, the problem decomposed intoone (stochastic mixed-integer) subproblem for each thermal unit. This allows to exploit the subproblems structure by specialized algorithms and to use a Lagrangianheuristic specialized for unit commitment problems.
Acknowledgements This work was supported by the DFG Research Center M ATHEON Mathematics for Key Technologies in Berlin (http://www.matheon.de) and by the German Ministry of Education and Research (BMBF) under the grant 03SF0312E.

References[BG+82]

Bank, B., Guddat, J., Kummer, B., Klatte, D., Tammer, K.: Non-linear Parametric
Optimization. Akademie-Verlag, Berlin (1982)

28[AEO03][ATS00][Ba98][Be82][Bi85][BJ82][CS99][CT97][CT98][CC+96][DR04][DGR03][EHR09][ER05][ER07][ERW05][EG+07][GK+02][GR05][HR07][HKR08][HKR09]Werner Römisch and Stefan Vigerske
Alonso-Ayuso, A., Escudero, L.F., Teresa Ortuńo, M.: BFC, a branch-and-fix coordination algorithmic framework for solving some types of stochastic pure and mixed
0-1 programs. Eur. J. Oper. Res. 151, 503–519 (2003)
Ahmed, S., Tawarmalani, M., Sahinidis, N.V.: A finite branch and bound algorithm
for two-stage stochastic integer programs. Math. Program. 100, 355–377 (2004)
Balas, E.: Disjunctive programming: Properties of the convex hull of feasible points.
Discrete Appl. Math. 89, 3–44 (1998). originally MSRR#348, Carnegie Mellon University (1974)
Bertsekas, D.P.: Constrained Optimization and Lagrange Multiplier Methods. Academic Press, New York (1982)
Birge, J.R.: Decomposition and Partitioning Methods for Multistage Stochastic Programming. Oper. Res. 33(5), 989–1007 (1985)
Blair, C.E., Jeroslow, R.G.: The value function of an integer program. Math. Progr.
23, 237–273 (1982)
Carøe, C.C., Schultz, R.: Dual decomposition in stochastic integer programming.
Oper. Res. Lett. 24(1-2), 37–45 (1999)
Carøe, C.C., Tind, J.: A cutting-plane approach to mixed 0 − 1 stochastic integer
programs. Eur. J. Oper. Res. 101(2), 306–316 (1997)
Carøe, C.C., Tind, J.: L-shaped decomposition of two-stage stochastic programs withinteger recourse. Math. Progr. 83(3, Ser. A), 451–464 (1998)
Carpentier, P., Cohen, G., Culioli, J.-C., Renaud, A.: Stochastic optimization of unit
commitment: A new decomposition framework. IEEE Trans. Power Systems 11,1067–1073 (1996)
Dentcheva, D., Römisch, W.: Duality gaps in nonconvex stochastic optimization.
Math. Progr. 101(3, Ser. A), 515–535 (2004)
Dupačová, J., Gröwe-Kuska, N., Römisch, W.: Scenarios reduction in stochastic programming: An approach using probability metrics. Math. Progr. 95, 493–511 (2003)
Eichhorn, A., Heitsch, H., Römisch, W.: Stochastic optimization of electricity portfolios: Scenario tree modeling and risk management. In: Kallrath, J., Pardalos, P.
(eds.) Power Systems Handbook. Springer (to appear)
Eichhorn, A., Römisch, W.: Polyhedral risk measures in stochastic programming.
SIAM J. Optim. 16, 69–95 (2005)
Eichhorn, A., Römisch, W.: Stochastic integer programming: Limit theorems andconfidence intervals. Math. Oper. Res. 32, 118–135 (2007)
Eichhorn, A., Römisch, W., Wegner, I.: Mean-risk optimization of electricity portfolios using multiperiod polyhedral risk measures. In: IEEE St. Petersburg PowerTech
Proceedings (2005)
Escudero, L.F., Garı́n, A., Merino, M., Pérez, G.: A two-stage stochastic integer programming approach as a mixture of branch-and-fix coordination and Benders decomposition schemes. Ann. Oper. Res. 152, 395–420 (2007)
Gröwe-Kuska, N., Kiwiel, K.C., Nowak, M.P., Römisch, W., Wegner, I.: Power management in a hydro-thermal system under uncertainty by Lagrangian relaxation. In:
Greengard, C., Ruszczyński, A. (eds.) Decision Making under Uncertainty: Energy
and Power, pp. 39–70. Springer, New York (2002).
Gröwe-Kuska, N., Römisch, W.: Stochastic unit commitment in hydro-thermal
power production planning. In Wallace, S.W., Ziemba, W.T. (eds.) Applications of
Stochastic Programming, pp. 633–653. MPS/SIAM Series on Optimization, SIAM,
Philadelphia (2005)
Heitsch, H., Römisch, W.: A note on scenario reduction for two-stage stochastic
programs. Oper. Res. Lett. 35, 731–738 (2007)
Henrion, R., Küchler, C., Römisch, W.: Discrepancy distances and scenario reduction
in two-stage stochastic mixed-integer programming. J. Ind. Manag. Optim. 4, 363–
384 (2008)
Henrion, R., Küchler, C., Römisch, W.: Scenario reduction in stochastic programming with respect to discrepancy distances, Comput. Optim. Appl. 43, 67–93 (2009)

Progress in two-stage mixed-integer stochastic programming[HKR09b][Ki90][KSV06][KV07][LaL93][LoS03][LuS04][Ni92][NR00][NSW05][NS05][NS08a][NS08b][NR02][PCW00][Ra91][RS01][RV08][Sch96][Sch03][SSV98][ST06][Se05]

29

Henrion, R., Küchler, C., Römisch, W.: A scenario reduction heuristic for two-stage
stochastic integer programs. in preparation
Kiwiel, K.C.: Proximity control in bundle methods for convex nondifferentiable optimization. Math. Progr. 46, 105–122 (1990)
Klein Haneveld, W.K., Stougie, L., van der Vlerk, M.H.: Simple integer recourse
models: convexity and convex approximations. Math. Progr. 108(2-3, Ser. B), 435–
473 (2006)
Küchler, C., Vigerske, S.: Decomposition of multistage stochastic programs with recombining scenario trees. Stoch. Progr. E-Print Series 9 (2007). www.speps.org
Laporte, G., Louveaux, F.V.: The integer L-shaped method for stochastic integer programs with complete recourse. Oper. Res. Lett. 13(3), 133–142 (1993)
Louveaux, F.V., Schultz, R.: Stochastic Integer Programming. In: Ruszczyński, A.,
Shapiro, A. (eds.) Stochastic Programming, pp. 213–266. Handbooks in Operations
Research and Management Science Vol. 10. Elsevier (2003)
Lulli, G., Sen, S.: A branch and price algorithm for multi-stage stochastic integer
programs with applications to stochastic lot sizing problems. Manage. Sci. 50, 786–
796 (2004)
Niederreiter, H.: Random Number Generation and Quasi-Monte Carlo Methods.
SIAM, Philadelphia (1992)
Nowak, M.P., Römisch, W.: Stochastic Lagrangian relaxation applied to powerscheduling in a hydro-thermal system under uncertainty. Ann. Oper. Res. 100, 251–272 (2000)
Nowak, M.P., Schultz, R., Westphalen, M.: A stochastic integer programming model
for incorporating day-ahead trading of electricity into hydro-thermal unit commitment. Optim. Eng. 6, 163–176 (2005)
Ntaimo, L., Sen, S.: The million-variable “march” for stochastic combinatorial optimization. J. Glob. Optim. 32(3), 385–400 (2005)
Ntaimo, L., Sen, S.: A comparative study of decomposition algorithms for stochastic
combinatorial optimization. Comput. Optim. Appl. 40(3), 299–319 (2008)
Ntaimo, L., Sen, S.: A branch-and-cut algorithm for two-stage stochastic mixedbinary programs with continuous first-stage variables. Int. J. Comp. Sci. Eng. 3, 232–
241 (2008)Nürnberg, R., Römisch, W.: A two-stage planning model for power scheduling in ahydro-thermal system under uncertainty. Optim. Eng. 3, 355–378 (2002)
Philpott, A.B., Craddock, M., Waterer, H.: Hydro-electric unit commitment subject
to uncertain demand. Eur. J. Oper. Res. 125, 410–424 (2000)
Rachev, S.T.: Probability Metrics and the Stability of Stochastic Models. Wiley,
Chichester (1991)
Römisch, W., Schultz, R.: Multistage stochastic integer programs: An introduction.
In: Grötschel, M., Krumke, S.O., J. Rambau, J. (eds.) Online Optimization of Large
Scale Systems, pp. 581–600. Springer, Berlin (2001)
Römisch, W., Vigerske, S.: Quantitative stability of fully random mixed-integer twostage stochastic programs. Optim. Lett. 2, 377–388 (2008)
Schultz, R.: Rates of convergence in stochastic programs with complete integer recourse. SIAM J. Optim. 6, 1138–1152 (1996)
Schultz, R.: Stochastic programming with integer variables. Math. Progr. 97, 285–
309 (2003)
Schultz, R., Stougie, L., van der Vlerk, M.H.: Solving stochastic programs withinteger recourse by enumeration: a framework using Gröbner basis reductions.
Math. Progr. 83(2, Ser. A), 229–252 (1998)
Schultz, R., Tiedemann, S.: Conditional value-at-risk in stochastic programs with
mixed-integer recourse. Math. Progr. 105, 365–386 (2006)
Sen, S.: Algorithms for stochastic mixed-integer programming models. In: Aardal,K., Nemhauser, G.L., Weismantel, R. (eds.) Handbook of Discrete Optimization, pp.
515–558. North-Holland Publishing Co. (2005)

30[SYG06][SeH05][SeS06][TBL96][TKW00][TiW81][vBL87][vdV04][vdV07][vdV05][VSW69][Wo81][ZG88]

Werner Römisch and Stefan VigerskeSen, S., Yu, L., Genc, T.: A stochastic programming approach to power portfolio
optimization. Oper. Res. 54, 55–72 (2006)
Sen, S., Higle, J.L.: The C3 theorem and a D2 algorithm for large scale stochastic
mixed-integer programming: set convexification. Math. Progr. 104(A), 1–20 (2005)
Sen, S., Sherali, H.D.: Decomposition with branch-and-cut approaches for two-stage
stochastic mixed-integer programming. Math. Progr. 106(2, Ser. A), 203–223 (2006)
Takriti, S., Birge, J.R., Long, E.: A stochastic model for the unit commitment problem. IEEE Trans. Power Systems 11, 1497–1508 (1996)
Takriti, S., Krasenbrink, B., Wu, L.S.-Y.: Incorporating fuel constraints and electricity spot prices into the stochastic unit commitment problem. Oper. Res. 48, 268–280(2000)
Tind, J., Wolsey, L.A.: An elementary survey of general duality theory in mathematical programming. Math. Progr. 21, 241–261 (1981)
van den Bosch, P.P.J., Lootsma, F.A.: Scheduling of power generation via large-scalenonlinear optimization. J. Optim. Theory Appl. 55, 313–326 (1987)
van der Vlerk, M.H.: Convex approximations for complete integer recourse models.
Math. Progr. 99(2, Ser. A), 297–310 (2004)
van der Vlerk, M.H.: Stochastic integer programming bibliography. http://mally.eco.rug.nl/spbib.html (1996-2007)
van der Vlerk, M.H.: Convex approximations for a class of mixed-integer recourse
models. Ann. Oper. Res. (to appear) and Stoch. Progr. E-Print Series (2005)
Van Slyke, R.M., Wets, R.: L-Shaped Linear Programs with Applications to OptimalControl and Stochastic Programming. SIAM J. Appl. Math. 17(4), 638–663 (1969)
Wolsey, L.A.: Integer programming duality: Price functions and sensitivity analysis.
Math. Progr. 20, 173–195 (1981)
Zhuang, G., Galiana, F.D.: Towards a more rigorous and practical unit commitment
by Lagrangian relaxation. IEEE Trans. Power Systems 3, 763–773 (1988)