Birge J.R., Louveaux F. Introduction to Stochastic Programming

Подождите немного. Документ загружается.

234 5 Two-Stage Recourse Problems

Schur complement approach in (5.45) is quick but can lead to numerical instabilities

(see Carpenter, Lustig, and Mulvey [1991]).

Carpenter et al. also propose the column splitting technique. The basic idea is to

rewrite problem (1.1) with explicit constraints on nonanticipativity. The formulation

then becomes:

min

∑

k=1

+ q

) (5.46)

s. t. Ax

=b , (5.47)

+Wy

, k = 1,...,K , (5.48)

−x

k+1

=0 , k = 1,...,K −1 , (5.49)

≥ 0 , y

≥0 , k = 1,...,K . (5.50)

The difference now is that the constraints in (5.47)and(5.48) separate into separate

subproblems k and constraints (5.49) link the problems together. Alternating con-

straints, (5.47), (5.48)and(5.49) for each k in sequence, the full constraint matrix

has the form:

A =

⎛

⎜

⎝

A 0000000

W 000000

I 0 −I 00000

00A 00000

00T

W 0000

00 I 0 −I 000

000 0 I 0 −I 0

000000A 0

000000T

⎞

⎟

⎠

. (5.51)

If we form

, then we obtain

⎛

⎜

⎝

A 00000

00000

+WW

2I −A

000 0

00−AAA

A 00

00T

00 0

+WW

000T

2I 00 0

000A

K−1

2I −A

00000−AAA

000000T

+WW

⎞

⎟

⎠

, (5.52)

5.5 Basis Factorization and Interior Point Methods 235

which is clearly much sparser than the original matrix in (5.18). It is, however, larger

than the matrix in (5.18)(seeExercise8) so there is some tradeoff for the reduced

density.

The third additional approach is to form the dual of (1.1) and to solve that prob-

lem using the same basic interior point method we gave earlier. (In the self-dual

form for (5.14), this corresponds to eliminating the primal variables ﬁrst and then

solving for the dual variables. The primal projective scaling method corresponds

to eliminating the dual variables and then solving for the primal variables. Another

alternative for the self-dual form is directly to solve the full system again taking ad-

vantage of the stochastic program constraint structure.) The dual approach considers

the problem:

max b

∑

k=1

(5.53)

s. t. A

∑

k=1

≤ c , k = 1,...,K , (5.54)

≤ q , k = 1,...,K , (5.55)

where the variables are not restricted in sign. For this problem, we can achieve a

standard form as in (5.12) by splitting the variables

and

into differences

of non-negative variables and by adding slack variables to constraints (5.54)and

(5.55)

. In this way the constraint matrix for (5.54)and(5.55) becomes A



⎛

⎜

⎝

−A

−T

0 T

−T

... T

−T

0 I

00W

−W

I 0000000

00 0 00W

−W

I 0000

00 0 000 00W

−W

I 0

⎞

⎟

⎠

. (5.56)

The matrix in (5.56) may again be much larger than the matrix in the original, but

the gain comes in considering A



T

which is now:

⎛

⎜

⎝

2(A

A +

∑

k=1

)+I 2T

W 2T

W ... 2T

W + I 000

02W

W + I 00

.00

0002W

⎞

⎟

⎠

, (5.57)

with an inherent sparsity of which an interior point method can take advantage. In

fact, it is not necessary to take the dual to use this alternative factorization form

by again using the Sherman-Morrison-Woodbury formula (see Birge, Freund, and

The dual problem may no longer have a bounded set of optima causing some theoretical difﬁcul-

ties for convergence results. In practice, bounds are placed on the variables to guarantee conver-

gence.

236 5 Two-Stage Recourse Problems

Vanderbei [1992]). In this way, the matrix in the form of (5.57) replaces the dense

matrix in (5.18).

In addition to reducing total computational effort, the factorizations described in

this section also allow signiﬁcant parallel processing for the computations involving

sub-matrices corresponding to the second-period subproblems (see, e.g., Yang and

Zenios [1997] and Gondzio and Grothey [2009]). The interior point method fac-

torization also can be extended to multistage problems using a recursive form (see

Pﬂug and Halada [2003]).

An additional strategy for interior point methods is to combine the outer lin-

earization and with an interior point method so that the interior point iterations are

taken with an increasingly constrained region as in the standard L-shaped method

to solve (1.2)–(1.4). This method may reduce the effort in solving subproblems to

optimality while still obtaining reﬁned information about the recourse function con-

straints without requiring full information as in the extensive form. Bahn, Gofﬁn,

du Merle, and Vial [1995] provide a description and computational results for this

approach.

Exercises

1. Use the matrix structure in Proposition 5 to complete the simplex iterations

starting from basis B

for Example 2 from Section 5.1.

2. Compare the number of operations to solve (5.5)using(5.6)and(5.7) compared

to solving (5.5) as an unstructured linear system of equations.

3. Give a similar basis factorization scheme to (5.6)and(5.7) to solve the backward

transformation, (

)B =(c

) , for a basis corresponding to columns B

from the constraint matrix of (1.1).

4. Describe an efﬁcient updating procedure for any possible combination of en-

tering and leaving columns in the basis matrix of (5.5) using the factorization

scheme in (5.6)and(5.7).

5. Find the number of arithmetic operations for a single step of the interior point

method using (5.22). Compare this to the number of arithmetic operations if no

special factorization is used.

6. Prove Proposition 6.

7. Assuming an initial lower bound

= 0 , follow the projective scaling al-

gorithm for Example 2 starting from the x

ext

=(3,7,1,3,1,2,2,1)

until

ext

−

= 0.01 . (Note: the number of iterations required here in com-

parison to the L -shaped method may surprise you, but the number of iterations

for interior point methods generally grows slowly as the problem size increases.)

8. Compare the sizes of the adjacency matrices in (5.18)and(5.51). Assuming

that each matrix A , T

,andW is completely dense, compare the number of

nonzero entries in these two matrices.

5.6 Inner Linearization Methods and Special Structures 237

5.6 Inner Linearization Methods and Special Structures

As mentioned earlier, the most direct alternative to an outer linearization, or cut

generation, approach is an inner linearization or column generation approach (see

Geoffrion [1970] for other basic approaches to large-scale problems). In fact, this

was the ﬁrst suggestion of Dantzig and Madansky [1961] for solving stochastic

linear programs. They observed that the structure of the dual in Figure 2 ﬁts the

prototype for Dantzig-Wolfe decomposition. In fact, we can derive this approach

from the L -shaped method by taking duals.

Consider the following dual linear program to (1.2)–(1.4).

max

b +

∑

=1



∑

=1



(6.1)

s. t.

A +

∑

=1



∑

=1



≤ c

, (6.2)

∑

=1



= 1 ,



≥ 0 ,= 1,...,r ,



≥ 0 ,= 1,...,s . (6.3)

The linear program in (6.1)–(6.3) includes multipliers



on extreme rays, or direc-

tions of recession, which cannot be produced with positive combinations of other

distinct recession directions, of the duals of the subproblems and multipliers



the expectations of extreme points of the duals of the subproblems. To see this, sup-

pose that (6.1)–(6.3) is solved to obtain a multiplier x

on constraint (6.2). Now,

consider the following dual to (1.9):

maxw =

−T

) s.t.

W ≤ q

. (6.4)

If (6.4) is unbounded for any k , we then must have some

such that

W ≤ 0

and

−T

) > 0,or(1.5)–(1.6) has a feasible dual solution (hence optimal

primal solution) with a positive value. So, Step 2 of the L -shaped method is equiv-

alent to checking whether (6.4) is unbounded for any k . In this case, we form D

r+1

and d

r+1

as in (1.7)and(1.8)ofthe L -shaped method and add them to (6.1)–(6.3).

Next, note, that if (6.4) is infeasible, the stochastic program is not well-formulated

(see Exercise 1). Consider when (6.4) has a ﬁnite optimal value for all k .Inthe L -

shaped method, if (1.9) was solvable for all k , then we formed E

s+1

and e

s+1

and

added them to (1.2)–(1.4). In this case in the inner linearization procedure, we again

use (1.10)and(1.11)toform E

s+1

and e

s+1

and add them to (6.1)–(6.3).

Solving the duals in Steps 1 to 3 of the L -shaped algorithm then consists of solv-

ing (6.1)–(6.3) as a master problem and problems (6.4) as subproblems. Formally,

this method is the following inner linearization method.

Inner Linearization Algorithm

Step 0.Set r = s =

= 0.

238 5 Two-Stage Recourse Problems

Step 1.Set

+ 1 and solve the linear program in (6.1)–(6.3). Let the solution

be (

) with a dual solution, (x

) .

Step 2.For k = 1,...,K ,solve(6.4). If any infeasible problem (6.4) is found, stop

and evaluate the formulation. If an unbounded solution with extreme ray

found for any k , then form new columns ( d

r+1

)

, D

r+1

)

), set

r = r + 1 , and return to Step 1.

If all problems (6.4) are solvable, then form new columns, E

s+1

and e

s+1

,asin

(1.10)and(1.11). If e

s+1

−E

s+1

−

≤0 , then stop; (

) and (x

)

are optimal in the original problem (1.2).

If e

s+1

−E

s+1

−

> 0,set s = s+ 1,andreturntoStep1.

Clearly, the inner linearization method follows the same steps as the L -shaped

method, except that we solve the duals of the problems instead of the primals.

Hence, convergence follows directly from the L -shaped method. We could also

view this approach directly as in Dantzig-Wolfe decomposition by stating that (6.1)–

(6.3) is an inner linearization of the dual of the basic L -shaped problem in (1.2)and

that the subproblems (6.4) generate new extreme points and rays to add to this inner

linearization (see Exercise 2).

If, as in many problems, n

>> m

, the primal version has smaller basis ma-

trices, at most of order m

+ m

, than the n

×n

bases for the dual. Hence, the

L -shaped implementation is usually preferred. Inner linearization can, however, be

applied directly to the primal by assuming T is ﬁxed using the form in (3.1.5),

which we repeat here:

min z = c

x +

(

) (6.5)

s. t. Ax = b ,

Tx−

= 0 ,

x ≥ 0 ,

where

(

)=E

(

)) and

(

)) = min{q(

)

y | Wy = h(

) −

,y ≥ 0}. Note that, in this form, we assume that T is ﬁxed but q and h may

still be functions of

. For this reason, we revert to the use of

for the recourse

function.

In this case, we wish to build an inner linearization of the function

(

) us-

ing the generalized programming approach as in Dantzig [1963, Chapter 24]. The

basic idea is to replace

(

) by the convex hull of points

(



) chosen in each

iteration of the algorithm. Each iteration generates a new extreme point of a region

of linearity for

, which is polyhedral as we showed in Theorem 3.6. Thus, ﬁnite

convergence is assured with ﬁnite numbers of realizations. The algorithm follows.

Generalized Programming Method for Two-Stage Stochastic Linear Programs

Step 0.Set s = t =

= 0.

5.6 Inner Linearization Methods and Special Structures 239

Step 1.Set

+ 1 and solve the linear program master problem:

min z

= c

x +

∑

i=1

(

∑

i=1

(

) (6.6)

s. t. Ax = b , (6.7)

Tx−

∑

i=1

−

∑

i=1

= 0 , (6.8)

∑

i=1

= 1 , (6.9)

≥ 0 , i = 1,...,r ,

≥ 0 , i = 1,...,s .

If (6.6)–(6.9) is infeasible or unbounded, stop. Otherwise, let the solution be

) with associated dual variables, (

) .

Step 2. Solve the subproblem:

min

(

)+(

)

−

, (6.10)

which we assume has value less than ∞ .

If (6.10) is unbounded, a recession direction

r+1

is obtained, such that for

some

(

αζ

r+1

)+(

)

(

αζ

r+1

) →−∞ as

→ ∞ . In this case, let

(

r+1

)=lim

→∞

(

αζ

r+1

)−

(

)

, r = r + 1,andreturntoStep1.

If (6.10) is solvable, let the solution be

s+1

.If

(

)+(

)

−

≥ 0,

then stop; (x

) corresponds to an optimal solution to (6.5). Otherwise, set

s = s+ 1 and return to Step 1.

This algorithm generates columns in (6.6)–(6.9) corresponding to new proposals

from the subproblem in (6.10). In the two-stage stochastic linear program form,

(6.10) can be recast as:

min

∑

k=1

)

−

(6.11)

s. t. Wy

= h

, k = 1,...,K ,

≥ 0 , k = 1,...,K .

This problem is not generally separable into different subproblems for each k .

Hence, for general problems, the L -shaped method has an advantage. In some

cases (notably simple recourse),

(

) is separable into components for each k ,

and (6.11) can again be divided into K independent subproblems. We discuss this

possibility further in Section 5.7.

To see how the generalized programming form of inner linearization can be ap-

plied to a stochastic program, we again consider Example 2 from Section 5.1.

240 5 Two-Stage Recourse Problems

Iteration 1:

Suppose we start with an initial solution of

= 1and

(

in (6.6),

which then takes the form:

min z

= 0x

(

) (6.12)

s. t. x

+ x

= 10 , (6.13)

− 0 ·

= 0 , (6.14)

= 1 , (6.15)

≥ 0 ,

which has an optimal solution (x

)=(0,10,1) with dual multipliers

(

)=(0, 0,

) .

Next, for Step 2, the solution is to ﬁnd the minimum value of (6.10)or

(

)−

over

, which is achieved at

= 2 with

(

)=1.Since

(

) −

= 1 −

= −

< 0 , the algorithm returns to Step 1 with

= 2.

Iteration 2:

The solution of (6.6)nowis (x

)=(1,9,0,1) with dual multipliers

(

)=(0, 0,1) . In Step 2, the minimum value for (6.10) occurs at

= 1 and the objective value

(

) −

= 0 , the termination condition.

The steps of this inner linearization algorithm can be viewed as taking the con-

vex hull of an increasing numbers of extreme points of the epigraph of the recourse

function. This can be seen for Example 2 in Figure 3. The solution starts at the point

on the function corresponding to x(=

)=0 and then moves directly to including

the point on the epigraph at x =

= 1 , where no further descent is possible. The

algorithm terminates virtually immediately for this example because the best candi-

date

directly yields an overall optimal solution. This circumstance of course does

not always occur, but the algorithm may be quite efﬁcient when T is ﬁxed and the

subproblems (6.10) can be solved efﬁciently.

To show that the generalized programming method also converges ﬁnitely, we

wish to demonstrate the property of generating extreme points on the epigraph of

byshowingthatanextremesolutionin(6.11) is an extreme value of linear regions

(

) . We do this for extreme points in the following proposition.

Proposition 18. Every optimal extreme point, ( ¯y

,..., ¯y

) , of the feasible region

in (6.11) corresponds to an extreme point

of {

(

},where

∑

k=1

, and each

is an extreme point of {

W ≤ q

Proof: Suppose ( ¯y

,..., ¯y

) is an optimal extreme point in (6.11). In this case,

we must have q

¯y

≤ q

for all Wy

−

. We must also have that ¯y

an extreme point of {y

|Wy

−

≥ 0} because, otherwise, we could take

¯y

=(1/2)(y

+ y

) for distinct feasible y

and y

.So, ¯y

has a complementary

dual solution,

, that is an extreme point of {

W ≤ q

} andsuchthat

−

W) ¯y

= 0.

5.6 Inner Linearization Methods and Special Structures 241

Now, suppose

is not an extreme point of the linearity region where

(

for

(

)−

with

∑

k=1

. In this case, there exists

and

such that

λχ

+(1−

)

where 0 <

< 1,for

(

and

(

.Wealsohavethat

(

∑

k=1

,where q

−

) for j = 1,2 , because, by

feasible in the k -th recourse problem, the only

other possibility is q

(

−

) , which would imply

(

) >

This also implies that

(

W −q

)(

+(1 −

)=0 , (6.16)

which implies that

+(1−

= ¯y

because ¯y

is an extreme point of the fea-

sible region in recourse problem k . In this case, ( ¯y

,..., ¯y

(1 −

)(y

,..., ¯y

) , with both terms feasible in (6.11). This contradicts that

( ¯y

,..., ¯y

) is an extreme point.

A similar argument shows that any extreme ray found in solving (6.11)isan

extreme ray of a region of linearity of

(

) (Exercise 3). Now, we can state the

generalized programming ﬁnite convergence result.

Theorem 19. The generalized programming applied to problem (6.5) with subprob-

lem (6.11) solves (6.5) in a ﬁnite number of steps.

Proof: At each solution of (6.11), a new linear region extreme value is generated.

First for a new extreme ray, we must have

(

r+1

)+(

)

(

r+1

) < 0 , while,

for 1 ≤i ≤s ,

(

) ≥−(

)

. For an extreme point, we only add that point if

(

s+1

)+(

)

s+1

−

< 0 , while, for 1 ≤i ≤s ,

(

)+(

)

−

≥0.

Because the number of such regions is ﬁnite and each has a ﬁnite number of extreme

points and rays, the algorithm converges ﬁnitely.

The solution found solves (6.5) because if we reach the termination condition,

then

(

)

b +

≤ (

)

b +

(

)+(

)

≤ (

A +(

)

T)x +

(

) , (x,

) feasible in (6.5) ,

≤ c

x +

(

) , (6.17)

for all (x,

) feasible in (6.5).

As with the L -shaped method, we can also modify the generalized linear pro-

gramming approach to consider only active columns so that s and t can be bounded

again by m

. Of course, this approach’s greatest potential is in simple recourse

problems as we mentioned earlier. It may also be advantageous if an algorithm can

take advantage of the special matrix structure in (6.11). The most direct approach

in this case is to construct a working basis and to try to perform most linear trans-

formations with submatrices chosen from W . In this case, the procedure becomes

quite similar to the procedures for directly attacking (3.1.2) that are given in the next

section.

242 5 Two-Stage Recourse Problems

The generalized programming approach is also useful in considering the stochas-

tic program as a procedure for combining tenders

(see Nazareth and Wets

[1986]) bid from the subproblems. In this case, the method may converge most

quickly if the initial set of tenders is chosen well. A method for choosing such an

initial set of tenders appears in Birge and Wets [1984]. This view of stochastic pro-

grams can also be quite useful for stochastic integer programs and is used to obtain

efﬁciencies in branch-and-bound algorithms as discussed in Section 7.3.

Exercises

1. Suppose Problem (6.4) is infeasible for some k . What can be said about the

original two-stage stochastic linear program? Find examples for these possible

situations.

2. Prove directly that the inner linearization method converges to an optimal solu-

tion to the two-stage stochastic linear program (3.1.2).

3. Show that any extreme descending ray in (6.11) corresponds to an extreme ray

of a linear piece of

(

) .

4. Describe the steps of the generalized programming method for a modiﬁcation

of Example 2 in which the ﬁrst period costs are

x ,where

= {−2,−0.5,

0.5,1,2}. What differs in the path of the algorithm as

changes?

5.7 Simple and Network Recourse Problems

In many stochastic programming problems, special structure provides additional

computational advantages. The most common structures that allow for further efﬁ-

ciencies are simple recourse and network problems. The key features of these prob-

lems are separability of any nonlinear objective terms and efﬁcient matrix computa-

tions.

Separability is the key to simple recourse computations. In Section 3.1 and Sec-

tion 5.6, we described how these problems involve a recourse function that separates

into components for each random variable. With simple recourse, the stochastic pro-

gram in (6.5) can then be written with a separable recourse function as:

min z = c

x +

∑

i=1

(

)

s. t. Ax = b ,

Tx−

= 0 ,

x ≥0 ,

(7.1)

5.7 Simple and Network Recourse Problems 243

where

(



≤

−

(

− h

)dF(h



−

)dF(h

) . Using this

form of the objective in

, we can substitute in (3.1.9) to obtain:

(

)=q

−(q

−q

(

))

−q



≤

dF(h

) , (7.2)

where

= E[h

] .

The separable objective terms in (7.1) offer advantages for computation. In gen-

eral, we can use nonlinear programming techniques that apply even when the ran-

dom variables are continuous. Linear programming-based procedures apply as well

when the random variables have a ﬁnite number of values. In this section, we will

ﬁrst show how to use linear programming structure, assuming that each h

takes on

the values, h

i, j

, j = 1,...,K

with probabilities p

i, j

. We then consider methods for

general nonlinear problems.

Wets [1983a] gave the basic framework for computation of ﬁnitely distributed

simple recourse problems as a linear program with upper bounded variables. The

idea is to split

into values corresponding to each interval, [h

i, j

i, j+1

] ,sothat

∑

j= 0

i, j

i,0

≤ h

i,1

, 0 ≤

i, j

≤ h

i, j+1

−h

i, j

, 0 ≤

i,K

. (7.3)

The objective coefﬁcients correspond to the slope of

(

) in each of these inter-

vals. They are:

i,0

= −q

i, j

= −q

+ q



∑

l=1

i,l



, j = 1,...,K

. (7.4)

The piecewise linear program with these objective coefﬁcients and variables is:

min z = c

x +

∑

i=1



∑

j= 0

i, j



+ q



s. t. Ax = b ,

Tx−

= 0 ,

x ≥ 0and(7.3).

(7.5)

The equivalence of (7.1)and(7.5) is given in the following theorem.

Theorem 20. Problems (7.1) and (7.5) have the same optimal values and sets of

optimal solutions, (x

∗

) .

Proof: We ﬁrst show any solution (x,

,...,

) to (7.1) corresponds to a so-

lution (x,

,...,

1,1

,...,

) to (7.5) with the same objective value. We

then also show the reverse to complete the proof. Suppose (x,

) feasible in (7.1).

If h

i, j

≤

< h

i, j+1

for some 1 ≤ j ≤ K

,thenlet

i,0

= h

i,1

i,l

= h

i,l+1

−h

i,l