Computational solution of capacity planning models under uncertainty

Parallel Computing 26 (2000) 511±538

www.elsevier.com/locate/parco

Computational solution of capacity planning models under uncertainty S.A. MirHassani, C. Lucas, G. Mitra *, E. Messina, C.A. Poojari Department of Mathematics and Statistics, Brunel University, Uxbridge, Middelesex UB8 3PH, UK Received 20 July 1998; received in revised form 2 October 1999

Abstract The traditional supply chain network planning problem is stated as a multi-period resource allocation model involving 0±1 discrete strategic decision variables. The MIP structure of this problem makes it fairly intractable for practical applications, which involve multiple products, factories, warehouses and distribution centres (DCs). The same problem formulated and studied under uncertainty makes it even more intractable. In this paper we consider two related modelling approaches and solution techniques addressing this issue. The ®rst involves scenario analysis of solutions to ``wait and see'' models and the second involves a two-stage integer stochastic programming (ISP) representation and solution of the same problem. We show how the results from the former can be used in the solution of the latter model. We also give some computational results based on serial and parallel implementations of the algorithms. Ó 2000 Elsevier Science B.V. All rights reserved. Keywords: Benders decomposition; Scenario analysis; Strategic planning; Two-stage stochastic programming; Parallel algorithm

1. Introduction The planning and utilisation of production capacity is a major strategic decision in manufacturing. The strategic decisions of a company are concerned with acquiring the resources needed to survive and prosper over the long term. Such decisions have to be made in the face of uncertainty, which arise in the future realisations of demand, price and technology data. When evaluating alternative decisions the

*

Corresponding author. Tel.: +44-1895-203304; fax: +44-1895-203303. E-mail address: [email protected] (G. Mitra).

0167-8191/00/$ - see front matter Ó 2000 Elsevier Science B.V. All rights reserved. PII: S 0 1 6 7 - 8 1 9 1 ( 9 9 ) 0 0 1 1 8 - 0

512

S.A. MirHassani et al. / Parallel Computing 26 (2000) 511±538

potential bene®ts of new resources in conjunction with the existing ones must be considered. All these challenging aspects of manufacturing capacity planning have long attracted the attention of economists as well as operational researchers and management scientists. In recent years the supply chain planning problem has been gaining importance due to the increasing competitiveness introduced by the market globalisation (for a review see [53]). Firms are obliged to maintain high customer service levels while at the same time they are forced to reduce cost and maintain pro®t margins. Strategic supply chain decisions are particularly important in the context of the corporate plan. Set against the requirement of reducing the inventory costs, a high quality customer service can be guaranteed only by adopting ¯exible con®gurations, which allow a prompt reaction to highly volatile customer demand. Traditionally this class of planning problems has been formulated as large mixed integer programs, which are dicult to solve (e.g., [8,10,25]). The solutions naturally depend on the data used to instantiate the model. Many aspects of this data are quite volatile and subject to wide variations, typically, national and global economic factors [29] which a€ect selling price and demand for goods as well as reliability (down time) of the production units. Some uncertainty is due to economic factors while others arise as a consequence of the introduction of new technologies. Thus, very large investment decisions are made on the basis of relatively unreliable forecasted data. This has prompted researchers to investigate ways of building planning models which can take into account wide variations in the data values which instantiate the model parameters. By taking such an approach, a robust optimisation model can be developed whose solution can provide decisions, which are near optimal for many realisations of model parameters. Stochastic programming with recourse models are ideally suited for analysing resource acquisition planning problems from two perspectives. They combine deterministic mathematical programming models for allocating resources optimally with decision analysis models that provide hedging strategies in an uncertain environment. Stochastic programs (SPs) which explicitly consider randomness are much more versatile and provide more information than a deterministic representation of the strategic plans and tactical operations. SP models have been successfully applied to strategic planning problems. In particular, two-stage stochastic programming models have been developed in a number of applications, e.g., electricity utility planning [5], goods distribution [9], capacity planning [18,44,54], communication network planning [21], transportation planning and vehicle routing [20,34]. The success of this modelling technique is due to its ability to combine ``here and now'' decisions with ``wait and see'' decisions. Indeed, in a co-ordinated supply chain network, decisions are made over time at an increasing level of detail. Strategic decisions such as production capacity levels are taken over a long planning horizon while operational decisions such as contingent production and distribution plans can be taken over a shorter time horizon. Furthermore, at the time when strategic decisions are made there is no detailed knowledge of many decision parameters such as future product demands and prices. At the time of making operational decisions, however, more information is available and this can be used to optimise the production plans with respect to contingent information. A detailed description of how

S.A. MirHassani et al. / Parallel Computing 26 (2000) 511±538

513

the hierarchical structure of planning systems can be exploited for building multistage stochastic decision models can be found in [7,13,14,51]. This hierarchical structure is well exploited by decomposition solution methods. Among decomposition methods, Benders decomposition [52], scenario aggregation, progressive hedging algorithms [46] and stochastic decomposition [28] have been successfully applied to solving a number of stochastic problems. An important class of stochastic models is the two-stage stochastic linear program with recourse. In the two-stage stochastic linear program with recourse, a ®rst-stage decision is made which minimises the cost of the ®rst-stage decisions (variables) as well as the expected cost of the second-stage decisions (variables) as a consequence of making the ®rst-stage decisions. This expected second-stage cost is referred to as the recourse function. Consider the deterministic optimisation problem P1. Let x 2 Rm1 and y 2 Rm2 represent the ®rst and second-stage decision variables and let k1 and k2 denote the number of ®rst- and second-stage constraints. P1 : min Z ˆ cx ‡ fy …1:1† s:t: Ax ˆ b; ÿ Bx ‡ Dy ˆ d;

…1:2† …1:3†

x; y P 0; where c; x 2 Rm1 ;

f ; y 2 Rm2 ;

b 2 Rk 1 ;

d 2 Rk 2 ;

A 2 Rk1 m1 ;

B 2 Rk2 m1 ;

D 2 Rk2 m2 : In the deterministic case the vectors and matrices c, f, A, b, B, D and d are known with certainty to the planner. In the stochastic case, the values of the second-stage parameters fs , Bs , Ds and ds are alternative realisations of the uncertain values of the parameters f, B, D and d of the model. Uncertainty may be included in the model in many ways. One such way is to represent the uncertain quantities as random variables, so that the ®rst-stage decisions are made while the second-stage parameters are known only by their probability distributions. Actual outcomes will be known later, at the time of taking the second-stage decisions. In this investigation we focus only on the case where the probability distributions are discrete. We denote a particular outcome of these random variables by a scenario s with corresponding probability ps , s 2 Sc, where Sc is the set of indices representing all possible outcomes and S ˆ jScj. The general formulation of a two-stage stochastic linear program with recourse is as follows: …1:4† P2 : min Z ˆ cx ‡ E…f s y s † s:t: Ax ˆ b; …1:5† …1:6† ÿ Bs x ‡ Ds y s ˆ d s 8s; x; y s P 0; s 2 Sc: Extensions of deterministic models to incorporate uncertainty in other ways have also been developed. Wagner and Berman [54] discuss di€erent ways of taking uncertainty into account.

514

S.A. MirHassani et al. / Parallel Computing 26 (2000) 511±538

In the strategic planning problem that we consider in this paper, uncertainty occurs in the right-hand side vector d and it is an example of a ``®xed recourse'' problem where the technological matrix D remains unchanged (deterministic) for each scenario s. Problem (P2) can be interpreted as follows. A decision-maker must select the activity levels for x before the actual values of the random vectors are known. After the values of the random vectors are known, the activity levels of y are chosen. Hence, the ®rst-stage activities are chosen so as to minimise the total cost of the ®rst stage and the expected cost of the second stage. The main drawback of SP formulations of capacity planning problems is that they are computationally demanding, in particular when discrete (integer) restrictions are imposed on some of the ®rst-stage and/or second-stage variables. When the second-stage variables y are restricted to be integer the SP problem loses some desirable properties such as convexity and continuity of the second-stage cost function, see [48,49]. Creating and investigating large-scale stochastic programming problems require developing tools for modelling and solving SP models [40,31]. Solution methods that are successful for stochastic linear programming cannot be simply adapted to the stochastic problem with integer second-stage variables. A number of approximation algorithms have been proposed in the literature, e.g., [13,14,51]. Exact optimisation algorithms have also been developed and in the case when the stochastic parameters follow a discrete distribution Lenstra et al. [33] proposed an algorithm based on dynamic programming. For a survey of these solution techniques see [50]. In a number of applications, such as ours, only the ®rst-stage variables are restricted to be integer. These problems often arise in capital investment, capacity planning, transportation, agricultural planning and plant location applications where decisions correspond to an investment policy and consequently can be represented as 0±1 variables. The second-stage problem stated in continuous variables retains its convexity property of the recourse function and ``conventional'' methods can still be applied that are successful in dealing with the convex second-stage costs. In our particular case only the right-hand side values are stochastic in the second stage. Special purpose algorithms exist for solving LP problems with multiple right-hand side values, see [23]. Although this provides a theoretical framework for solving the problem, the computational complexity still remains the major diculty. Decomposition methods are the most widely applied solution techniques for solving these problems because they can simultaneously exploit the stochastic programming and mixed integer programming model structures. Wollmer [55] proposed a cutting plane algorithm for the solution of the integer stochastic programming (ISP) problem where second-stage variables are continuous and the random variables are discrete. Multi-stage programs have been addressed by Louveaux [37] and Bienstock and Shapiro [5] that exploit their special structure to formulate and solve an equivalent two-stage problem with integer variables in the ®rst stage only. However, these approaches are well suited only when the resulting master problem (®rst-stage model) is tractable and when there are a limited number of realisations of the random variables, i.e., only a small manageable number of scenarios are considered. Moreover, Laporte and Louveaux [35] provide an algorithm for problems

S.A. MirHassani et al. / Parallel Computing 26 (2000) 511±538

515

with ®rst-stage binary decision variables where the random vectors and the secondstage decisions variables may be discrete or continuous. In particular, decomposition methods lend themselves to creating parallel computational algorithms. This aspect has been substantially explored by a number of workers, e.g., [6,15,45]. Parallel versions of progressive hedging algorithms have also been developed in [2,3,27,41,42]. In respect of useful applicable models of reasonable size not much computational work has been undertaken. One important aspect of parallelising solution algorithms for these problems is that it increases the reach to the user community involved in practical decisionmaking under uncertainty. The structure of these problems can be exploited for designing specialised solution procedures based on decomposition methods (e.g., [5,35]). Decomposition methods when applied to ISP problems are often supported and guided by heuristic procedures in order to speed up the convergence. When a large number of scenarios are to be considered or when the subproblems are dicult to solve, then heuristic procedures may be applied to reduce the number of subproblems solved at each iteration of the decomposition algorithm [5]. When the computational diculty arises in the solution of the master problem, as in our case, a heuristic must be designed in order to speed up the convergence and reduce the number of passes through the two stages. In order to design e€ective heuristics an accurate analysis of the speci®c problem must be carried out. The ®rst step to gain insight into the problem is the solution of the deterministic ``wait and see'' model for each scenario. By analysing the optimal solutions obtained against all available scenarios one can understand the behaviour of the model under di€erent realisations of the stochastic parameters. The complete probability distribution of the optimal solution value to the stochastic programming problem could also be computed. A discussion on the importance of the distribution problem and its related solution method in the case of 0±1 problems can be found in [56]. We use the compact notation a/b/c adopted by Laporte and Louveaux [35] in order to summarise some of the computational work that has taken place in recent years, where: a b c

describes the ®rst-stage variables and is either C (continuous), B (binary) or M (mixed); describes the second-stage variables and is either C, B or M; describes the probability distribution of the unknown parameters and is either Dd (discrete distribution ) or Cd (continuous distribution).

Refs.

ISP

Scale

Parallel

[5] [6] [30] [15] [19] [18]

M/M/Dd C/C/Dd C/C/Dd C/C/Dd C/C/Dd C/C/Dd

Small Large Large Large Large Large

No Yes Yes Yes Yes No

516

[22] [28] [33] [35] [36] [41] [43] [47] [55]

S.A. MirHassani et al. / Parallel Computing 26 (2000) 511±538

C/C/Dd C/C/Cd B/C/Dd B/C/Cd M/B/Dd C/C/Dd M/C/Dd C/C/Dd B/C/Dd

Large Large Small Small Large Large Large Large Small

No No No No No Yes Yes Yes No

The structure of this paper is as follows. In Section 2, we describe a capacity planning model whose ®rst stage contains integer decision variables with the second stage containing stochastic continuous variables. An analysis of the ``wait and see'' models is presented in Section 3. The results of this analysis are later incorporated into the stochastic representation of the model. BenderÕs decomposition method which has gained considerable acceptance as a solution technique for two-stage SPs is introduced in Section 4. The customisation of this method to solve two-stage integer SPs is also discussed in this section. In Section 5, two alternative decompositions of the problem are introduced and the computational results of solving the prototype problem on a serial computer are presented. In Section 6, the parallel adaptation of the solution algorithms and the corresponding computational results are set out. Finally, our conclusions are reported in Section 7. 2. The network model structure The model under consideration is formulated so as to represent the entire manufacturing chain, from the acquisition of raw material to the delivery of ®nal products. The four major stages of the manufacturing process are production and packing at plants (sites), and transportation through (DCs), to customer zones (CZs). (see Fig. 1). The ®rst-stage strategic decisions are concerned with the opening and closing of sites and DCs and setting their capacity levels, in terms of how many and of what

Fig. 1.

S.A. MirHassani et al. / Parallel Computing 26 (2000) 511±538

517

type of line technology to have operating both in sites and DCs. All these strategic decisions are represented by discrete variables. The second-stage operational decisions, such as production quantities, packing quantities and transportation amounts are represented as continuous variables. These account for the vast majority of variables in the model. The model itself is made up of logical and operational constraints. The logical constraints ensure that lines can only exist on sites that are open and ensure physical capacity limitations at sites are obeyed. The operational constraints representing the vast bulk of the model maintain material balances and measure the shortage of demand. A full algebraic description of the model can be found in [38]. The only uncertain component in this problem is represented by the demand for products. The problem owner (decision-maker) has considered a reasonable number (100) of different demand scenarios. The resulting deterministic model even with one scenario is very large and contains many integer variables. Because of the size of the problem and the number of integer variables it is impossible to solve the deterministic equivalent model to optimality or even verify optimality of the decomposition method. It is also very dif®cult to solve a single occurrence of one scenario of demand. In order to investigate the performance of di€erent solution techniques we consider a simpli®ed model hereafter referred to as the prototype model, which maintains the same model structure. The statistics of the fully detailed planning model and the reduced prototype model are given in Table 1.

Table 1 Network dimensions The The The The The The The The

number of sites, I types of packing line technology, YC types of production line technology, YR number of distribution centres, J types of DC line technology types, YD number of customer zones, H number of products, P number of time periods, T

Model statistics Logical constraints: sites, DCs opening and closing, limit on number of sites, DCs, and lines Continuous constraints: production, packing, ordering, transportation, balance, demand, and also production and packing capacities Discrete decision variables: sites, DCs, production lines, packing lines, DC lines Continuous variables: production, packing, ordering, transportation, and shortage quantities Non-zeros Density Scenarios

Prototype model

Planning model

2 1 1 3 1 4 2 6

8 4 2 15 2 30 13 6

268

1460

155

4520

197

1195

218

44 021

2265 1.3% 50

118 884 0.044% 100

518

S.A. MirHassani et al. / Parallel Computing 26 (2000) 511±538

If we ignore the last row of Table 1, the rest of the table summarises a single scenario deterministic mixed integer linear programming model presented as a prototype model in an aggregated form or as a detailed planning model of much larger dimension. These two models may be considered as a single scenario twostage SP in which discrete variables specify the ®rst-stage decisions and the logical constraints the ®rst-stage model restrictions. The deterministic mixed integer linear program is initially used to validate logical relationships and usefulness of the model. Subsequently repeated solutions of the model for all scenario realisations provide some insight into the robustness of the planning decisions and the distributions of the optimum solution values in respect of the stochastic demand pattern. The computation of an exact optimum solution for a two-stage integer SP of the size of the detailed planning model is pretty much intractable. Hence, we have developed heuristic procedures which take the results of the ``wait and see'' model and use these results to provide acceptable near-optimal solutions to the two-stage integer SPs. 3. Investigation of ``wait and see'' models The stochastic version of the prototype model set out in Table 1 is still very large when we consider the number of scenarios that can be generated. In order to analyse a large-scale stochastic programming problem it is opportune to solve the ``wait and see'' problems for di€erent scenarios. This analysis may be used to achieve a reduction in the features of the deterministic (one scenario) model making the stochastic model more tractable. When performing this scenario analysis we ®rst ®nd the discrete solution xs of the MIP model for one particular scenario s, that will be referred to as ``con®guration''. Hence, a con®guration is a solution to the ®rst-stage discrete variables, x. Since we do not prove optimality of this solution, it is necessary to carry out evaluation of con®gurations against all the scenarios in order to assess their performance. The mean and the variance of the LP optimal values of each con®guration for all the scenarios can be taken as performance indices. The expected value of the objective function obtained for each con®guration represents a lower bound on the stochastic optimal objective function, e.g., [54]. 3.1. A strategy for scenario analysis Let PMIP (s) denote a MIP problem for the data instance of scenario s and xs the (sub)optimal integer solution. Therefore, for each s ˆ 1; . . . ; S we de®ne: PMIP …s† :

min

Z ˆ cx ‡ fy

…3:1†

s:t: Ax ˆ b; ÿ Bx ‡ Dy ˆ d s ;

…3:2† …3:3†

x 2 f0; 1g;

y P 0:

S.A. MirHassani et al. / Parallel Computing 26 (2000) 511±538

519

Fig. 2.

Let PLPR (s) denote the LP relaxation of the problem PMIP (s). Let PLP …s; j† denote the LP problem when the ®rst-stage variables, x, are ®xed to xs under scenario j. Hence, for s ˆ 1; . . . ; S; j ˆ 1; . . . ; S we de®ne: PLP …s; j† :

min

Z ˆ cxs ‡ fy

s:t: Dy ˆ d j ‡ Bxs ; y P 0:

…3:4† …3:5†

It is easy to see that the problems set out in (3.1)±(3.5) lead to S, PMIP problems and S  …S ÿ 1†; PLP problems. The major steps for scenario analysis are called phases. These phases are outlined below: Phase 1. Finding an integer solution for each scenario (computing a con®guration). Phase 2. Evaluating the con®gurations. Phase 3. Carrying out solution matching and integer ®xing. The overall structure of the algorithm is given in Pseudo-code 1. The optimum value of the LP model solved for scenario j and con®guration s represents the e€ect that con®guration s would have if scenario j occurred. The results of the

520

S.A. MirHassani et al. / Parallel Computing 26 (2000) 511±538

computations carried out in phase 1 and phase 2 create a vector of S con®gurations and a pay-o€ table of S  S objective function values, respectively, as shown in Fig. 2. The pay-o€ table illustrates for each scenario j the consequence of implementing con®guration s. This is an important analysis since many of the con®gurations found in phase 1 are suboptimal. Thus, in this analysis we determine the Min PLP …s; j† 8j. This identi®es the con®guration s that leads to the best solution for scenario j. If PMIP (j) could be solved to optimality then Min PLP …s; j† 8j would equal PLP …j; j†. Then for every con®guration s, s ˆ 1; . . . ; S, we determine the number, ns , of times this con®guration leads to the best solution for scenario j; j ˆ 1; . . . ; S.

S.A. MirHassani et al. / Parallel Computing 26 (2000) 511±538

521

In general, we use scenario analysis to: · assess the performance of a con®guration taking into consideration the extreme scenarios (low demand and high demand); · identify di€erent scenarios which lead to similar solutions; · identify good con®gurations which are then used to carry out variable ®xing; · compute the expected value and the variance of the pro®t for each con®guration. 3.2. Solution matching and integer ®xing The scenario analysis leads to a subset of ``good'' con®gurations, Sgc, which have a positive value of nj (j ˆ 1; . . . ; S) (con®guration with nj ˆ 0 is no longer considered). We next analyse the solution values of this subset of con®gurations establishing which variables, if any, have the same solution value for all these con®gurations. These variables are then ®xed to this matching value and are removed from the problem in further investigations. Thus, given the subset of ``good'' con®gurations xs , with s 2 Sgc ˆ fs1 ; s2 ; s3 ; . . .g  Sc and a solution vector of the ®rst-stage variables, v, for the realisation of scenario s1 2 Sgc, we analyse individual components of the solution values and if xks are the same for all s1 2 Sgc, we then ®x xks to this value. The steps of the solution matching and integer ®xing algorithms are set out in Fig. 3, where xks denotes the kth component of the con®guration for scenario s. 3.3. Summary of computational results Computational tests were carried out in serial on one processor of a 12 processor PARSYTEC CC12 computer, in which each distributed unit is a power PC-604,

Fig. 3.

522

S.A. MirHassani et al. / Parallel Computing 26 (2000) 511±538

Fig. 4.

Fig. 5.

133 MHz and has 64 MB of local memory. We also used the optimisation solver FortMP [16] developed by our group. The three phases of scenario analysis have been carried out on both the prototype model and the real planning model with 50 and 100 scenarios, respectively. In the ®rst phase we identify a set of (sub)optimal con®gurations. In the second phase we evaluate the performance of the con®gurations against scenarios. Other analyses can also be carried out on the results. For example, computing the average value of a con®guration for all scenarios. The results of analysing the planning model are shown in Fig. 4, which shows a histogram of nj that is the frequency distribution of the best con®guration. The expected best and worst optimum values for all con®gurations are shown in Fig. 5. The processing times for all these phases of scenario analysis are displayed in Table 2. Table 2 Models

Phases

Prototype

Phase Phase Phase Phase

Planning

1 2 1 2

Time (s) (IP) (LP) and phase 3 (IP) (LP) and phase3

710.21 169.7 92549.6 115540.7

S.A. MirHassani et al. / Parallel Computing 26 (2000) 511±538

523

4. Benders decomposition 4.1. De®nition of master and subproblems Our main interest is to solve the two-stage ISP problem. We naturally start with Benders decomposition since this has become a widely accepted method of solving twostage stochastic LP problems. This decomposition method was originally introduced [4] to solve mixed integer programming problems. Since that time it has been extended to process stochastic programming problems and we ®nd it a natural approach to solve the integer stochastic problem at hand. We reformulate the problem P2 as P3: P3 :

min Z ˆ cx ‡ Q…x† s:t: Ax ˆ b; x P 0;

…4:1† …4:2†

where Q(x) is called the ``recourse function'' and is de®ned as Q…x† ˆ Ef min … fy s † j Dy s ˆ d s ‡ Bx; y s P 0; s 2 Scg:

…4:3†

Benders decomposition method is applied to solve problem P3. This involves solving a ®rst-stage master problem in tandem with a sequence of second-stage subproblems P5. The master problem is set out in (P4): P4 :

min Z ˆ cx ‡ h s:t: Ax ˆ b; x P 0; X ps Ps
…d s ‡ Bx†; hP s2Sc

Ps
…d s ‡ Bx† 6 0;

s 2 Sc;

…4:4† …4:5† …4:6† …4:7†

where the restrictions (4.6) or (4.7) are introduced only after solving the second-stage problems. There are s ˆ 1; . . . ; S second-stage subproblems P5 which are de®ned using the optimum solution values x ˆ x of the ®rst-stage problem. In (4.6) Ps is the vector of dual variables of a typical second-stage problem (P5) and ps is the corresponding probability for scenario s: P5 :

min Z s ˆ fy s s:t: Dy s ˆ d s ‡ Bx;

y s P 0:

…4:8† …4:9†

After solving the second-stage subproblems, if subproblem (P5) has no feasible solution (i.e., its dual is unbounded), then a feasibility cut (4.7) is added to the master problem. If all subproblems are feasible or optimal (dual bounded), then an optimality cut (4.6) is added. A recourse function that is ®nite (P5 is feasible) for any choice of x is said to be complete. Therefore, after each pass either an optimality cut (4.6) or a feasibility cut (4.7) is introduced into the master problem. Let ubl be the upper bound, and lbl be the lower bound in the lth pass. These are calculated as follows:   …4:10† lbl ˆ C x ‡ h ; n    o ubl ˆ min ublÿ1 ; C x ‡ Z x ; …4:11†

524

S.A. MirHassani et al. / Parallel Computing 26 (2000) 511±538

P where Z…x† ˆ s2Sc ps Z s and …x; h † is the solution of the master problem in the lth pass. When the relative gap, i.e., …ubl ÿ lbl †=…jlbl j ‡ 1†, between the upper bound and the lower bound falls within a given tolerance then the problem is said to be solved with sucient accuracy. 4.2. Algorithmic framework In the formulation of the planning problem we have introduced appropriate shortage variables, as a result all second-stage subproblems are feasible for nonnegative demand quantities. The Pseudo-code 2 describes a simpli®ed Benders decomposition algorithm that does not take into account feasibility cuts and only generates optimality cuts. Throughout our work in solving subproblems in serial, we utilise the fact that the only stochastic parameters are presented by demand levels, which are the right-hand side values. Thus, after solving an initial subproblem the subsequent subproblems are ``warm started'' by specifying their starting basis to be the optimal basis of the previous subproblem.

S.A. MirHassani et al. / Parallel Computing 26 (2000) 511±538

525

5. Alternative decompositions 5.1. The standard representation We ®rst apply Benders decomposition as described in Section 4 to a standard representation of the two-stage problem. We observe that the algorithm starts by ®nding a poor initial solution to the master problem and converges slowly to an optimal solution. At each pass in this process an extra row is introduced to the master problem which usually makes it more dicult to solve. When dealing with integer variables normally much more computational time is required to obtain an optimal or at least a good suboptimal solution. In many cases it is impractical to search for an optimal solution in a reasonable time. Therefore, it is important to ®nd ways to decrease the total number of passes as well as the computational time for the master problem and subproblems. We next apply Benders decomposition algorithm to our problem while ®xing the appropriate strategic decision variables to the values found from the scenario analysis and integer matching procedure discussed in Section 3.2. The results are reported in Table 3. This table illustrates that by using the results of the solution analysis and integer matching leads to an e€ective decrease in the computational time while not a€ecting the solution quality. 5.2. The expanded representation In order to reduce the number of total passes that are needed to obtain an optimal solution, it is desirable to have a good initial solution to the master problem. Usually a ``good'' solution is one that is not in con¯ict with the subproblems. This means the master problem should be more representative of the subproblems. By including a scenario, say k, in the master problem we obtain the following expanded representation (the problem of choosing the most e€ective scenario k is discussed later): X ps fy s …5:1† min Z ˆ CX ‡ s6ˆk

s:t:

 AX ˆ b; ÿ BX ‡ Dy s ˆ d s ; X ; y s P 0; s 2 Sc0 ˆ Sc n fkg;

…5:2† …5:3†

Table 3 Standard master and subproblems Models

LP

IP

Time (s)

No. passes

Equ. deterministic Standard Standard and ®x

7.0849E + 7 ± ±

7.0865E + 7 7.0865E + 7 7.0865E + 7

3708.24 5910.7 273.35

± 146 74

526

S.A. MirHassani et al. / Parallel Computing 26 (2000) 511±538

where  X ˆ x " A ˆ

y

 k T

B ˆ ‰ B

;

0 Š;

  C ˆ c; pk f ;

#

A

0

ÿB

D

" b ˆ

b dk

# and

8s 2 Sc0

Therefore, the master problem is: min

Z ˆ CX

…5:4†

s:t:

 AX ˆ b;

…5:5†

X P 0 and a typical subproblem for s 2 Sc n fkg is: min s:t:

Z ˆ fy s

…5:6†

Dy s ˆ d s ‡ BX ;

…5:7†

y s P 0;

s 2 Sc0 ˆ Sc n fkg:

This formulation changes the master problem from an IP to a larger MIP problem while reducing the number of subproblems by one, without any change in the structure of the LP subproblem. By including a scenario we change the master problem to re¯ect a weighted realisation of a typical scenario. This makes the master problem more re¯ective of each subproblem and the initial master strategic decisions now take into account information about future observations of the environment. We call this the expanded master problem. The selection of which scenario to use in the expanded master problem is an outcome of our investigations on all the scenarios. Fig. 6 shows the structure of the standard master and subproblems versus the expanded master problem and subproblems.

Fig. 6.

S.A. MirHassani et al. / Parallel Computing 26 (2000) 511±538

527

The expanded model has been solved using the same algorithmic framework as given in Section 4.2. The results in Table 4 show that the expanded decomposition is more ecient in terms of the number of passes and computational time compared with the standard representation. Fig. 7 shows the movement of the objective value under di€erent decomposition methods after each pass until optimality is reached. It is easy to see that by adding a scenario to the standard master problem considerably improves the solution quality of the master problem in the ®rst few passes. The e€ect of choosing di€erent scenarios for inclusion in the master problem is seen in Table 5. All problems are solved with the same approach. It is easy to see that the number of passes when including scenario 3 in the master problem is the most e€ective in terms of computational time and passes.

Table 4 Expanded master and subproblems Models

LP

IP

Time (s)

No. passes

Equ. deterministic Standard Expanded

7.0849E + 7 ± ±

7.0865E + 7 7.0865E + 7 7.0865E + 7

3708.24 5910.7 151.9

± 146 20

Fig. 7.

Table 5 Scenario K

Obj. value (IP)

No. passes

Time (s)

Demand

3 14 29

7.0865E + 7 7.0865E + 7 7.0865E + 7

20 156 226

151.9 919.6 1625.91

342718.15 235984.55 122877.80

528

S.A. MirHassani et al. / Parallel Computing 26 (2000) 511±538

Since the only uncertain component in this problem is represented by the demand for products, the total demand for each scenario can be considered as a parameter to explain di€erent passes for di€erent scenario choices in the master problem. Scenario 3 represents the biggest demand in our model. The results shown in Table 5 highlight how the total demand e€ects the total number of passes needed to reach the solution. 5.3. The weighted representation After considering the results in Table 5 we conclude that a strategy with high demand, which calls for large amounts of capacity, provides a good initial solution to the master problem in the Benders decomposition solution process. We therefore investigated another decomposition in which the direction of optimisation for the master problem is reversed. Thus, in the ®rst-stage the master problem is changed to a maximisation problem. This is achieved by introducing a strictly negative weight a to multiply the objective cx, of the master problem. Therefore, we rewrite problem (P2) as follows: min s:t:

Z ˆ …ÿacx† ‡ f…1 ‡ a†cx ‡ E… fy s †g

…5:8†

Ax ˆ b;

…5:9†

ÿ Bx ‡ Dy s ˆ d s x; y s P 0;

s 2 Sc;

8s;

…5:10†

0 < a 6 1:

The objective of the weighted master problem becomes min…ÿacx ‡ h†, or equivalently max…acx ÿ h†. Thus, the master problem is stated as a maximisation problem: max s:t:

Z 0 ˆ acx ÿ h

…5:11†

Ax ˆ b;

…5:12†

hP

X

ps Ps
…d s ‡ Bx†;

…5:13†

s2Sc

Ps
…d s ‡ Bx† 6 0 x P 0;

8s 2 Sc;

…5:14†

S.A. MirHassani et al. / Parallel Computing 26 (2000) 511±538

529

where 0 < a 6 1. A typical subproblem for s 2 Sc is: min s:t:

Z ˆ …1 ‡ a†cx ‡ fy s

…5:15†

Dy s ˆ d s ‡ Bx;

…5:16†

y s P 0: Fig. 8 shows the structure for the standard master and subproblems versus the weighted master and subproblems. We have solved the weighted model using the algorithmic framework outlined in Section 4.1. The results of solving the weighted master problem and subproblems are reported in Table 6. In the weighted master approach the number of passes to

Fig. 8.

Table 6 Weighted master and subproblems Models

LP

IP

Time (s)

No. passes

Equ. deterministic Standard Weighted

7.0849E + 7 ± ±

7.0865E + 7 7.0865E + 7 7.0865E + 7

3708.24 5910.7 21.99

± 146 5

Table 7 Weight (a)

IP

Time (s)

No. passes

0.2 0.4 0.8 1.0

7.0865E + 7 7.0865E + 7 7.0865E + 7 7.0865E + 7

21.99 85.3 32.1 52.7

5 11 5 6

530

S.A. MirHassani et al. / Parallel Computing 26 (2000) 511±538

achieve the optimum solution as well as the computational time decrease in comparison with the standard method. This is because the master problem is much easier to solve at each pass. The number of passes are reduced since ``maximisation'' ®nds a good initial solution for the master. By varying the weight over a range of values we have investigated the e€ect on the number of passes required by the Benders decomposition algorithm. The results set out in Table 7 illustrate that for the given prototype model, in¯uence of the weight parameter is limited and the number of passes range between 5 and 11. 5.4. Computational results The results set out in Table 8 illustrate the importance of ®nding a good decomposition representation. Clearly, the expanded method compared to the standard method leads to a reduction in the number of passes but there is more computational time in each pass. The weighted approach is better on both counts. It has a fairly fast computational time in each pass while its convergence to the solution is quicker than the expanded approach. This is particularly important in solving the real life application where it is crucial to have as small a master problem as possible. In Fig. 9 we have described in a ``schematic form'' the nature and behaviour of the three alternative decompositions. The standard decomposition has an initial master solution, which starts with a minimum amount of strategic resources (sites, DCs, and lines serving the supply chain) and with proposing cuts arrives at the optimum resource allocation. The Table 8 Models

LP

IP

Time (s)

No. passes

Equ. deterministic Standard Expanded Weighted

7.0849E + 7 ± ± ±

7.0865E + 7 7.0865E + 7 7.0865E + 7 7.0865E + 7

3708.24 5910.7 151.9 21.99

± 146 20 5

Fig. 9.

S.A. MirHassani et al. / Parallel Computing 26 (2000) 511±538

531

expanded model due to the inclusion of a representative scenario in the master gives an initial strategic resource solution, which is near to the optimum allocation. The weighted model on the other hand starts with a maximum (``greedy'') allocation of strategic resources but then rapidly reduces these resources through the decomposition process. 6. Investigation of parallel algorithms 6.1. Scope of parallel computing Multi-processor computational platforms and parallel algorithms have been deployed to address a number of intractable problems, which would otherwise require substantial numerical processing time even on a fast serial computer. Essentially the motivation for using parallel processing is to: I. obtain a good speed up of execution time [1]; II. scale up the data instances of a given problem class [26]. Both these aspects have the common objective of achieving solutions of acceptable instances of practical applications in a reasonable time. Dantzig [12] has long been a proponent of applying parallel algorithms to SP problems and over the last 10 years a number of researchers [15,30,32,42,45,47] have presented parallelization of SP solution Algorithms. Our approach is to not only parallelise the solution algorithm, but also to achieve performance gains on most readily available hardware and software platforms. Thus, we have used both special purpose parallel computers PARSYTEC CC-12 (discussed in Section 3.3) and have also adapted: (a) a collection of high speed PCs; (b) these are connected by a 100 MB hub; (c) and PVM [24] is used as the message passing interface to achieve parallel processing. Each of the ®ve PC stations is a 333 MHz Pentium II, 128 MB RAM, VIGLEN computer and the operating system is windows NT4. The cluster is con®gured as a single program multiple data (SPMD) distributed processing environment. In the algorithms presented in this section we have used a master/slave tree topology parallel programming model. Another aspect of our approach is to present these solution algorithms essentially as black boxes which may be used as a scalable optimisation server. Using client/server [11] software tools our aim is to scale up from a prototype model to realistic applications with very little additional e€orts by the end users and modelling experts. In scenario analysis, of the two main solver tasks, phase 1 requires an independent collection of similar models to be solved; therefore each scenario is solved on one slave processor. Whereas in phase 2 each slave processor evaluates one con®guration against all scenarios. The two-stage SP also displays substantial subproblem independence where the solver is con®gured for a two-stage Benders decomposition algorithm.

532

S.A. MirHassani et al. / Parallel Computing 26 (2000) 511±538

6.2. Scenario analysis in parallel The natural independence of problems in scenario analysis suggests that much of the computational e€ort of the scenario analysis algorithm can run in parallel. Let N be the number of available processors, which are made up of 1 master processor and N ÿ 1 slave processors.

S.A. MirHassani et al. / Parallel Computing 26 (2000) 511±538

533

Table 9 Models

Phase

Prototype Prototype Planning Planning

Phase Phase Phase Phase

1 2 1 2

(IP) (LP) and phase3 (IP) (LP) and phase 3

Serial time (s)

Parallel time (s)

710.21 169.7 92549.6 115540.7

317.60 38.81 19435.4 23570.3

As before (see Section 3) we apply a three-phase algorithm which is shown in Pseudo-code 3. It is obvious from the computational times shown in Section 3.1 that it is not worth confusing the algorithm logic by parallelizing phase 3. In Table 9 we have shown the processing times for the prototype as well as the full planning problem. Note: (a) In phases 1 and 2 it was observed that a static distribution of the scenariobased problem led to uneven loads and poor speed up. Hence, an ``on the ¯y'' dynamic load distribution was implemented which gave better load balancing and near linear speed up. See Table 9. (b) Starting with the prototype model, scenario analysis can be scaled up for the planning model.

6.3. Benders decomposition in a parallel framework To solve the ISP problem in a parallel environment we modify our previous algorithm. We consider the following parallel implementation for a master/slave strategy for solving the problem. An outcome of our investigations into solving subproblems in serial, was that by using the previous basis as the initial basis for the next subproblem led to the creation of deep cuts in the master problem. When we carry out this task in a parallel environment we solve unique groups of subproblems on each processor. Initially the scenarios were chosen in a random fashion and the subproblems were spread over the master and slave processors. The resulting cuts were less deep and hence less e€ective than those created in the serial approach. This approach led to many more passes in the Benders decomposition when compared with the number of passes required in the serial algorithm. As a result of the analysis of the pay-o€ table from phase 2 we are able to identify an important scenario (high demand scenario) and solve this on the master processor in serial. In order to restore the bene®cial performance of the serial approach, we broadcast the basis of the ®rst subproblem from the master processor to all slave processors. Hence, each slave processor warm starts from the same basis and works on its own unique group of subproblems.

534

S.A. MirHassani et al. / Parallel Computing 26 (2000) 511±538

6.4. Computational results A parallel implementation of the Benders decomposition is shown in Pseudo-code 4. The master problem is solved on the master processor using a serial algorithm. All experiments were carried out on a PARCYTEC CC12 computer (see Section 3.3) Table 10 Models Standard Expanded Weighted

IP, Obj. 7.0865E + 7 7.0865E + 7 7.0865E + 7

Serial

Parallel

Time (s)

No. passes

Time (s)

No. passes

5910.7 151.9 21.99

146 20 5

5119.4 112.6 7.18

173 18 5

S.A. MirHassani et al. / Parallel Computing 26 (2000) 511±538

535

using the prototype model. Table 10 illustrates the time taken in solving the twostage problem in serial and parallel.

7. Conclusions Supply chain planning in the face of future uncertainties is an important as well as a computationally challenging decision problem. The problem can be well analysed applying scenario analysis [17] as well as two-stage SP. However, practical instances of such problems for a multinational with a global sourcing requirement (D. Gregg, private communication) lead to instances of models which are intractable as SPs with todayÕs solver technologies. In this paper we have shown how scenario analysis can be easily scaled up for practical problems and be used to provide what-if answers to sensitivities which an end user (decision-maker) may wish to investigate. Although two-stage SP can be parallelised with considerable simplicity it may not converge as well as the serial Benders approach. We have shown how a parallel Benders algorithm can be modi®ed to achieve a convergence rate comparable with the serial Benders approach. For the given problem we have investigated alternative decomposition representations which take into consideration the nature of the master problem. We have shown that by reversing the sense of the master problem from minimisation to maximisation we can speed up the solution process. Finally, we have identi®ed a number of ways in which we might be able to improve the computational performance of the two-stage SP solution algorithm. We believe that the introduction of a Lagrangean method for solving the master problem may improve the solution time since the master is a dicult discrete optimisation problem. We can also apply parallel branch and bound [39] to solve the master problem which is currently solved in serial. These are topics of our future investigations.

Acknowledgements The research reported in this paper was supported by EPSRC ROPA grant awarded number GR/K92443. Also Mr. MirHassani has been supported by the Iranian government. The work of Dr. A. Nagar in carrying out computational tests is gratefully acknowledged.

References [1] G.M. Amdahl, Validity of the single processor approach to achieving large-scale computing capabilities, AFIPS Conference Proceedings 30 (1967) 483±485.

536

S.A. MirHassani et al. / Parallel Computing 26 (2000) 511±538

[2] K.A. Ariyawansa, D.D. Hudson, Performance of a benchmark parallel implementation of the Van Slyke and Wets algorithm for two-stage stochastic programs on the sequent/balance, Concurrency Practice and Experience 3 (1991) 109±128. [3] N.J. Berland, K.K. Haugen, Mixing stochastic dynamic programming and scenario aggregation, Annals of Operation Research 64 (1996) 1±19. [4] J.F. Benders, Partitioning procedures for solving mixed variables programming problems, Number Math 4 (1962) 238±252. [5] D. Bienstock, J.F. Shapiro, Optimising resource acquisition decisions by stochastic programming, Management Science 34 (2) (1984). [6] J.R. Birge, C.J. Donohue, D.F. Holmes, O.G. Svintsiski, A parallel implementation of the nested decomposition algorithm for multi-stage stochastic linear programs, Mathematical Programming 75 (2) (1996) 327±352. [7] G.R. Bitran, D. Tirupati, Hierarchical production planning, in: S.C. Graves, A.H.G. Rinnoy Kan, P.H. Zipkin (Eds.), Handbook on Operations Research and Management Science, vol. 4, NorthHolland, Amsterdam, 1993, pp. 523±568. [8] G.G. Brown, G.W. Graves, M.D. Honczarenko, Design and operation of a multi-commodity production/distribution system using primal goal decomposition, Management Science 33 (11) (1987) 1469±1480. [9] R.K.M. Cheung, W.B. Powell, Models and algorithms for distribution problems with uncertain demands, Transportation science 30 (1) (1996) 43±59. [10] M.A. Cohen, H.L. Lee, Strategic analysis of integrated production±distribution systems: models and methods, Operations Research 36 (2) (1988) 216±228. [11] J. Cruise, G. Mitra, D. Hague, Client/Server integration tools for mathematical programming systems, Presented at Informs National Meeting, San Diego, 4±7 May, 1997. [12] G.B. Dantzig, Planning under uncertainty using parallel computing, in: R.R. Meyer, S.A. Zenios (Eds.), Parallel Optimization on Novel Computer Architectures (special issue), Annals of Operations Research 14 (1988). [13] M.A.H. Dempster, M.L. Fisher, L. Jansen, B.J. Lageweg, J.K. Lenstra, A.H.G. Rinnoy Kan, Analytical evaluation of hierarchical planning systems, Operations Research 29 (4) (1981) 707±716. [14] M.A.H. Dempster, M.L. Fisher, L. Jansen, B.J. Lageweg, J.K. Lenstra, A.H.G. Rinnoy Kan, Analysis of heuristics for stochastic programming: results for hierarchical scheduling problems, Operations Research (1984). [15] M.A.H. Dempster, R. Thompson, EVPI-based importance sampling solution procedures for multistage stochastic linear programmes on parallel MIMD architecture, Annals of OR 90 (1999) 161±184. [16] E.F.D Ellison, M. Hajian, R. Levkovitz, I. Maros, G. Mitra, D. Sayers, FortMP manual, Brunel University, West London and NAG Ltd, 1996. [17] G.D. Eppen, R.K. Martin, L. Schrage, A scenario approach to capacity planning, Operations Research 37 (4) (1989) 517±525. [18] L.F. Escudero, P.V. Kamesam, A.J. King, R.J.B. Wets, Production planning via scenario modelling, Annals of Operations Research 43 (1993) 311±335. [19] L.F. Escudero, J.L. de la Fuente, C. Garcia, F.J. Prieto, Hydropower generation management under uncertainty via scenario analysis and parallel computation, IEEE Transaction on Power Systems 11 (2) (1996) 683±689. [20] A. Fefergruen, P. Zipkin, A combined vehicle routing and inventory allocation problem, Operations Research 32 (5) (1984) 1019±1037. [21] F. Fantauzzi, A. Gaivornski, E. Messina, Decomposition methods for network optimisation problems in the presence of uncertainty, in: P.M. Pardalos, D.W. Hearn, W. Hager (Eds), Lecture Notes on Economics and Mathematical Systems: Network Optimization, vol. 450, Gainesville, FL, USA, 1996, pp. 234±248. [22] H.I. Gassmann, MSLIP: a computer code for the multi-stage stochastic linear programming problem, Mathematical Programming 47 (1990) 407±423. [23] H.I. Gassman, S.W. Wallace, Solving linear programs with multiple right-hand sides: pricing and ordering schemes, Annals of Operations Research 64 (1996) 237±259.

S.A. MirHassani et al. / Parallel Computing 26 (2000) 511±538

537

[24] Al. Geist, A. Begulin, J. Dongarra, W. Jiang, R. Mancheck, V. Sunderam, PVM 3 UserÕs guide and reference manual, 1994. [25] A.M. Geo€rion, G.W. Graves, Multi-commodity distribution system design by Benders decomposition, Management Science 20 (5) (1974) 822±844. [26] J.L. Gustafson, Re-evaluating AmdahlÕs law, Communications of the Association for Computimg Machinery 31 (5) (1988) 532±533. [27] T. Helgason, S.W. Wallace, Approximate scenario solutions in the progressive hedging algorithm, Annals of Operations Research 31 (1991) 425±445 1991. [28] J.L. Higle, S. Sen, Stochastic Decomposition: A Statistical Method for Large-scale Stochastic Linear Programming, Kluwer Academic Publisher, Dordrecht, 1996. [29] A. Huchzermeier, M.A. Cohen, Valuing operational ¯exibility under exchange rate risk, Operations Research 44 (1) (1996). [30] G. Infanger, Planning under Uncertainty: solving large scale stochastic linear program, Boyd and Fraser, Danvers, MA, 1994. [31] P. Kall, J. Mayer, SLP-IOR: an interactive model management system for stochastic linear programs, Mathematical Programming 75 (1996) 221±240. [32] A. King, SP/OSL version 1.0. Stochastic programming interface library: userÕs guide, Research Report RC, 19757 (87525) 9/26/94 Mathematics, T.J Watson Research Center, Yorktown Heights, NY. [33] B.G. Lageweg, J.K. Lenstra, A.H.G. Rinnooy Kan, L. Stougie, Stochastic integer programming by dynamic programming, Statistica Neerlandica 39 (2) (1985). [34] G. Laporte, F.V. Louveaux, H. Mercure, The vehicle routing problem with stochastic travel times, Transportation Science 26 (1992) 161±170. [35] G. Laporte, F.V. Louveaux, The integer L-shape method for stochastic integer programs with complete recourse, Operation Research Letters 13 (1993) 133±142. [36] C.M.A. Leopoldino, M.V.F. Pereira, L.M.V. Pinto, C.C. Riberio, A constraint generation scheme to probabilistic linear problems with an application to power system expansion, Annals of operations research 50 (1994) 367±385. [37] F.V. Louveaux, Multi-stage stochastic programs with block-separable recourse, Mathematical Programming Studies 28 (1986) 48±62. [38] C. Lucas, E. Messina, G. Mitra, Risk and return analysis of a multi-period strategic planning problem, in: L. Thomas, A. Christers (Eds.), Stochastic Modelling in Innovative Manufacturing, Springer, Berlin, 1996. [39] G. Mitra, I. Hai, M.T. Hajian, A distributed processing algorithm for solving integer programs using a cluster of workstations, Parallel Computing 23 (1997) 733±753. [40] E. Messina, G. Mitra, Modelling and analysis of multi-stage stochastic programming problems: a software environment, European Journal of Operational Research 101 (1997) 343±359. [41] J.M. Mulvey, A. Ruszczynski, A new scenario decomposition method for large-scale stochastic optimization, Operations Research 43 (3) (1995) 477±490. [42] J.M. Mulvey, H. Vladimirou, Evaluation of a parallel heading algorithm for stochastic network programming, in: R. Sharda, B. Golden, E. Wasil, O. Balci, W. Stewart (Eds.), Impacts of recent computer advances on operations research, Elsevier, Amsterdam, 1989. [43] J.M. Mulvey, H. Vladimirou, Parallel and distributed computing for stochastic network programming, Technical Report SOR-90-11, Department of Civil Engineering and Operations Research, Princeton University, Princeton, NJ, 1990. [44] E.M. Modiano, Derived demand and capacity planning under uncertainty, Operations Research 35 (2) (1987) 185±197. [45] S.S. Nielsen, S.A. Zenios, Solving multi-stage stochastic network programs on massively parallel computers, Mathematical Programming 73 (3) (1996) 227±250. [46] R.T. Rockafellar, R.J.B. Wets, Scenario and policy aggregation in optimization under uncertainty, Mathematics of Operations Research 16 (1991) 119±147. [47] A. Ruszczynski, Parallel decomposition of multi-stage stochastic programming problems, Mathematical Programming 58A (2) (1993) 201±228. [48] R. Schultz, Continuity and stability in two-stage stochastic integer programming, in: K. Marti (Ed.), Lecture Notes in Economics and Mathematical Systems, vol. 379, Springer, Berlin, 1992, pp. 81±92.

538

S.A. MirHassani et al. / Parallel Computing 26 (2000) 511±538

[49] R. Schultz, Continuity properties of expectation functions in stochastic integer programming, Mathematics of Operations Research 18 (1993) 578±589. [50] R.L. Schultz Stougie, M.H. Van der Vlerk, Two-stage stochastic integer programming: a survey, Statistica Neederlandica 50 (3) (1996) 404±416. [51] R. Stougie, Design and analysis of algorithms for stochastic integer programming, CWI Tract, vol. 37, Centrum voor Wiskunde en Informatica, Amsterdam, 1987. [52] van Slyke, R.J.-B. Wets, L-shaped linear programs with application to optimal control and stochastic programming, SIAM Journal on Applied Mathematics 17 (1969) 638±663. [53] D.J. Thomas, P.M. Grin, Coordinated supply chain management, European Journal of Operational Research 94 (1996) 1±15. [54] J.M. Wagner, O. Berman, Models for planning capacity expansion of convenience stores under uncertain demand and the value of information, Annals of Operations Research 59 (1995) 19±44. [55] R.D. Wollmer, Two-stage linear programming under uncertainty with 0±1 integer ®rst-stage variables, Mathematical programming 19 (1980) 279±288. [56] M.A. Zimmermann, M.A. Pollatschek, The probability distribution function of the optimum of a linear program with randomly distributed coecients of the objective function and the right-hand side, Operations Research 23 (1) (1975) 137±149.