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

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

404 9 Monte Carlo Methods

Exercises

1. Consider Example 1 in Section 5.1. Find the projection of a point onto X = {x |

0 ≤ x ≤ 10}. Solve this problem using the stochastic quasi-gradient method

until 20 consecutive iterations are within 1% of the optimal solution.

2. Consider Example 1, where both x

and x

can be chosen instead of x =

= x

. Follow the stochastic quasi-gradient method again until 20 consecutive

iterations are within 1% of the optimal solution.

3. Prove Theorem 2.

4. Consider Example 1 in Section 5.1. Find the projection of a point onto X = {x |

0 ≤ x ≤ 10}. Solve this problem using the stochastic quasi-gradient method

until three consecutive iterations are within 1% of the optimal solution.

9.4 Sampling Methods for Probabilistic Constraints and

Quantiles

Monte Carlo sampling methods can also be quite useful for general types of proba-

bilistic constraints, such as:

P{g

(x,ξ) ≤ 0,i = 1,...,m}≥

, (4.1)

where the functions g

,i = 1,...,m are all convex in x ∈ ℜ

. In this case, a

simple procedure is to select a random independent sample of

realizations,

{

,...,

},of ξ and, with a linear objective, to solve the sample problem:

min c

x (4.2)

s. t. g

(x,

) ≤ 0,i = 1,...,m;k = 1,...,

As shown in Caliﬁore and Campi [2005], we have the following:

Theorem 4. If

≥

εβ

−1 for any

∈ (0,1 −

] and

∈ (0,1) , then with prob-

ability at least 1 −

, the solution x

to (4.2) also satisﬁes the probabilistic con-

straint (4.1).

Other results for general conditions on x , such as in de Farias and Van Roy

[2005] can also be useful in this context. To see how these approximations can work,

consider Example 8.3 for a single period where each A

is a Bernoulli random

variable that each loan has a value of one at maturity if not in default and has value

zero otherwise. To include correlation, we suppose the Merton [1974] model where

default occurs for loan j if (the natural logarithm of) its underlying asset value sat-

isﬁes the inequality,

√

1 −

≤ d

for some default point d

,where

is the correlation between any pair of underlying asset values, and ξ

and ξ

are

independent and normally distributed random variables with zero mean with unit

9.4 Sampling Methods for Probabilistic Constraints and Quantiles 405

variance. In terms of the probabilistic constraint, A

= 1

{

√

1−

}

where

1 is the indicator with value one on the given set. The probability of a default for

loan j(= 1,...,n) is p = p

= 1 −E[A

) ,where

is the standard nor-

mal cumulative distribution function. (In this development, which follows Vasiˇcek

[1987, 1991, 2002], the time to maturity is normalized to one and the information

about the drift and volatility of the underlying asset values are subsumed in d

We can evaluate the probability of satisfying the constraint with liability level h

P{Ax ≥h} = P{

∑

j= 1

≥

}, (4.3)

and then use the probability mass function of

∑

j= 1

given by (Exercise 1):

∑

j= 1

= n −k} =







∞

−∞



√

1 −

(

−1

(p) −

√





√

1 −

(

−1

(p) −

√



n−k

(s). (4.4)

Using (4.4)tosolve:

minb s. t. P{Ax ≥h}≥

,x =

, (4.5)

may not be practical for large n (e.g., when n = 125 as in the example). Instead,

we can use the sample approximation in (4.2). Exercise 2 asks for this computation

for typical default probabilities and correlations.

The sampling approximation in (4.2) can be relaxed with some allowable frac-

tion of constraint violations to achieve more precise approximations. Luedtke and

Ahmed [2008] use Hoeffding’s inequality to achieve bounds from this sampling ap-

proach. Other approximations for probabilistic constraints appear in De´ak [1980],

Gassmann [1988], Sz´antai [1986], and elsewhere. We brieﬂy describe Sz´antai’s

method here. The basic idea is to use Bonferroni-type inequalities to write the prob-

ability of a set with many constraints in terms of sums and differences of integrals of

subsets of the constraints, as we described in Section 8.5. In sampling procedures,

these alternative estimates allow for signiﬁcant variance reduction.

For Sz´antai’s approach, suppose we wish to ﬁnd

p = P[A = A

∩···∩A



dF(ξ) . (4.6)

Sz´antai takes three estimates of p :

1. ˆp

—a direct Monte Carlo sample;

2. ˆp

—ﬁnding the ﬁrst-order Bonferroni terms, 1 −

∑

i=1

) , directly and sam-

pling from higher-order terms;

406 9 Monte Carlo Methods

3. ˆp

—Calculating the ﬁrst- and second-order terms explicitly, 1 −

∑

i=1

∑

i< j

∩

) , and sampling from higher order terms.

Sampling from all higher order terms may be difﬁcult, but Sz´antai shows that the

effort may be reduced at each sample

to ﬁnding ˆn( j) deﬁned as the number

of constraints violated by

, i.e., ˆn( j)=

∑

i=1

{

∈A

}

. With this quantity deﬁned,

we can deﬁne unbiased estimates, i.e., estimates whose expectations have no error,

using the following:

∑

j= 1

max{0,1 − ˆn(j)} , (4.7)

∑

j= 1

max{0, ˆn( j) −1} , (4.8)

and

∑

j= 1

max{0, ˆn( j) −1}ˆn( j) −2)

· (4.9)

These quantities are then used to form unbiased estimates:

ˆp

, (4.10)

ˆp

= 1 −

∑

i=1

, (4.11)

ˆp

= 1 −

∑

i=1

∑

i< j

∑

∩

] −

. (4.12)

These three estimators are combined to form

ˆp

+(1 −

−

) ˆp

, (4.13)

where the weights

and

are chosen to minimize the variance of ˆp

.Theyare

calculated using the sample covariance matrix of (

) , which we denote as

C =[c

] . In this case,

,where

= c

−c

)+c

−c

)+c

−c

) , (4.14)

= c

−c

)+c

−c

)+c

123

−c

) , (4.15)

and

= c

−c

)+c

−c

)+c

−c

) . (4.16)

The result is that ˆp

can have signiﬁcantly lower variance than standard Monte

Carlo. In fact, Sz´antai obtains efﬁciencies (variance ratios) of 100 and higher, im-

plying that the same error can be obtained with ˆp

in 1% of the number of samples

for using ˆp

alone.

9.4 Sampling Methods for Probabilistic Constraints and Quantiles 407

This approach combines analytical techniques with simulation to produce lower

variance. Another approach is to use empirical sample information. This is the area

studied in Jagganathan [1985], where some sample information can be used in a

Bayesian framework to determine probabilities of underlying distributions. These

may be used for probabilistic constraints, for recourse functions, or for both.

As an example, consider the basic two-stage model in (3.1.1), where the distribu-

tion function of ξ is F(ξ,η) ,where η is a k -vector of unknown parameters with

prior distribution function, G(·) . Given an observation,

,...,

) , we can

deﬁne a posterior distribution, G

(·|

) . Using this, we may obtain an improved

solution.

Without sample information, we would have the solution to (3.1.1)as

R(G)=min

x∈K



x +



Q(x,ξ)F(ξ,η)G(dη)



. (4.17)

However, using

, which we assume has a conditional distribution given by

) for some value

of η , we obtain a value with sample information as

(G)=





min

x∈K





Q(x,ξ)F(dξ,η)G

(dη |

)



W(d

,η)G(dη) . (4.18)

The difference R

(G) −R(G) is the expected value of sample information.This

represents the additional expected value from observing the sample information.

This type of analysis can also be extended to problems with probabilistic constraints.

A different use of sample information is for dynamic problems that may change

over time. In these cases, future characteristics, such as product demand, may not be

known with certainty but they can be predicted roughly using past experience. These

problems were examined by Cipra [1991], who also considered the possibility that

more recent information might be more valuable than older information.

For example, consider the news vendor problem in Section 1.1. Suppose the

demand occurs as ξ

for periods t = 1,...,H . At time H , suppose that

(

,...,

) have been observed. The news vendor wishes to place an order based

on these observations. One solution might be to use a discount factor,

∈(0,1) ,to

choose x(H) to

min

x≥0



H−1

∑

i=0

((a −s)x +(s −r)(x −

H−i

)



. (4.19)

The solution of this problem is straightforward (Exercise 5). Alternative perspec-

tives on the value of empirical observations can also be introduced, as could

Bayesian approaches as in (4.18). For another view of decisions made over time,

refer to Jagganathan [1991].

408 9 Monte Carlo Methods

Exercises

1. Derive (4.4)(Vasiˇcek [1987]).

2. Use a sufﬁcient number of samples

from Theorem 4 for the sample approxi-

mation, (4.2), to ensure that the probabilistic constraint in (4.5) is satisﬁed with

a conﬁdence level of 1 −

= 0.99 with target reliability level,

= 0.95 ,

h = 0.95 , n = 125 , p = 0.01 for the probability of default on any single loan,

and correlation coefﬁcient

= 0.5. (Since b is the only decision parameter

to consider when all the loans are symmetric, you can assume the dimension to

use in computing

is one.) Find the minimum b

∗

for 100 different sample

problems. Verify the result in Theorem 4 empirically by constructing a sample

of 10,000 sample portfolios and solving (4.5) for this large sample. What would

you expect to happen if the problem is interpreted as making 125 separate deci-

sions x

on the initial size of loan j ?

3. In Exercise 2, you should have noticed that Theorem 4 provides an estimate

of b

∗

that appears quite conservative. An alternative approximation is to use a

limiting distribution when a portfolio contains many loans. Suppose F

is the

cumulative distribution of the fractional loss of a portfolio of n loans so that

(

∑



n

k=1

∑

j= 1

= n −k} . Substituting y for

(

√

1−

(

−1

(p) −

√

s)) and dG(y)=d

(s) , F

(

) can be written as:

(



n

∑

k=1







(1 −y)

n−k

dG(y). (4.20)

Take the limit n →∞ in (4.20)toshow:

∞

(

)=G(

), (4.21)

and ﬁnd G (Vasiˇcek [1991]). Compare this approximation to your simulation

results in Exercise 2.

4. Show that ˆp

are unbiased estimators of the probability p in (4.6).

5. Suppose that γ

, i = 1,2,3 , are independent standard gamma random variables

with parameters,

, i = 1, 2,3.Let x

= γ

+γ

i+1

for i = 1, 2 . Give a one di-

mensional integral that represents P[x

≤w

,i = 1, 2] using cumulative gamma

distribution functions in the integrand.

6. The result in Exercise 2 allows calculations of ˆp

. For example, suppose that

, i = 1, 2,3,4inExercise2andx

= y

+ y

i+1

for i = 1,2, 3.Find ˆp

for

p = P [x

≤ z

,i = 1,2,3] when z

= 6, i = 1, 2,3,and

= 3, i = 1, 2,3,4.

Also, ﬁnd sample variances for increasing sample sizes and compare to the sam-

ple variances for ˆp

7. Suppose that ξ is known to take on a ﬁnite number K of possible values but

the probabilities

of these values are not known but have a Dirichlet prior

distribution. Show how to ﬁnd R(G) and R

(G) in this case.

9.5 General Results for Sample Average Approximation and Sequential Sampling 409

8. Find the solution to (4.19). (Hint: Order the observed demands.)

9.5 General Results for Sample Average Approximation and

Sequential Sampling

We will give a brief overview of general sampling results. For this analysis, we

consider a stochastic program in the following basic form:

inf

x∈X



g(x,ξ)P(dξ) , (5.1)

where X ⊂ℜ

and ξ is now deﬁned on the probability space (

,B,P) so that we

can work directly with

instead of through

. Suppose that (5.1) has an optimal

solution, x

∗

,andvalue, z

∗

A direct sampling approach to solving (5.1) is to consider an approximate prob-

lem derived by taking

samples from ξ . The discrete distribution with these sam-

ples could be P

, which would allow us to apply the results in Chapter 9 to obtain

convergence of the

problem optimal solutions to the optimal solution in (5.1). It

can be even more valuable to describe distributional properties of these solutions so

that we can construct conﬁdence intervals in place of the (probability one) bounds

found in Chapter 8.

We therefore wish to consider a sample {

} of independent observations of ξ

that are used in the general sample average approximation (SAA) problem:

= inf

x∈X

∑

i=1

g(x,

) . (5.2)

Suppose that x

is the random vector of solutions to (5.2) with independent random

samples, ξ

, i = 1,...,

. The general question considered in King and Rockafellar

[1993] is to ﬁnd a distribution u such that

√

−x

∗

) converges to u in distri-

bution. Properties of u can then be used to derive conﬁdence intervals for x

∗

from

an observation of x

We give the main result without proof. The interested reader can refer to King

and Rockafellar [1993] and, for the statistical origin, Huber [1967].

Theorem 5. Suppose that g(·,

) is convex and twice continuously differentiable,

X is a convex polyhedron, ∇g :

×ℜ

→ ℜ

i. is measurable for all x ∈ X;

ii. satisﬁes the Lipschitz condition that there exists some a :

→ ℜ ,



|a(ξ)|

P(dξ) < ∞ , |∇g(x

) − ∇g(x

)|≤a(

)|x

− x

|, for all x

∈ X;

410 9 Monte Carlo Methods

iii. satisﬁes that there exists x ∈ X such that



|g(x,ξ)|

P(dξ) < ∞ ; and, for

∗



∇

g(x

∗

,ξ)P(dξ) ,

iv. (x

−x

)

∗

−x

) > 0 , ∀x

= x

∈X.

Then the solution x

to (5.2) satisﬁes:

√

−x

∗

) → u , (5.3)

where u is the solution to:

min

∗

u + c

s. t. A

i·

≤ 0 ,i ∈I(x

∗

) ,u

∇¯g

∗

= 0 , (5.4)

X = {x | Ax ≤ b} , (x

∗

) solve ∇



g(x

∗

,ξ)P(dξ)+(

∗

)

A = 0 ,

∗

≥ 0 ,

∗

≤ b, I(x

∗

)={i | A

i·

∗

= b

}, ∇¯g

∗



∇g(x

∗

,ξ)P(dξ) , and c is distributed

normally N(0,

∗

) with

∗



(∇g(x

∗

,ξ) −∇¯g

∗

)(∇g(x

∗

,ξ) −∇¯g

∗

)

P(dξ) .

Example 2

Suppose that X =[a,∞) , ξ is normally distributed N(0,1) ,and g(x,

)=(x −

)

. Problem (5.1) then becomes:

inf

x≥a



(x −

)

√

−

, (5.5)

where we substituted for P the standard normal density with mean zero and unit

standard deviation.

Because the expectation in (5.5)isjust x

+ 1,for a ≥ 0 , the clear solution is

∗

= a .For a < 0, x

∗

= 0 . In this case, ∇g(x

∗

)=2(x

∗

−

) , G

∗

= 2, A =

[−1] ,and ∇¯g

∗

= 2x

∗

. The variance of c is

∗

= E

[(2ξ)

]=4 . The asymptotic

distribution u then solves:

min u

+ c

s. t. u ≥ 0ifx

∗

= a ,u(2x

∗

)=0 . (5.6)

For a > 0 , the solution of (5.6)is u

∗

= 0 so that asymptotically

√

−x

∗

) →0

in distribution. If a = 0 , then note that because c/2isN(0,1) , the overall result is

that asymptotically the estimate,

√

,for(5.5) approaches a distribution with a

probability mass of 0.50 at 0 and the density of the normal distribution, N(0,1) ,

over (0,∞) . Exercise 1 asks the reader to ﬁnd the asymptotic distribution for a <

0 . In each case, the actual distribution of x

can be found and compared to the

asymptotic result (see Exercise 2).

9.5 General Results for Sample Average Approximation and Sequential Sampling 411

Many other results along these lines are possible (see, e.g.,

Dupaˇcov´a and Wets [1988]). They often concern the stability of the solutions with

respect to the underlying probability distribution. For example, one might only have

observations of some random parameter but may not know the parameter’s distri-

bution. This type of analysis appears in Dupaˇcov´a [1984], R¨omisch and Schultz

[1991a], and the survey in Dupaˇcov´a [1990].

Another useful result is to have asymptotic properties of the optimal approxi-

mation value. For this, suppose that z

∗

is the optimal value of (5.1)and z

is the

random optimal value of (5.2). We use properties of g and

so that each g(x,

)

is an independent and identically distributed observation of g(x, ξ) ,and g(x,ξ)

has ﬁnite variance, Var(g(x)) =



|g(x,ξ)|

P(dξ) −(Eg(x,ξ))

. We can thus ap-

ply the central limit theorem to state that

√





∑

i=1

g(x,

) −



g(x,ξ)P(dξ)



converges to a random variable with distribution, N(0, Var (g(x))) . Moreover, with

the condition in Theorem 5, the random function on x deﬁned by

√





∑

i=1

g(x,

) −



g(x,ξ)P(dξ)



is continuous. We can then derive the fol-

lowing result of Shapiro [1991, Theorem 3.3].

Theorem 6. Suppose that X is compact and g satisﬁes the following conditions:

i. g(x,·) is measurable for all x ∈ X;

ii. there exists some a :

→ ℜ ,



|a(ξ)|

P(dξ) < ∞ , |g(x

) −g(x

)|≤

)|x

−x

|, for all x

∈X;

iii. for some x

∈ X,



g(x

,ξ)P(dξ) < ∞ ;

and Eg(x) has a unique minimizer x

∈X.Then

√

−z

∗

] converges in distri-

bution to a normal N(0, Va r g(x

)) .

Further results along these lines are possible using the speciﬁc structure of g

for the recourse problem as in (3.1.1). For example, if K

is bounded and Q has

a strong convexity property, R¨omisch and Schultz [1991b] show that the distance

between the optimizing sets in (5.1)and(5.2) can be bounded.

Given the results in Theorems 5 and 6 and some bounds on the variances and

covariances, one can construct asymptotic conﬁdence intervals for solutions using

(5.2). We discuss this use in a sequential sampling method below. In addition, note

that all previous discrete methods can be applied to (5.2) to obtain solutions as

increases. Various procedures can be used to increment

and solving the resulting

approximation (5.2) using a previous solution.

Stronger results than Theorem 6 are possible when the minimum in (5.1)isa

sharp minimum in the following sense:

Eg(x, ξ) ≤ Eg(x

∗

,ξ)+kx −x

∗

, (5.7)

for some k > 0forallx ∈ X . In this case, with probability one, x

= x

∗

for all

sufﬁciently large, i.e., the convergence is exact, and, for two-stage stochastic linear

programs with relatively complete recourse, the rate of convergence is exponentially

fast, i.e., the probability of not converging in

iterations is bounded by

−

βν

for

some constants

> 0and

> 0 (Shapiro and Homem-de-Mello [2000]). Similar

412 9 Monte Carlo Methods

convergence results in general cases are possible using large-deviation theory such

that the probability of error in the objective value and in the solution (under certain

conditions) is greater than any tolerance decreases exponentially fast in the number

of samples. We state these results in the following theorem (see Theorems 3.1 and

3.2 in Dai, Chen, and Birge [2000] for the proof).

Theorem 7. Assume that there exist a > 0 ,

> 0 ,

(·) : ℜ

→ ℜ

such that

|g(x,

)|≤a

(

), E[e

θη

(

)

] < ∞

for all x ∈ X and for all 0 ≤

≤

; then, for any

> 0 , there are

> 0,

> 0

such that

P[E[z

−z

∗

)] ≥

] ≤

−

βν

, (5.8)

for all

> 0 , and, if x

∗

is unique,

P[||x

−x

∗

|| ≥

] ≤

−

βν

(5.9)

for all

≥ 1 .

Theorem 7 provides the possibility of some stopping criteria for sampling meth-

ods to achieve approximate optimality with some conﬁdence. Exercise 4 asks for es-

timates of the parameters

and

for Example 1. In some cases (with a quadratic

objective) discussed in Dai, Chen, and Birge [2000], these parameters can be found

explicitly (although these analytical values of the parameters often result in loose

bounds). These results provide asymptotic results that may be used in an algorithm

to obtain convergence within some tolerance of the optimal solution value with a

given level of conﬁdence. A key aspect of these procedures is that they need to

include increasingly large sample sizes to ensure that the algorithm does not ter-

minate prematurely. We give a basic algorithm from Bayraksan and Morton [2009]

that follows earlier results in Morton [1998] and Mak, Morton, and Wood [1999].

We wish to use convergence for the two-stage, sample-average linear program

with

samples. In this case, x

∗

may not be unique. In that case, we let the set

of optima be X

∗

and let x

∗

min

(x)=argmin



∈X

∗

Var[g(x, ξ) −g(x



,ξ)] .Forthetwo-

stage model, we have g(x,ξ)=c

x+{min q

y|W y = h−Tx,y ≥0}, which means

that we can write the sample-average approximation problem (5.2) as follows:

= min c

x +

∑

i=1

(5.10)

s. t. Ax = b,

x + Wy

= h

x ≥0, y ≥0,

with optimal solution (x

,...,y

) . We would like to estimate the gap to opti-

mality,

)=E[g(x,ξ)] −z

∗

, (5.11)

9.5 General Results for Sample Average Approximation and Sequential Sampling 413

and the smallest variance of the objective difference among the optimal solutions,

(x)=Var [g(x,ξ) −g(x

∗

min

(x),ξ)]. (5.12)

We then deﬁne an optimality gap estimator and its sample variance estimator as

follows:

(x)=

∑

i=1

(g(x,

) −g(x

)) (5.13)

(x)=

−1

∑

i=1

[(g(x,

) −g(x

)) −G

(x)]

. (5.14)

The goal in sequential sampling is to obtain, for any given conﬁdence level

∈

(0,1) and tolerance

> 0, x

after

samples such that

liminf

l↓l



P(E[g(x

,ξ) −g(x

∗

,ξ)] ≤lb+



) ≥ 1 −

(5.15)

for some parameters l > l



> 0, b > 0and



. For deﬁning an algorithm, we

use additional parameters: k

, the frequency of re-sampling, and p > 0,whichis

used in determining a minimum sample size for k iterations as follows:

≥



l −l







max



2ln



∞

∑

j= 1

−pln j

√

πα



+ 2pln



(5.16)

(following Bayraksan and Morton [2009] who use p ≈2×10

−1

≈2×10

−7



≈

−7

,l ≈ 0.045,l



≈ 0.015 ). There are also two sequences of sample numbers:

for checking optimality and

(e.g.,

= 2

) for choosing the next candidate.

The candidate solution is x

that solves (5.10) with

samples.

Sequential Sampling Method (SSM)

Step 0. Initialize with k = 1,

from (5.16).

Step 1. Generate

samples to obtain ˆx

= x

. (These can start with the pre-

vious

k−1

samples to use the previous solution as a starting point, but they are

independent from the

samples for gap estimates.)

Step 2. Generate

samples (IID) to form G

= G

( ˆx

) and s

= s

( ˆx

) .

Step 3.If s

> b or G

> l



b +



, then: set k = k + 1, ﬁnd new

,re-sample

if k is a multiple of k

, and return to Step 1. Else, x

= ˆx

(with, for

= 2

∑

i=1

total samples, including re-samples).

The convergence result for this method is contained in the following theorem,

where P refers to the probability measure over the sampling distribution.