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

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

254 5 Two-Stage Recourse Problems

infz = f

(x)+E

[ f

(y(

)] (8.3)

s. t. g

(x) ≤ 0 ,i = 1,...,m

(x,y(

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

, a. s.

If we let (

) be a multiplier vector associated with the constraints, then we can

form a dual problem to (8.3)as:

max

(

)≥0

w =

(

) , (8.4)

where

(

)=inf

x,y

z = f

(x)+E

[ f

(y(

)]

+ E

[

∑

i=1

(

)

(x,y(

))] (8.5)

s. t. g

(x) ≤ 0 , i = 1,...,m

We show duality in the ﬁnite distribution case in the following theorem.

Theorem 21. Suppose the stochastic nonlinear program (8.1) with all functions

convex has a ﬁnite optimal value and a point strictly satisfying all constraints,

and suppose

= {1,...,K} with P{

= i} = p

.Then z≥ w for every feasi-

ble x,y

,...,y

in (8.1)–(8.2) and

,...,

feasible in (8.4), and their optimal

values coincide, z

∗

= w

∗

Proof: From the general optimality conditions (see, e.g., Bazaraa and Shetty [1979,

Theorem 6.2.1]), the result follows by noting that we may take x satisfying the ﬁrst-

period constraints as a general convex constraint set X so that only the second-

period constraints are placed into the dual. We also divide any multipliers on the

second-period constraints in (8.3)by p

if they correspond to

= i .Inthisway,

the expectation over

in (8.5) is obtained.

Now, we can follow a dual ascent procedure in (8.4). This takes the form of a

subgradient method. We note that

∂θ

(

)=co{(

,...,

)

,...,(

,...,

)

} , (8.6)

where again “ co ” denotes the convex hull,

= g

( ¯x, ¯y

,k) , (8.7)

and ( ¯x, ¯y

,..., ¯y

) solves the problem in (8.5)given

. This again follows

from standard theory as in, for example, Bazaraa and Shetty [1979, Theorem 6.3.7].

We can now describe a basic gradient method for the dual problem. For our

purposes, we assume that (8.5) always has a unique solution.

5.8 Methods Based on the Stochastic Program Lagrangian 255

Basic Lagrangian Dual Ascent Method

Step 0.Set

≥ 0,

= 0 and go to Step 1.

Step 1.Given

in (8.5), let the solution be (x

,...,y

) .Let

= 0if

= 0andg

,k) ≤ 0,and

= g

,k) ,otherwise.If

= 0for

all k , stop.

Step 2.Let

maximize

(

) over

≥ 0,

≥ 0. Let

+ 1,andgotoStep1.

Assuming the unique solution property, this algorithm always produces an as-

cent direction in

. The algorithm either converges ﬁnitely to an optimal solution

or, assuming a bounded set of optima, produces an inﬁnite sequence with all limit

points optimal (see Exercise 1). For the case of multiple optima for (8.5), some

nondifferentiable procedure must be used. In this case, one could consider ﬁnding

the maximum norm subgradient to be assured of ascent or one could use various

bundle-type methods (see Section 5.9).

The basic hope for computational efﬁciency in the dual ascent procedure is that

the number of dual iterations is small compared to the number of function evalua-

tions that might be required by directly attacking (8.1)and(8.2). Substantial time

may be spent solving (8.2) but that should be somewhat easier than solving (8.1)be-

cause the linking constraints appear in the objective instead of as hard constraints.

Overall, however, this type of procedure is generally slow due to our using only

a single-point linearization of

. This observation has led to other types of La-

grangian approaches to (8.1) that use more global or second-order information.

Rockafellar and Wets [1986] suggested one such procedure for a special case of

(8.5)where f

(x)=c

x +

Cx and y(

) can be eliminated so that the second

and third objective terms become

(

,x) and the dual problem in (8.4)isthen

max

≥0

inf

{x|g

(x)≤0}

x +

Cx+

(

,x)] . (8.8)

Their approach is not to restrict the search to a single search direction but to al-

low optimization over a low dimensional set. Implementation of this method, called

the Lagrangian ﬁnite-generation method for linear-quadratic stochastic programs,

is described in King [1988a] and its application to solve practical water manage-

ment problems concerning Lake Balaton in Hungary appears in Somly´ody and Wets

[1988].

A similar method based on inner linearization approaches in nonlinear program-

ming is restricted simplicial decomposition (Ventura and Hearn [1993]). This proce-

dure replaces the line search in the Topkis-Veinott [1967] feasible direction method

with a search over a simplex. The ﬁnite generation algorithm is analogously an en-

hancement over basic Lagrangian dual ascent methods that consider only gradient

or subgradient steps. Both the ﬁnite-generation and restricted simplicial decompo-

sition methods tend to avoid the zigzagging behavior that often occurs in methods

based on single-point linearizations.

256 5 Two-Stage Recourse Problems

Another method for accelerating convergence is to enforce strictly convex terms

in the objective. Rockafellar and Wets discussed methods for adding quadratic terms

to the matrices C and D(

) so that these matrices become positive deﬁnite. In

this way, the ﬁnite generation method becomes a form of augmented Lagrangian

procedure. We next discuss the basic premise behind these procedures.

In an augmented Lagrangian approach, one generally adds a penalty

rg

( ¯x, ¯y

,k))



(

) and performs the iterations including this term. The

advantage (see the discussion in Dempster [1988]) is that Newton-type steps can be

applied because we would obtain a nonsingular Hessian. The result should gener-

ally be that convergence becomes superlinear in terms of the dual objective without

a signiﬁcantly greater computational burden over the Lagrangian approach.

The computational experience reported by Dempster suggests that few dual it-

erations need be used but that a more effective alternative was to include explicit

nonanticipative constraints as in (3.5.4) and to place these constraints into the ob-

jective instead of the full second-period constraints. In this way,

becomes



(

)=infz = f

(x)+

∑

k=1

[ f

,k)]

∑

k=1

[

(x −x

)+r/2x −x



] (8.9)

s. t. g

(x) ≤ 0 , i = 1,...,m

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

k = 1,...,K .

Notice how in (8.9) the only links between the nonanticipative x decision and the

scenario k decisions are in the (x−x

) objective terms. Dempster suggests solving

this problem approximately on each dual iteration by iterating between searches

in the x variables and search in the x

variables. In this way, the augmented

Lagrangian approach of solving (8.9) to ﬁnd a dual ascent Newton-type direction

achieves superlinear convergence in dual iterations. The only problem may come in

the time to construct the search directions through solutions of (8.9).

This method also resembles the progressive hedging algorithm of Rockafellar

and Wets [1991]. This method achieves a full separation of the individual scenario

problems for each iteration and, therefore, has considerably less work in each itera-

tion; however, the number of iterations as we shall see, may be greater. The method

can offer many computational advantages, particularly for structured problems (see

Mulvey and Vladimirou [1991a]). The key to this method’s success is that individual

subproblem structure is maintained throughout the algorithm. Related implemen-

tations by Nielsen and Zenios [1993a, 1993b] on parallel processors demonstrate

possibilities for parallelism and the solution of large problems.

The basic progressive hedging method begins with a nonanticipative solution ˆx

and a multiplier

. The nonanticipative (but not necessarily feasible) solution is

used in place of x in (8.9). The ﬁrst-period constraints are also split into each x

5.8 Methods Based on the Stochastic Program Lagrangian 257

In this way, we obtain a subproblem:

infz =

∑

k=1

[ f

)+ f

,k)+

− ˆx

)+r/2x

− ˆx



]

s. t. g

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

, k = 1,...,K ,

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

, k = 1,...,K .

(8.10)

Now (8.10) splits directly into subproblems for each k so these can be treated sep-

arately.

Supposing that (x

) solves (8.10). We obtain a new nonanticipative de-

cision by taking the expected value of x

as ˆx

and step in

+(x

− ˆx

) .

The steps then are simply stated as follows.

Progressive Hedging Algorithm (PHA)

Step 0. Suppose some nonanticipative x

, some initial multiplier

,and r > 0.

Let

= 0.GotoStep1.

Step 1.Let (x

) for k = 1,...,K solve (8.10). Let ˆx

=(ˆx

+1,1

, ...,

ˆx

+1,K

)

where ˆx

+1,k

∑

l=1

+1,l

for all k = 1,...,K .

Step 2.Let

+ r(x

+1,k

− ˆx

) .If ˆx

= ˆx

and

, then,

stop; ˆx

and

are optimal. Otherwise, let

+ 1 and go to Step 1.

The convergence of this method is based on Rockafellar’s proximal point method

[1976a]. The basis for this approach is not dual ascent but the contraction of the pair,

( ˆx

) , about an optimal point. The key is that the algorithm mapping can be

described as (

/r)=(I −V)

−1

(

/r) ,whereV is a maximal

monotone operator and

is the diagonal matrix of probabilities corresponding to

and

,i.e,where

(k−1)n

+i,(k−1)n

= p

for i = 1,...,n

and k = 1,...,K .

To describe this approach we ﬁrst deﬁne a maximal monotone operator at V

(see Minty [1961] for more general details) such that for any pairs (w,z) where

z ∈V(w) and (w



) for z



∈V (w



) ,wehave

(w−w



)

V(z −z



) ≥ 0 . (8.11)

The key point here is that if we have a Lagrangian function l(x,y) that is convex in

x and concave in y , then the subdifferential set of l(x,y) at ( ¯x, ¯y) deﬁned by

{(

) |

(x − ¯x)+l( ¯x, ¯y) ≤ l(x, ¯y), ∀x ;

(y − ¯y)+l( ¯x, ¯y) ≥ l( ¯x,y), ∀y} (8.12)

yields a maximal monotone operator by

V( ¯x, ¯y)={(

)} (8.13)

258 5 Two-Stage Recourse Problems

for (

,−

) ∈

∂

l( ¯x, ¯y) (Exercise 3).

The second result that follows for maximal monotone operators is that a contrac-

tion mapping can be deﬁned on it by taking (I −V)

−1

(x,y) to obtain (x



) ,or,

equivalently, where (x



−x,y



−y) ∈ V(x



) . The contraction result (Exercise 4)

is that, if V is maximal monotone, then, for all (x



)=(I −(1/r)V)

−1

(x,y) and

( ¯x



, ¯y



)=(I −V)

−1

( ¯x, ¯y) ,

(x



− ¯x



− ¯y



)

≤ (x − ¯x,y − ¯y)



− ¯x



− ¯y



) . (8.14)

These results then play the fundamental role in the following proof of convergence.

Theorem 22. The progressive hedging algorithm, applied to (8.1) with the same

conditions as in Theorem 14, converges to an optimal solution, x

∗

, (or termi-

nates ﬁnitely with an optimal solution) and, at each iteration that does not terminate

in Step 2,

(

ˆx

/r) −(

∗

/r)< (

ˆx

/r) −(

∗

/r) . (8.15)

Proof: As stated, the key is to ﬁnd the associated Lagrangian and to show that the

iterations follow the mapping as in (8.14). For the Lagrangian, deﬁne

l( ¯x,

)=inf

(1/r)z(x)+

x (8.16)

s. t. J

x − ¯x = 0 ,

where z(x) is deﬁned as

∑

k=1

[ f

)+Q(x

,k)] for feasible x

and as +∞ oth-

erwise,

is deﬁned as the diagonal probability matrix, and J is the matrix corre-

sponding to column sums, J

r,s

equal one if r (mod n

)= s (mod n

) and zero

otherwise. We want to show that

(

( ˆx

− ˆx

),(

−

)/r) ∈

∂

ˆx

/r);

so, we can use the contraction property in (8.14) from the maximal monotone oper-

ator deﬁned on

∂

ˆx

/r) .

Note that, for ¯x =

ˆx

and

/r =

∑

i=1

−ˆx

) ,¯x

= ˆx

(

∑

i=1

−

ˆx

)) = (x



)

(

∑

i=1

− ˆx

)) for (x



)

=(1/K) ˆx

. Because J

= ˆx

,we

have ¯x

= 0 . We can thus add the term, ¯x

to the objective in (8.16) without

changing the problem. We then obtain:

∈

∂

l( ¯x,

) ⇔−

∈ (1/r)

∂

−1

(−

)+ ¯x)+

, (8.17)

where J

(

−1

(−

)) = ¯x and

is some multiplier. For

∂

¯x

l( ¯x,

) ,

= −

and some

∈

∂

¯x

l( ¯x,

) ⇔

−

∈ (1/r)

∂

Z(x



) , (8.18)

for some J



= ˆx . We combine (8.17)and(8.18) to obtain that (

) ∈

∂

l( ¯x,

)

5.8 Methods Based on the Stochastic Program Lagrangian 259

−

∈ (1/r)

∂

−1

(−

)+ ¯x) . (8.19)

We wish to show that

( ˆx

− ˆx

) −

Πρ

/r ∈ (1/r)

∂

−1

(

−

)/r + ˆx

) . (8.20)

From the algorithm,

−

Πρ

∈

∂

z(x

)+r

− ˆx

) . (8.21)

Substituting,

+ r(x

− ˆx

) , we obtain from (8.21),

−

Πρ

+ r

− ˆx

) ∈

∂

z(x

)+r

− ˆx

) , (8.22)

which, after eliminating r

from both sides, coincides with (8.20).

By the nonexpansive property, there exists (

∗

/r) , a ﬁxed point of this

mapping. By substituting into (8.14), with (

∗

/r)=(I −V)(

∗

/r) and

(

ˆx

/r)=(I −V)(

ˆx

/r) , we have (Exercise 5):

(

ˆx

/r) −(

∗

/r) < (

ˆx

/r) −(

∗

/r) . (8.23)

Our result follows if (x

∗

) is indeed a solution of (8.1). Note that in this case,

we must have 0 = x

− ˆx

= x

− ˆx

; so, from (8.21), −

Πρ

∗

∈

∂

z(x

∗

) .

From Theorem 3.2.5, optimality in (8.1) is equivalent to

∈

∂

z(x

∗

) for some

,where J

Πρ

= 0 , which is true because J

(−

∗

)=−

∑

−x

)=0.

Hence, we obtain optimality. The method converges as desired.

We note that Rockafellar and Wets obtained these results by deﬁning an inner

product as 

,x =

x and using appropriate operations with this deﬁnition.

They also show that, in the linear-quadratic case, the convergence to optimality is

geometric.

Variants of this method are possible by considering other inner products and

projection operators. For example, we can let

¯x

be the standard orthogonal pro-

jection of x

into the null space of J

. This value is the simple average of x

values, so that

¯x

(i)=(1/K)

∑

k=1

(i) for all k = 1,...,K . The multiplier

update is then:

+ r

−1

−

¯x

) . (8.24)

One can again obtain the maximal monotone operator property, and, observing that

= J

¯x

, obtain J

Πρ

∗

= 0 and optimality.

Example 3

The algorithm’s geometric convergence may require many iterations even on small

problems as we show in the following small example. Suppose we can invest

260 5 Two-Stage Recourse Problems

$10,000 in either of two investments, A or B. We would like a return of $25,000,

but the investments have different returns according to two future scenarios. In the

ﬁrst scenario, A returns just the initial investment while B returns 3 times the initial

investment. In the second scenario, A returns 4 times the initial investment and B

returns twice the initial investment.The two scenarios are considered equally likely.

To reﬂect our goal of achieving $25,000, we use an objective that squares any return

less than $25,000. The overall formulation is then:

min z = 0.5(y

+ y

)

s. t. x

+ x

≤ 10 ,

+ 3x

+ y

≥ 25 ,

+ 2x

+ y

≥ 25 ,

≥ 0 .

(8.25)

Clearly, this problem has an optimal solution at x

∗

= 2.5andx

∗

= 7.5 with an

objective value z

∗

= 0 . A single iteration of Step 1 in the basic Lagrangian method

is all that would be required to solve this problem for any positive

value. A

single iteration is also all that would be necessary in the augmented Lagrangian

problem in (8.9). The price for this efﬁciency is, however, the incorporation of all

subproblems into a single master problem. Progressive hedging on the other hand

maintains completely separate subproblems. We will follow the ﬁrst two iterations

of PHA for r = 2 here.

Iteration 0:

Step 0. Begin with a multiplier vector of

= 0 , and let x

=(x

)=(0,10)

and let x

=(x

)=(10, 0)

. The initial value of ˆx

=(5,5)

Step 1.Wewishtosolve:

min(1/2)[y

+ y

+(x

−5)

+(x

−5)

+(x

−5)

+(x

−5)

]

s. t. x

+ x

≤ 10 ,

+ x

≤ 10 ,

+ 3x

−y

≥ 25 ,

+ 2x

−y

≥ 25 ,

≥ 0 .

(8.26)

This problem splits into separate subproblems for x

, x

, y

and x

, y

, as mentioned earlier. For x

, x

, y

feasible in (8.26), the K-K-T

conditions are that there exist

≥ 0,

≥ 0 such that

2(x

−5)+

−

≥ 0 ,

2(x

−5)+

−3

≥ 0 ,

5.8 Methods Based on the Stochastic Program Lagrangian 261

≥ 0 ,

(2(x

−5)+

−

= 0 ,

(2(x

−5)+

−3

= 0 ,

(2y

= 0 ,

+ x

−10)

= 0 ,

+ 3x

−y

−25)

= 0 , (8.27)

which has a solution of (x

)=(10/3,20/3,5/3) and (

(20/3,10/3) . Similar conditions exist for the second subproblem, which has a so-

lution (x

)=(5, 5,0) .Wethenlet ( ˆx

, ˆx

)=(4

) for i = 1, 2.

Step 2. The new multiplier is

)

= 2((10/3−25/6), (20/3−

35/6),(5 −25/6),(5 −35/6))

=(−5/3, 5/3, 5/3,−5/3)

Iteration 2:

Step 1. The ﬁrst subproblem is now

min y

−(5/3)(x

−25/6)+(5/3)(x

−35/6)

+(x

−25/6)

+(x

−35/6)

s. t. x

+ x

≤ 10 ,

+ 3x

−y

≥ 25 ,

≥ 0 ,

(8.28)

which again has an optimal solution, (x

)=(10/3,20/3,5/3) . Curiously,

we also have the second subproblem solution of (x

)=(10/3,20/3,0) .In

this case, ( ˆx

, ˆx

)=(10/3,20/3) for i = 1,2.

Step 2. Because the subproblems returned the same solution,

. We continue

because the x values changed, even though we took no multiplier step.

The full iteration values are given in Table 1. Notice how the method achieves

convergence in the x values before the

values have converged. Also, notice how

the convergence appears to be geometric. This type of performance appears to be

typical of PHA. It should be noted again, however, that the iterations are quite simple

and that little overhead is required.

Exercises

1. Show that the basic dual ascent method converges to an optimal solution under

the conditions given.

262 5 Two-Stage Recourse Problems

Table 1 PHA iterations for Example 3.

k ˆx

ˆx

= −

0 5.0 5.0 0.0 0.0 3.33 6.67 5.0 5.0

1 4.17 5.83 -1.67 1.67 3.33 6.67 3.33 6.67

2 3.33 6.67 -1.67 1.67 3.06 6.94 2.50 7.50

3 2.78 7.22 -1.11 1.11 2.78 7.22 2.41 7.59

4 2.59 7.41 -0.74 0.74 2.65 7.35 2.41 7.59

5 2.53 7.47 -0.49 0.49 2.59 7.41 2.43 7.57

6 2.50 7.50 -0.33 0.33 2.56 7.44 2.45 7.55

7 2.50 7.50 -0.22 0.22 2.54 7.46 2.46 7.54

8 2.50 7.50 -0.15 0.15 2.53 7.48 2.48 7.52

9 2.50 7.50 -0.10 0.10 2.52 7.48 2.48 7.52

10 2.50 7.50 -0.07 0.07 2.51 7.49 2.49 7.51

11 2.50 7.50 -0.04 0.04 2.51 7.49 2.49 7.51

12 2.50 7.50 -0.03 0.03 2.50 7.50 2.50 7.50

2. Show that (8.4) can be reduced to (8.8)when g

(y(

)=T(

)x+Wy(

)−

) , f

(y(

)=q(

)

)y(

) ,and D is positive

deﬁnite.

3. Show that V as deﬁned in (8.13) is a maximal monotone operator.

4. Prove the contraction property in (8.14).

5. Use (8.14) to obtain (8.23).

6. Apply the dual ascent method and the augmented Lagrangian method with prob-

lem (8.9) to the example in (8.25). Start with zero multipliers (

= 0or1,

and positive penalty r . Show that each obtains an optimal solution in at most

one iteration.

5.9 Additional Methods and Complexity Results

In the previous sections, we considered cutting plane methods and Lagrangian meth-

ods for problems with discrete random variables and simple recourse-based tech-

niques for problems with continuous random variables. Other nonlinear program-

ming procedures can also be applied to stochastic programs, although these other

procedures have not received as much attention in stochastic programming prob-

lems. A notable exception is No¨el and Smeers’ [1987] multistage combined inner

linearization and augmented Lagrangian procedure, which we will describe in more

detail in the next chapter.

A difﬁculty with discrete random variables is that

or Q generally loses

differentiability. In this case, derivative-based methods cannot apply. As we saw,

5.9 Additional Methods and Complexity Results 263

the L -shaped method and other cutting plane approaches are a standard approach

that requires only subgradient information. We also saw that augmented Lagrangian

techniques can smooth nondifferentiable functions.

Explicit nondifferentiable methods include the nonmonotonic reduced subgradi-

ent procedure considered by Ermoliev [1983]. Another possibility is to use bundles

of subgradients as in Lemar´echal [1978] and Kiwiel [1983]. Results by Plambeck

et al. [1996], for example, show good performance for bundle methods in practical

stochastic programs.

Nonsmooth generalizations of the Frank-Wolfe procedure are also possible.

These and other options are described in detail in Demyanov and Vasiliev [1981].

With general continuous random variables or with large numbers of discrete ran-

dom vector realizations, direct nonlinear programming procedures generally break

down because of difﬁculties in evaluating function and derivative values. In these

cases, one must rely on approximation. These approximations either take the form

of bounds on the actual function values or are in some sense statistical estimates of

the actual function values. We present these approaches in Chapters 8 to 10.

While models with discrete random variables inherit the complexity results

of their deterministic equivalent forms with possible improvements due to prob-

lem structure as shown for interior point methods in Section 5.5, general dis-

tributions can present difﬁculties even in the two-stage case. For the common

mean-variance objective, for example, the two-stage stochastic program is NP-hard

(Ahmed [2006]). While exact solutions to general stochastic programs are difﬁcult

in general, bounds may be obtained efﬁciently using the methods in Chapter 8 and

other approaches that can achieve a priori bounds on error in special cases. For

example, Dye, Stougie, and Tomasgard [2003] consider a problem of a central re-

source serving facilities with random demands; Gupta, et al. [2007] provide bounds

on the related stochastic Steiner tree problem to connect a source node to termi-

nal nodes that are randomly revealed in the second period; Ravi and Sinha [2006]

provide results for the stochastic shortest path version with generalizations to other

combinatorial problems; and Flaxman, Frieze, and Krivelevich [2005] give a so-

lution for a two-stage stochastic spanning tree problem, where instead of random

demand, uncertainty is in the cost of edges which can be purchased for known costs

in the ﬁrst period and for random costs in the second period. Swamy and Shmoys

[2006] provide a survey of these and other approaches including sampling methods

which are discussed in Chapter 9.