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

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

414 9 Monte Carlo Methods

Theorem 8. For a two-stage stochastic linear program (3.1.1) with relatively com-

plete recourse, almost surely ﬁnite second-stage value, and compact non-empty fea-

sible region X ,

1. for



> 0 ,p> 0 , and 0 <

< 1 ﬁxed values, if the method stops at

iteration

,then

liminf

l↓l



) ≤ls

) ≥1 −

; and, (5.17)

2. for ﬁxed



> 0 and l > l



> 0 where the sequential sampling method stops at

iteration

, P (

< ∞)=1 .

Proof: The proof follows Bayraksan and Morton [2009], Theorem 3 and Proposi-

tion 4, with additional observationsabout characteristics of the estimates for sample-

average approximations of two-stage stochastic linear programs (see, e.g, R¨omisch

[2003]).

As a ﬁnal note, we should mention that analogous procedures can be built around

quasi-random sequences that seek to ﬁll a region of integration with approximately

uniformly spaced points. The result is that errors are asymptotically about of the

order log(

instead of 1/

√

(see Niederreiter [1978]). The difﬁculty is in

the estimation of the constant term but quasi-Monte Carlo appears to work quite

well in practice (see Fox [1986] and Birge [1994]). In terms of expected perfor-

mance over broad function classes, quasi–Monte Carlo performs with the same or-

der of complexity (Wo´zniakowski [1991]). For the methods used in this chapter, we

may substitute quasi-random sequences for pseudo-random sequences for practical

implementations. Other generalizations known as sparse grid (from Smolyak

[1960]) can also be used to obtain efﬁcient characterizations of the integrals in

stochastic programs (see Chen and Mehrotra [2007]).

The sample average approximation method has also been used for stochastic in-

teger programs. An SAA problem for an SIP is similar to the program (1.2), with

integer requirements in the ﬁrst- and/or the second-stage programs. An SAA with

a moderate sample size can be solved using classical deterministic MIP techniques.

The process can be repeated with different samples to obtain candidate solutions

along with statistical estimates of their optimality gaps. These various candidate

solutions cannot be combined (or averaged) as they would produce non-integer so-

lutions. Instead, a new and independent large sample is created to form an estimated

objective function. The various ﬁrst-stage candidate solutions are evaluated using

this estimated objective solution. These evaluations still require several second-stage

optimizations. The computational burden remains low as the ﬁrst-stage is given. At

the end, the best candidate ﬁrst-stage solution is selected. A detailed computational

study of the application of the SAA method to solve three classes of stochastic rout-

ing problems can be found in Verweij et al. [2003].

9.5 General Results for Sample Average Approximation and Sequential Sampling 415

Exercises

1. For Example 2, ﬁnd the asymptotic result from Theorem 5 for

√

−x

∗

) for

a < 0.

2. For Example 2, derive the actual distribution of

√

−x

∗

) for a feasible

region x ≥ a in each case of a , a < 0, a = 0anda > 0 . Find the limits of

these distributions and verify the result from Theorem 5.

3. Consider a news vendor problem as in Section 1.1. Suppose this problem is

solved using a sampling approach. The sampled problem with continuous cu-

mulative distribution function F

has a solution at (F

)

−1



s−a

s−r



= x

.Find

the distribution of this quantile and show how to construct a conﬁdence interval

around x

∗

4. Consider Example 1. First, verify the assumptions in Theorem 7. Solve 100

samples each for

= 10 + 10i for i = 0,1,...,10 . Use these observations to

estimate values for

and

in Theorem 7 for

= 0.03 .

5. Implement the sequential sampling method for the continuous distribution ver-

sion of the two-stage stochastic linear program for the farming example in Sec-

tion 1.1. Start with the parameter recommendations above. Vary them to observe

the impact of the parameters on the convergence behavior.

Chapter 10

Multistage Approximations

Most decision problems involve effects that carry over from one time to another.

Sometimes, as in the power expansion problem of Section 1.3, random effects can be

conﬁned to a single period so that recourse is block separable. In other cases, how-

ever, this separation is not possible. For example, power may be stored by pumping

water into the reservoir of a hydroelectric station when demand is low. In this way,

decisions in one period are inﬂuenced by decisions in previous periods.

Problems with this type of linkage among periods are the subject of this chapter.

We again wish to derive approximationsthat can be used to bound the error involved

in any decision based on the approximate problem solution. In Chapter 9,wesaw

that the number of random variables can lead to rapidly growing problems. In this

chapter, we have the additional effect that the number of periods leads to exponential

increases in problem size even if the number of realizations in each period remains

constant (see Figure 3.4).

We can again construct bounds based on the properties of the multistage recourse

functions. These analogues of the basic Jensen and Edmundson-Madansky bounds

are given in Section 10.1. They correspond to ﬁxing values at means or extreme

values of the support of the random vectors in each period.

Keeping the number of periods ﬁxed may not lead to sufﬁcient reductions in

problem size, especially if no time is clearly the end of the problem. This case would

mean facing either an uncertain or an inﬁnite horizon decision problem. These prob-

lems can also be approximated by aggregating several periods together. Section 10.2

describes this procedure to obtain both upper and lower bounds.

Sampling methods that apply generally to multistage methods are described in

Section 10.3. Section 10.4 then describes methods based on decomposition ap-

proaches.

The bounds of Sections 10.1 and 10.2 and the sampling methods used in in Sec-

tions 10.3 and 10.4 can be viewed as discretization procedures. We can also con-

struct separable bounds that do not require discretization as in Chapter 8.These

bounds correspond to separable responses to any changes in the problem and are

part of a general approach to value-function approximation known as approximate

dynamic programming. They are described in Section 10.5. In multistage problems,

J.R. Birge and F. Louveaux, Introduction to Stochastic Programming, Springer Series 417

in Operations Research and Financial Engineering, DOI 10.1007/978-1-4614-0237-4

10,

 Springer Science+Business Media, LLC 2011

418 10 Multistage Approximations

speciﬁc problem forms and structures can also lead to substantial savings. These

structures are particularly valuable for approximations of the value function. We

also describe such special cases for network revenue management, production, and

vehicle allocation in this concluding section.

10.1 Extensions of the Jensen and Edmundson-Madansky

Inequalities

The basic Jensen and Edmundson-Madansky inequalities can be extended to mul-

tiple periods directly. The principle is to use Jensen’s inequality (or a feasible dual

solution) to derive the lower bound and construct a feasible primal solution using

extreme points to construct the upper bound. To present these results, we consider

the linear case ﬁrst, although extensions to nonlinear, convex problems are directly

possible. We use concepts from measure theory in the following discussion. Readers

without this background may skip to the declarations to ﬁnd the major results for

actual implementations.

The multistage stochastic linear program is to ﬁnd x =(x

,...,x

) (where

we suppress transposes as earlier when they can be implied from the context) in the

following:

min c

+ E

+ ···+ c

]

s. t. W

= h

t−1

= h

, t = 2,...,H , a.s.,

−E

]=0 , t = 2,...,H , a.s.,

≥ 0 , t = 1,...,H , a. s.,

(1.1)

where we have used explicit nonanticipativity constraints as in (3.5.11). We have

also assumed that the recourse within each period W

is known and not random.

The basic Jensen bound again follows by assuming a partition of

, the support

vector of all random components. Here, we write

×···×

.We

suppose that

= {

,...,

) |

∈

,i = 1,...,t}. In this way, we can

characterize all events up to time t by measurability with respect to the Borel ﬁeld

deﬁned by

. We assume that

is partitioned as

= S

∪···∪S

and

that S

= ∪

j∈D

t+1

(i)

{

| (

t+1

) ∈ S

t+1

} so that the partitions are consistent

from one period to another. We construct measurable decisions at time t if they are

constant over each S

t, j

Next, assume that p

= P[S

] , c

= c

,andthat E

[(h

)] = (

) for all t

and i . The problem then is to ﬁnd:

min c

∑

i=2

+ ···+

∑

i=1

10.1 Extensions of the Jensen and Edmundson-Madansky Inequalities 419

s. t. W

= h

t−1

, t = 2,...,H , i = 1,...,

t−1

j ∈ D

t+1

(i) ,

≥ 0 , i = 1,...,

, t = 1,...,H .

(1.2)

The ﬁrst result is that (1.2) provides a lower bound on the optimal solution in (1.1)

provided the expectations of (

) are independent of the past. If not, then the

conditional expectation form in (1.2) may not actually achieve a bound.

Theorem 1. Given that E

[(h

)] = (

)=E

[(h

)] for all S

and S

that have a common outcome at time t , i.e., such that (

t−1

) ∈ S

if and only

if there exist some (

t−1

) ∈ S

. The optimal value of (1.2) with the deﬁnitions

given earlier provides a lower bound on the optimal value of (1.1).

Proof: Suppose an optimal solution x

∗

to (1.2) with optimal dual variables

t∗

corresponding to constraints in (1.2) with right-hand sides,

.By dual feasibility in

(1.2),

≥

t∗

∑

j∈D

t+1

(i)

t+1∗

t+1

, (1.3)

for every (t, i) .Let

(

∑

i=1

{

∈S

}

[

t∗

] .We also have

(

)=−

∑

i=1

{

∈S

}

[

t∗

(

)/p

∑



∈D

t−1

(i))

[

t∗



t−1,A

t−1

(i)

Note how the

variables represent nonanticipativity. The condition for dual feasi-

bility from the multistage version of Theorem 3.13 (see Exercise 1) is that

(

) −

(

−

t+1

(

t+1

(

) −

t+1

(

) ≥0 , a.s., (1.4)

and

[

t+1

(

)] = 0 . (1.5)

Substituting in the right-hand side of (1.4) yields:

−(

t∗

−[

t+1∗

t+1

(

)



[

t+1∗

t+1

(

) −

∑

j|j∈D

(i)

[

t+1∗

t+1

]



(1.6)

for each S

and j ∈ D

(i) , which is non-negative from (1.3). Also, by their def-

inition and the assumption that integration of T

t+1

(

) over varying S

does not

change its conditional outcome,

420 10 Multistage Approximations

[

t+1

(

)] =

∑

j∈D

(i)



(

t+1∗

t+1

)

t+1

−

∑

j|j∈D

(i)

(

t+1∗

t+1

)



= 0

(1.7)

yielding (1.5). Hence, we have constructed a dual feasible solution whose objective

value is a lower bound on the objective value of (1.1) by the multistage version

of Theorem 3.13. Because this value is the same as the optimal value in (1.2), we

obtain the result.

Thus, lower bounds can be constructed in the same way for multistage problems

as for two-stage problems, provided the data have serial independence. Such inde-

pendence is not necessary if only right-hand sides vary because the dual feasibility

is not affected in that case. The key procedure is in developing a dual feasible so-

lution (lower bounding support function). Upper bounds can follow as before by

constructing primal feasible solutions. These bounds can also be used in conjunc-

tion with the lower bounds to obtain bounds when objective coefﬁcients ( c

)are

also random.

To develop the upper bounds, the basic result is an extension of Theorem 8.2. We

assume the following general form in which the decision variables x are explicit

functions of the random outcome parameters, ξ :

inf

x∈N



∑

t=0

(ξ),x

t+1

(ξ),ξ

t+1

)



, (1.8)

where we use the convention for the general nonlinear objective formulation that

subscript t corresponds to decisions or parameters within period t while super-

script t refers to all periods from 1 to t . The random vector ξ =(ξ

,...,ξ

) has

an associated probability space, (

,P) , N is the space of nonanticipative deci-

sions, f

is convex, and ξ

t+1

is measurable with respect to

t+1

and ξ

t+1

∈

t+1

which is compact, convex, and has extreme points, ext

t+1

, with Borel ﬁeld, E

t+1

In this representation, x nonanticipative means that x

(

)) is

-measurable

for all t . It could also be described in terms of measurability with respect to

the Borel ﬁeld deﬁned by the history process ξ

=(ξ

,...,ξ

) .

Suppose that e =(e

,...,e

)

where each e

∈ ext

. The set of all such ex-

treme points is written ext

. Suppose x



=(x



,...,x



) ,where x



: ext

→ ℜ

We say that x



is extreme point nonanticipative,or x



∈ N



,if x



is measur-

able with respect to the Borel ﬁeld, E

,on ext

,deﬁnedby (e

,...,e

) ,where

∈ ext

(for t = 1 , this will be with respect to {/0,ext

}). With these deﬁni-

tions, we obtain the following result.

Theorem 2. Suppose that

→ f

t+1

) is convex for t = 0,...,H,

compact, convex, and has extreme points, ext

. For all

∈

,let

(

,·) be a

probability measure on (ext

,E ) where E is the Borel ﬁeld of ext

, such that



e∈ext

(

,de)=

, (1.9)

10.1 Extensions of the Jensen and Edmundson-Madansky Inequalities 421

and

→

(

,A) is measurable with respect to

for all A ∈ E

. Then there

exists, x ∈N , such that x

(



e∈ext



(e)

(

,de) ,



∑

t=0

t+1

,ξ

t+1

)



≤



e∈ext

∑

t=0

((x



)



t+1

)

(de) , (1.10)

where x



is extreme point nonanticipative and

is the probability measure on E

deﬁned by

(A)=



(

,A)P(d

) . (1.11)

Proof: We must ﬁrst show that x as deﬁned in the theorem is nonanticipative, or

that x

(

) is

-measurable. This follows because x



(e) is E

-measurable, and,

for any A ∈E

(

,A) is

-measurable. Because each f

is convex, for any

(

),x

t+1

(

t+1

)

= f





e∈ext



)

(e)

(

,de),



e∈ext



t+1

(e)

(

,de),



e∈ext

t+1

(

,de)



≤



e∈ext

((x



)

(e),x



t+1

(e),e

t+1

)

(

,de) . (1.12)

Integrating with respect to P , the result in (1.10) is obtained.

As in Chapter 8, we implement the result in Theorem 2 by ﬁnding an appropriate

and then solving the following approximation problem.

inf

x∈N





ext



∑

t=0

(e),x

t+1

(e),e

t+1

)



(de) (1.13)

to ﬁnd an upper bound on the value in (1.8). One can also reﬁne these bounds by

taking partitions of

The simplest type of bound from Theorem 2 is the extension of the Edmundson-

Madansky bound on rectangular regions with independent components. For this

bound, we assume that all components,

(i) , are stochastically independent and

distributed on [a

(i),b

(i)] . In this case, we can deﬁne

EM−I

(

,e)=

∏

t=1

∏

i=1

(i) −e

(i)|

(i) −a

(i))

, (1.14)

so that

EM−I

(e)=

∏

t=1

∏

i=1

(i) −e

(i)|

(i) −a

(i))

. (1.15)

It is easy to check that this

meets the nonanticipative measurability requirements.

Problem (1.13) now can be written as:

422 10 Multistage Approximations

inf



∑

t=0



∑

···

t+1

∑

t+1



t+2

∑

t+2

+···+

H+1

∑

H+1

,...,e

H+1

)



,...,i

),x

t+1

,...,i

t+1

),e

t+1

)



···



, (1.16)

where x

,...,i

) corresponds to the t th-period decision depending on the out-

comes in extreme point combination e

from each period s = 1,...,H . This places

the nonanticipativity back into the problem implicitly.

Example 1

To see how this bound might be implemented, consider Example 1 in Section 6.1.

Suppose that demand is uniformly and independently distributed on [1,3] in each

period. In this case, we obtain a decision vector (x

) in period t for scenario

s = 2

+ i

for i

and i

in {1,2}. Problem (7.1.7) is, therefore, the upper bound-

ing problem (1.16) for this uniform distribution case, yielding an upper bound of

6.25 . In this case, the lower bound using the expected demand value of two in each

period is three. In Exercise 2, you are asked to reﬁne these bounds until they are

within 25% of each other.

Other extreme point combinations are clearly also possible in multiperiod prob-

lems as they are in single-period problems. Extensions to dependent random vari-

ables and f

concave in some arguments can also be made.

The bounds given in this section so far apply only to ﬁxed numbers of periods.

When periods are combined, we call the resulting problem an aggregated problem.

These problems are described in the next section.

Exercises

1. Consider the multistage stochastic linear program in the form of (1.1). Prove the

multistage version of Theorem 3.13.

2. Reﬁne the extreme point (Edmundson-Madansky) and conditional expectation

(Jensen) bounds on partitions for Example 1 from Section 6.1 until the upper

bound is within 25% of the lower bound.

10.2 Bounds Based on Aggregation

The main motivation for aggregation bounds is to deal with problems with many

(perhaps an inﬁnite number of) periods by combining periods to obtain a simpler

10.2 Bounds Based on Aggregation 423

approximate problem with fewer periods. The basic procedures in this chapter ap-

pear in Birge [1985a] and Birge [1984]. They follow the general aggregation results

in Zipkin [1980a, 1980b]. Similar methods, especially for dealing with inﬁnite hori-

zon problems, appear in Grinold ([1976, 1983, 1986]). Generalizations appear in

Wright [1994] and Kuhn [2008].

To derive both upper and lower bounds in this framework, we consider a special

form for the multistage problem in (3.4.1). We allow feasibility by adding a penalty

variable y

that can achieve feasibility in each period. This notion of model robust-

ness is quite common, although the penalty parameter q may be quite high. The

form of the multistage stochastic linear program in this case is:

minz = c

+ E

[

∑

t=2

t−1

(ξ

,...,ξ

)+q

(ξ

,...,ξ

))]

s. t. Wx

≥ h

t−1

(ξ

,...,ξ

t−1

)+Wx

(ξ

,...,ξ

)+y

(ξ

,...,ξ

) ≥ξ

t = 2,...,H ,

≥ 0; x

(

) ≥ 0 , a.s., t = 2,...,H ,

(

) ≥0 , a.s., t = 2,...,H ,

(2.1)

where superscript t again represents the variables or parameters for period t (i.e.,

not the full history), c is a known vector in ℜ

, h

is a known vector in ℜ

(

)=h

(

) is a random m -vector deﬁned on (

,P) (where

⊂

t+1

)

for all t = 2,...,H ,and T and W are known m ×n matrices. We also suppose

that

is the support of

. The parameter

is a discount factor.

Note that in (2.1), we assume that the parameters T , W , c ,and q are all con-

stant across time (with objective coefﬁcients varying only with the discount factor).

This assumption is basically made to simplify the following presentation. Varying

parameters are possible with little additional work.

The key observation for these bounds is that an optimal solution in (2.1)isno

lower than

+ E



∑

t=2

(

(ξ

,...,ξ

))



(2.2)

for any (

,...,

(ξ

,...,ξ

),...,

(ξ

,...,ξ

)) ≥ 0 a.s. that satisﬁes

(

)

W + E

[

(ξ

)]

T ≤ c

(ξ

,...,ξ

)

W + E

ξ|(ξ

,...,ξ

[

t+1

(ξ

,...,ξ

t+1

)]

T ≤

t−1

t = 2,...,H −1 ,

(ξ

,...,ξ

)

W ≤

H−1

(ξ

,...,ξ

)

W ≤

H−1

. (2.3)

You are asked to show that (2.2) subject to (2.3) provides a bound in Exercise 1.