Strategic Capacity Allocation under Uncertainty by Using a Two-Stage

2812

Ind. Eng. Chem. Res. 2010, 49, 2812–2821

Strategic Capacity Allocation under Uncertainty by Using a Two-Stage Stochastic Decomposition Algorithm with Incumbent Solutions Sergio Frausto-Hernandez,*,† Vicente Rico-Ramirez,‡ and Ignacio E. Grossmann§ Departamento de Ingenierıá Quı´mica y Bioquı´mica, Instituto Tecnolo´gico de Aguascalientes, AVe. Adolfo Lo´pez Mateos No. 1801 Ote., Fracc. Bona Gens, Aguascalientes, Ags., C. P. 20256, Me´xico, Departamento de Ingenierıá Quı´mica, Instituto Tecnolo´gico de Celaya, AVe. Tecnolo´gico y Garcıá Cubas S/N, Celaya, Gto., C.P.38010,Me´xico,andDepartmentofChemicalEngineering,CarnegieMellonUniVersity,Pittsburgh,PennsylVania15213

This paper addresses the optimization of strategic capacity allocation problems under uncertainty. Uncertain parameters include the availabilities and the demands of chemicals. The optimization problem is formulated as a two-stage mixed-integer stochastic program (first-stage integer) that maximizes the expected net present value of the project over the time horizon. The net present value includes the discounted investment costs, operating costs, and revenues from sales in the markets. The decisions about the selection of processes and about the plant capacities are defined as the first stage decision variables; given such first stage decisions, recourse is then provided by calculating the operating costs in the second stage. The paper presents the twostage model and the solution framework. A stochastic decomposition algorithm (SD) based on incumbent solution is used as the solution algorithm. The use of incumbent solutions as convergence criterion allows us to solve the problem in reasonable computational time. Two capacity allocation problems are used as case studies: the first one is a small size problem and is used to compare the stochastic optimal solution against other solution methods; the second one is a larger problem, with an increased number of variables and stochastic parameters. The results show that the stochastic decomposition algorithm with incumbent solutions is stable and robust and greatly reduces the number of samples required for convergence. 1. Introduction Supply chain optimization is attracting a great deal of attention from industry and academia.1 While the applications are aimed mostly at commodities, there are also several applications in specialties, pharmaceuticals,2 and in the food industry.3 There is not a systematic way of defining the scope of a supply chain problem. Chopra and Meindl proposed categorizing the supply chain problem as design problems, planning problems, and operational or scheduling problems.4 The characterization of the problems depends on the time frame over which the decisions made apply. In this work we deal with the strategic capacity allocation problem. Supply chain strategy decisions are typically made for the long-term and are expensive to change once they have been established. Consequently, when companies make these decisions, they must take into account uncertainty over the future. Strategic decisions include location and capacities of production and warehouse facilities, products to be manufactured or stored, modes of transportation, and the type of information system to be shared among the chain. Major challenges in this area include the development of models for strategic and tactical planning for process networks5 that often requires the solution of large-scale multiperiod optimization problems. Furthermore, these models must be eventually integrated with scheduling models. The incorporation of uncertainty in strategic models through stochastic optimization still remains a great challenge due to the very large computational requirements that are needed.6 However this is clearly an area that might be ripe for significant * To whom correspondence should be addressed. E-mail: serfraher@ yahoo.com.mx. † Instituto Tecnolo´gico de Aguascalientes. ‡ Instituto Tecnolo´gico de Celaya. § Carnegie Mellon University.

progress, including ways for better quantifying financial risk (e.g., see the work of Barbaro and Bagajewicz7). Uncertainty occurs in both the objective function and the constraints. The objective function includes uncertainty in future costs, prices, and actions. The process model constraints can include both endogenous (e.g., yield coefficients) and exogenous (e.g., future product demands) stochastic parameters.8 Failure to account for the uncertainty of key parameters (e.g., future supplies and demands) in decisions problems can lead to sub optimal or infeasible decisions (e.g., see the work of Birge and Louveaux9). The most widely used models in stochastic optimization to solve strategic and planning problems under uncertainty can be credited to earlier works such as those of Dantzig10 and Beale.11 However, it was until the decade of the 1990s that the first research works on the solution of stochastic linear programs were reported (see, for example, the works of Kall and Wallace12 and Haneveld and Vander Vlerk13). In particular, the stochastic decomposition (SD) algorithm developed by Higle and Sen14 and the L-shaped method described in the work of Birge and Louveaux9 are the two main algorithms for solving two-stage stochastic linear programming problems. The SD algorithm is suitable for problems involving uncertain parameters represented by continuous distribution functions. The L-shaped method is suited for discrete distribution functions of the uncertain parameters. Among recent developments in chemical engineering, PonceOrtega et al.15 developed a computational framework for solving two-stage stochastic linear programs with recourse. The Hammersley sequence sampling technique was used within the basic stochastic decomposition algorithm.14 The authors used small size planning problems to test the implementation and analyze the computational effort needed for convergence. Recently, Rico-Ramı´rez et al.16 extended the implementation of the basic SD algorithm to address two-stage stochastic models with integer first-stage decisions; a formulation representing

10.1021/ie901021b  2010 American Chemical Society Published on Web 02/16/2010

Ind. Eng. Chem. Res., Vol. 49, No. 6, 2010

2813

Figure 1. Flows within a supply chain.

security in municipal water networks against a chemical attack was used as a case study. That problem involves a large number of continuous variables, binary variables, stochastic parameters, and scenarios. It was thus shown that the basic SD algorithm is suitable for solving large scale-problems. However, the number of iterations (or samples) needed for convergence and the CPU time greatly increases with the size of the problem. Motivated by these results, the goals of this work are as follows: (i) First, we propose to reformulate the strategic capacity allocation model of Iyer and Grossmann17 as a stochastic program. Thus, as our first contribution, the problem is now represented as a two-stage mixed-integer stochastic problem that maximizes the expected net present value of the project over the horizon. (ii) To solve the two-stage stochastic problem, a computational implementation of the SD algorithm is proposed. We have implemented the SD algorithm with a convergence criterion based on incumbent solutions14 in order to reduce the large computational requirements of the basic SD algorithm. A stochastic two-stage supply chain model with the uncertainties of future supplies and demands characterized by continuous probability distribution functions is used as a case study. (iii) Finally, the last goal of the paper is to compare the performance of the improved implementation of SD algorithm with other solution techniques suitable for discrete distributions of uncertain parameters. Section 2 describes the problem statement and the two-stage stochastic reformulation. The solution approach is reviewed in section 3. The results for two illustrative examples are presented in section 4. Finally, conclusions are summarized in section 5.

2. Strategic Capacity Allocation in Supply Chain Management A supply chain can be thought of as a network of trading partners converting raw materials into finished product for customers (see Figure 1). Each entity provides some activity necessary for this transformation. These entities include, but are not limited to, customers, retailers, wholesalers, distributors, feedstock suppliers, and manufacturing facilities. Interactions among them can take the form of materials, information, or monetary flow. Typically, there is a forward flow of goods and backward flow of information throughout the supply chain. However, modern and efficient supply chains require a forward flow of information as well. Supply chain management is then the design, planning, execution, control, and monitoring of supply chain activities with the objective of creating net value; this task involves the management of flows among the stages in a supply chain to maximize profitability. The objective of any supply chain is thus to minimize the total supply chain cost while meeting customer demand within a target level of customer service. The major costs areas affected by supply chain decisions are the following: inventory holding costs, inventory shortage costs, and transportation costs.18 A simple production network can be used to illustrate the strategic supply chain problem.8 This network is shown in Figure 2. Material balance nodes are represented as circles, and processing units, as rectangles. Raw materials and intermediates are purchased on the markets subject to prevailing prices, contracts, and availability. Products are sold on the market based on demands, prevailing prices, contracts, and production capacity. As described above, the main goal here is to address strategic capacity allocation problems in the context of supply chain management. The following subsections provide the problem

2814


Figure 2. Simple production network.

statement and the reformulation of the strategic problem as a two-stage linear program with integer decisions in the first stage. 2.1. Problem Statement. Consider the problem proposed by Iyer and Grossmann.17 The problem consists in a production superstructure of NP chemical processes and NC chemicals (including raw materials, intermediates, and products). There are NK streams in the network, and the material balances for each process is expressed linearly based on the production rate of the main product. The demands and prices of chemicals may vary for each time period over the time horizon, but in this model, it is assumed that they are known with precision. The objective function to be maximized is the expected net present value (NPV) of the project over the horizon consisting of NT time periods. The NPV includes the discounted investment costs, operating costs, and revenues from sales in NM markets. 2.2. Two-Stage Stochastic Model for Strategic Capacity Allocation. As our first contribution, we propose the reformulation of the model of Iyer and Grossmann17 as a two-stage stochastic model. Uncertainties in the availability and the demand of the products are explicitly incorporated into the formulation. An assumed continuous probabilistic distribution function is assigned to each of these parameters which makes the problem a stochastic program. Also, we reformulate the model constraints into two-stages, following the concept of recourse. We have conserved the notation of the original reference.17 Hence, the indices, parameters, and variables defined in the model are as follows: Indices i ) process (i ) 1, ..., NP) j ) chemical (j ) 1, ..., NC) t ) time period (t ) 1, ..., NT) k ) stream in network (k ) 1, ..., NK) l ) market (1, ..., NM) Parameters I(j) ) index set of streams of chemical j produced in the complex O(j) ) index set of streams of chemical j consumed in the complex Li ) index set of streams corresponding to inputs/outputs of process i mi ) stream corresponding to the main product of process i (mi ∈ L i) µik ) material balance coefficient for process i and stream k NEXP(i) ) maximum number of allowable expansions of process (i) CI(t) ) capital investment limit in period t

Rit ) variable investment cost coefficient adjusted to the NPV of process i ($/unit capacity) R j it ) variable investment cost coefficient not adjusted to the NPV of process i ($/unit capacity) βit ) fixed investment cost coefficient adjusted to the NPV of process i ($) j it )fixed investment cost coefficient not adjusted to the NPV β of process i ($) δit ) unit operating cost of process i ($/unit production amount) γjlt ) sales price of chemical j in market l ($/unit sold) Γjlt ) purchase cost of chemical j in market l ($/unit purchased) θjlt ) purchase cost of shortfall chemical j for market l ($/unit shortfall purchased) u1jlt ) availability of chemical j in market l ($/sold) u2jlt ) demand of chemical j in market l ($/purchased) Variables yit ) decision variable for expansion of process i in period t, equal to 1 if there is an expansion, 0 otherwise Qit ) capacity of the process i in period t QEit ) capacity expansion of process i in period t Fjlt ) shortfall of chemical j in market l Pjlt ) amount of product j purchased in period t in market l Sjlt ) amount of product j sold in period t in market l Wkt ) amount of flow of stream k in period t Superscripts ()L )lower bound ()U )upper bound The first stage includes the decisions about the selection of processes and production capacities through the binary variables yit. Also, the total capacity Qit and the expansion capacity QEit are continuous variables in the first-stage. The objective function in eq 1 represents the maximization (through minimizing the negative value) of the NPV, which includes the investment costs of the process and their expansions (Rit and βit) as well as the expected revenues by the sum of operating costs and cost of purchase of chemicals, EuR(Qit, ujlt). Equation 2 bounds the capacity expansion in each time period. Equation 3 defines the total capacity in each time period t, and eq 4 limits the number of expansions NEXP(i) for process i over the time horizon. Equation 5 limits the amount of capital investment CIt in each j it are coefficients of costs time period t. The parameters R j it and β not adjusted to the NPV. First-stage NT

min -

NP

∑ ∑ (R QE it

t)1 i)1

it

+ βityit) + EuR(Qit, ujlt)

(1)


subject to yitQELit e QEit e yitQEUit

time periods are lumped into the second stage (e.g., see the work of Dantzig19 and Bienstock and Shapiro20).

i ) 1, ..., NP;t ) 1, ..., NT

(2) Qit ) Qi,t-1 + QEit

i ) 1, ..., NP;t ) 1, ..., NT

(3)

NT

∑y

it

i ∈ I' ⊆ {1, 2, ..., NT}

e NEXP(i)

(4)

t)1

NP

∑ (R¯ QE it

it

+ β¯ ityit) e CIt

t ) 1, ..., NT

(5)

i)1

QEit, Qit g 0 yit ∈ {0, 1}

(6)

In the second-stage, the recourse function R(Qit, ujlt) includes operating costs, sales costs, purchase costs, and shortfall costs by demand not satisfied (eq 7). Pjlt and Sjlt are the amounts of purchases and sales of chemicals j in the market l; Wkt is the amount of flow of the stream k. All these variables are assumed as second-stage decision variables. Equation 8 imposes production limits based on existing capacity, and eq 9 is a mass balance for each chemical in the complex. Equation 10 represents the mass balance equations of chemicals in each process i, and eqs 11 and 12 incorporate uncertainties in the availabilities (u1jlt) and the demands of the chemicals (u2jlt), respectively. In the model, we have assumed a continuous probabilistic distribution function to define the availabilities and demands of the chemicals. Notice that, in eq 12, total sales available for chemical j might not be able to satisfy the demand for the market l; for this reason, we consider a shortfall production Fjlt, so that the second-stage problem always leads to a feasible solution (complete recourse condition). This last term is going to be penalized in the objective function by a cost θjlt for a potential loss of sales or for a breaking a contract. Second-stage NT NC NM

R(Qit, ujlt) ) min - (

NT

∑ ∑ ∑γ

jltSjlt

t)1 j)1 l)1 NT NC NM

∑ ∑ ∑Γ

-

NP

∑ ∑δ W it

t)1 i)1 NT NC NM

mit

-

∑ ∑ ∑θ

(7)

i ) 1, ..., NP;t ) 1, ..., NT

(8)

jltPjlt

-

t)1 j)1 i)1

jltFjlt)

t)1 j)1 i)1

subject to Wmit e Qit NM

∑ j)1

Pjlt +

∑

NM

Wkt )

k∈I(j)

∑S

jlt

j)1

+

∑

Wkt

j)

k∈O(j)

1, ..., NC;t ) 1, ..., NT

Wkt ) µikWmitk ∈ Li /mi ;

Sjlt + Fjlt ) u2jlt

(9)

i ) 1, ..., NP;t ) 1, ..., NT (10)

l ) 1, ..., NM;t ) 1, ..., NT

Pjlt e u1jlt

l ) 1, ..., NM;t ) 1, ..., NT

Wkt, Pjlt, Sjlt, Fjlt g 0

2815

(11) (12) (13)

We should note that another alternative for this formulation is the development of a multistage mixed-integer linear program with recourse (MSMILPwR). However, in this paper, all future

3. Solution Approach This section describes an extension to the basic SD method that helps to improve the computational performance of the algorithm and stabilize the approximations and iterates generated by the method. The extensions discussed in this paper rely on both the introduction of a sequence of incumbent solutions (extracted from the sequence of master program solutions) and its use as stopping criterion. 3.1. Incorporation of Incumbent Solutions. An incumbent solution is a vector of first-stage decision variables that are found to have the lowest objective function value than all the other vectors found at some iteration in the search process of the basic SD algorithm.12 The first value assigned to the incumbent solution is to the one that solves the problem with the mean values of the stochastic parameters. Within the basic SD algorithm, any solution of the master program may become the incumbent solution. If the estimated objective value associated with the solution is sufficiently low, the solution becomes the incumbent. Because the objective function approximations used in SD are based on statistical estimates, incumbent solutions are not guaranteed to have the lowest objective value among those encountered during the iterative process. Nevertheless, as iterations proceed, the estimated objectives associated with an incumbent solution become more accurate and, therefore, asymptotic optimality of the incumbents can be ensured. Kall and Wallace12 illustrated a procedure to identify and verify if an incumbent should be changed during the iterative procedure (see Figure 3). Following the incumbent, there will be an index iν that shows in which iteration the current jxν was found, and we shall use the function φν(x) to define the actual estimate of φ(x) ) cTx + Q(x). In Figure 3, we have ν ) 3. When the algorithm is at iteration ν, the approximation of φ(x) is originally given by φν-1(x) ) φ2(x) (see lines AB and BC in Figure 3a). Also, the incumbent solution is jxν-1 ) jx2 and the iterate xν ) x3. The position of jx2 is somewhat arbitrary, since we have no information from the previous iteration. Hence, φν-1(xν) - φν-1(xjν-1) ) φ2(x3) - φ2(xj2) (segment AD in Figure 3a) is the current approximation of how much one should gain by making xν ) x3 our new incumbent. However, xν ) x3 might be in an area where φν-1(x) ) φ2(x) is a bad approximation of φ(x) because it was obtained by using two cuts (AB and BC). The function φν(x) ) φ3(x), on the other hand, was developed around xν ) x3 (see Figure 3a) and should therefore be a good estimate around that area, the same way that φν-1(x) ) φ2(x) is a good estimate around jxν-1 ) jx2, because they were obtained using three and two cuts, respectively (see lines EF, FG, and GH in Figure 3b and lines AB and BC in Figure 3a). Hence φν(xν) - φ(xjν-1) ) φ3(x3) - φ3(xj2) is a measure of how much one would actually gain if we the incumbent solution is changed (see segment EI in Figure 3b). If φν(xν) - φν(xjν-1) < r[φν-1(xν) - φν-1(xjν-1)] ) φ3(x3) φ3(xj2) < r[φ2(x3) - φ2(xj2)] ) EI < rAI (14) one would gain at least a fraction r ∈ (0,1) of what is it expected and could define xν as the new incumbent solution; that is, jxν ) xν and iν ) ν. If eq 14 is not satisfied, the change is not satisfactory either, and the incumbent solution should not change jxν ) jxν-1 and iν ) iν-1. Once the incumbent has been updated,

2816


Figure 3. Graphical illustration of the incumbent solution update.

a new master problem (see section 3.3) is solved to obtain xν+1 and the iteration is repeated. 3.2. Termination Criterion. Termination criteria based on the convergence of the incumbent objective values are easily envisioned. Let γν represent some statistical measure of the incumbent objective values obtained during the first ν iterations. Furthermore, suppose that γν satisfies the condition |fjν(xjν-1) - γν| f 0. If the sequence {fjν(xjν-1)} converges, γν is a statistical mean, γν )

ν 2σˆ jf ν 1 jf k(xjk) + V k)1 √ν

∑

(15)

where the term (2σˆ jfν)/(ν) is a penalization error.21 A termination criterion based on the convergence of the incumbent objective value would test whether the following condition is satisfied.

|

|

γν - γν-1 eε γν

(b) Determine the coefficients of the optimality cut: V

max πT(hν - TVxν) s.t. πTW e q to find the values of the vector of the second stage dual variables (Lagrange multipliers of primal restrictions) π, πνν, and set Vν ) Vν-1 ∪ πνν.

1 (πV)Th V k)1 k k

βνν )

1 (πV)TT V k)1 k k

V

∑

where πνk is the solution to the problem (for all k|k * ν) max πT(hk - Tkxν) s.t. π ∈ Vν Determine the cut that will be added to the first stage problem, as well as the linear approximation of the problem f. θ g Rνν + βννx

(16)

where ε > 0 is an acceptable tolerance level. Otherwise, the number of the samples ν must be increased. 3.3. SD Algorithm with Incumbent Solutions. The steps of the SD algorithm that uses incumbent solutions as termination criterion to obtain a faster convergence is similar to the one proposed by Elshorbagy et al.22 and can be summarized as follows:13 • Step 0 Set the iteration counter ν ) 0, V0 ) {L}, u0 ) E[u], the initial values of the first decision variables, x1, are assumed as given. jx0 ) x1, iν ) 0, f0 ) 0, jf0 ) 0, F0 ) 0. Let ∆1 (a sufficiently large number) and r ∈ (0,1) be given. • Step 1 Set ν ) ν + 1 and sample to generate an observation uν independent of any previous observation. • Step 2 Determine the coefficients of a piecewise linear approximation to the recourse function, Q(x) (optimality cut): (a) Solve the dual program of the second stage problem:

∑

Rνν )

f ν ) Rνν + (c + βνν)xν (c) Update the coefficients of previous cuts: RVk )

V - 1 V-1 Rk , V

k ) 1, ..., V - 1

βVk )

V - 1 V-1 βk , V

k ) 1, ..., V - 1

(d) If ν ≡ iν + 1, then re-evaluate the iνth incumbent cut as follows: Solve the dual program of the second stage problem max π¯ T(hν - TVjxν) s.t. π¯ TW e q to find the values of the vector of incumbent variables of the second stage dual problem π j, π j νν, and make Vν ) j νν. Vν-1 ∪ π Determine the coefficients of the optimality cut of the incumbent solution V

Riνν )

∑

1 (π¯ V)Th V k)1 k k


2817

V

βiνν )

∑

1 (π¯ V)TT V k)1 k k

Determine the cut of the incumbent solution that will be added to the first stage problem, as well as, the linear approximation of the objective jfν. θ g Riνν + βiννx jf ν ) Riν + (c + βiν )xjν ν ν (e) Test the incumbent. If f ν - jfν < rην then jxν ) xν, iν ) ν, ∆ν ) 2∆ν-1. γν ) jf ν Otherwise, jxν ) jxν-1, iν ) iν-1, ∆ν ) ∆ν-1/2. Figure 4. Solution SD incumbent algorithm. ν

γν )

2σˆ jf ν-iν

∑

1 jf (xjk) + (V - iν) k)i k √ν - iν ν

• Step 3 Solve the master problem with the additional constraint -∆ν e x - jxν e ∆ν min cTx + θν s.a. Ax ) b θν + βνk x g Rνk ,

k ) 1, ..., ν

-∆ν e x - jx e ∆ν ν

to obtain xν+1. Then calculate η ) cxν+1 + max1ekeν{Rνk + βνk xν+1} η¯ ) cxjν+1 + max1ekeν{Rνk + βνk jxν} F ) η - η¯ • Step 4 The algorithm stops if the estimated value of the objective function of the incumbent solution is less or equal to a specified tolerance.

|

|

γν - γν-1 eε γν 3.4. Computational Implementation. Figure 4 shows a diagram of the solution steps. The computational implementation of the algorithm integrates the sampling technique (coded in FORTRAN), the GAMS modeling environment,23 and a C++ manager program. The C++ program executes the main steps of the algorithm, transfers the control of execution, verifies the convergence and generates the appropriate MILP and LP problems for each of the SD iterations. Figure 5 depicts the implementation.

Figure 5. Computational implementation.

composition, L-shaped, and deterministic equivalent, to compare their solution to that obtained by the algorithm proposed here. Example 2 is an industrial size problem; it includes 38 processes, 28 chemicals, an international market and a local market. The problem involves 38 integer variables, 267 constraints, 202 continuous variables, and 33 stochastic parameters. All models were formulated and solved on an hp workstation xw6000 using GAMS with OSL and BDMLP as the LP and MILP solvers, respectively.

4. Illustrative Examples and Results The goals of the examples are to determine the best structure for a network of processes and the scheduling of capacity expansions given a limited budget and an increase in demand. The stochastic problems take into consideration uncertainties in the availability of raw material and in the demand of products. Example 1 is a small example with three chemicals and three processes. We used different solution methods, Benders de-

Figure 6. Network of example 1.

2818


Figure 7. Convergence results for example 1 with different sampling techniques. Table 1. Availabilities and Demands for Example 1

Table 6. Scenarios for Example 1

demand

chemical 1

chemical 2

chemical 3

scenarios

chemical 1

chemical 2

chemical 3

probability

lower upper

5400 6000

18000 20000

27000 30000

1 2 3 4 5 6 7 8

6000 6000 6000 6000 5400 5400 5400 5400

20000 20000 18000 18000 20000 20000 18000 18000

30000 27000 30000 27000 30000 27000 30000 27000

1/8 1/8 1/8 1/8 1/8 1/8 1/8 1/8

Table 2. Cost Parameters for Example 1 process

R in 102$/ton

β in 105$

δ in 102$/ton

1 2 3

1.38 2.72 1.76

85 73 110

0.40 0.60 0.50

chemical

Γ in 102$/ton

1 2

4.0 9.6

chemical

γ in 102$/ton

3

26.20

Table 3. Mass Balance Coefficients and Initial Capacities for Example 1 process

µ

Qi,0

1 2 3

1.0 1.0 1.0

0.0 10000 0.0

Table 4. Investment Data, Expansion Capacities, and Some Data to Calculate the NPV for Example 1 maximum investment capacity minimum expansion capacity maximum expansion capacity project lifetime interest rate tax rate working capital factor salvage value factor

$200,000 5000 ton 50000 ton 20 y 0.10 0.45 0.15 0.10

4.1. Example 1. The model was taken from the work of Sahinidis et al.24 Consider the superstructure shown in Figure 6 consisting of three processes that can be used to produce chemical 3 from raw materials chemical 1 and chemical 2. Given

only one time period with a 20-year horizon, the goal is to select the process(es) and the best schedule for capacity expansion to use in that time period, with uncertain demand for chemical 3 as well as uncertain availability for chemicals 1 and 2. All these uncertainties were represented by a continuous distribution of uniform probability. The initial capacity for process 2 is 10 000 ton/y with a spending limit of $200,000. The values for the uncertain parameters can be found in Table 1, the cost parameters are given in Table 2, mass balance coefficients and initial capacities are reported in Table 3. Finally, Table 4 provides the investment data, expansion capacities, and the additional data to calculate the NPV. First, the convergence results of the SD algorithm, with incumbent solutions using different sampling techniques are reported in Figure 7. The error upper bound was set to 1 × 10-6. Using latin hypercube (LHS) and Monte Carlo (MC) sampling techniques, the algorithm converges in 112 and 128 iterations, respectively. After 300 iterations, convergence was not achieved by using median latin hypercube (MLHS) and Hammersley sequence sampling techniques (HSS). Furthermore, the use of HSS does not result in asymptotic behavior. Table 5 shows the results obtained by the different solution methods. Note that to represent the uncertainties in the methods that use a discrete probability distribution, scenarios in terms of upper and lower bounds of the availability of the raw materials and the demand of the products were used (see Table 6), whereas for the incumbent SD algorithm a continuous

Table 5. Results by Different Methods for Example 1 variables

deterministic equivalent (8 scenarios)

Benders (8 scenarios)

L-shaped (8 scenarios)

incumbent SD (128 samples by MC)

iterations computational time (s)

-1.02322 × 108 Q1,1 ) 5000 Q2,1 ) 22904.412 QE1,1 ) 5000 QE2,1 ) 12904.412 y1,1 ) 1 y2,1 ) 1 1 0.10

-1.023 × 108 Q1,1 ) 5000 Q2,1 ) 22904.412 QE1,1 ) 5000 QE2,1 ) 12904.412 y1,1 ) 1 y2,1 ) 1 1 0.25

-1.02321 × 108 Q1,1 ) 5000 Q2,1 ) 22904.411 QE1,1 ) 5000 QE2,1 ) 12904.411 y1,1 ) 1 y2,1 ) 1 5 31

-1.01099 × 108 Q1,1 ) 5000.02 Q2,1 ) 22904.4 QE1,1 ) 5000.02 QE2,1 ) 12904.4 y1,1 ) 1 y2,1 ) 1 128 292

NPV Qit QEit yit


2819

Figure 8. Capacity of process 1 with respect to the number of iterations for Example 1. B-SD is the basic SD algorithm, and I-SD is the SD algorithm with incumbent solutions. Table 7. Results of the Optimal Configuration of Example 2 variables

incumbent SD algorithm (209 samples by MC)

NPV Qit

-21694.0572 Q1,1 ) 75.639, Q2,1 ) 64.361, Q3,1 ) 130.235, Q5,1 ) 27.059, Q6,1 ) 24.0, Q8,1 ) 24.0, Q12,1 ) 39.9, Q13,1 ) 25.0, Q14,1 ) 32.941, Q16,1 ) 300.0, Q17,1 ) 120.0, Q20,1 ) 76.0, Q26,1 ) 64.0, Q28,1 ) 77.96, Q29,1 ) 80.0, Q32,1 ) 216.0, Q34,1 ) 40.0, Q35,1 ) 76.458, Q38,1 ) 200.00 QE1,1 ) 75.639, QE2,1 ) 64.361, QE3,1 ) 130.235, QE5,1 ) 27.059, QE6,1 ) 24.0, QE8,1 ) 24.0, QE14,1 ) 32.941, QE16,1 ) 300.0, QE17,1 ) 120.0, QE20,1 ) 76.0, QE26,1 ) 64.0, QE28,1 ) 77.96, QE29,1 ) 80.0, QE32,1 ) 216.0, QE34,1 ) 40.0, QE35,1 ) 76.458 y1,1 ) 1, y2,1 ) 1, y3,1 ) 1, y5,1 ) 1, y6,1 ) 1, y8,1 ) 1, y14,1 ) 1, y16,1 ) 1, y17,1 ) 1, y20,1 ) 1, y26,1 ) 1, y28,1 ) 1, y29,1 ) 1, y32,1 ) 1, y34,1 ) 1, y35,1 ) 1

QEit

yit

Table 8. Comparison of Methods for Example 1 method

L-shaped

incumbent SD

continuous variables binary variables stochastic parameters iterations scenarios subproblems computational time (s)

6 3 3 5 8 40 31

6 3 3 128 128 512 292

Table 9. Comparison of Methods for Example 2 method

L-shaped

incumbent SD

continuous variables binary variables stochastic parameters iterations scenarios subproblems cmputational time (s)

202 38 33 not calculated 233 a (233)(iterations)a not calculated

202 38 33 209 209 836 993

a Assuming two values for discrete probabilities of uncertain parameters.

distribution of uniform probability was used to represent the uncertainties. It is interesting that, although the SD algorithm uses continuous distributions and the other techniques use discrete distributions, all the methods predict the same optimal configuration as well as same expansions capacities. The variation in the distributions, however, causes a difference in the estimated value for the NPV of about 1.2% between SD incumbent algorithm and the other methods. An additional issue is that the SD method has a stopping criterion that is statistical in nature, while the other has exact stopping criteria. Finally, we want to add that, when comparing the performance of the basic SD algorithm (B-SD) with the SD algorithm with incumbent solutions (I-SD), the I-SD is clearly advanta-

geous in terms of convergence. By using the same starting point and the same sequence of values of the uncertain parameters (obtained through the Monte Carlo sampling technique for each iteration), the B-SD algorithm did not achieve convergence after 300 iterations, whereas the I-SD converged in 128 iterations (as mentioned above). As an illustration, Figure 8 shows the optimal value of the capacity for process 1 obtained for both of the algorithms. I-SD rapidly reaches the optimal value of the capacity, whereas the B-SD algorithm fluctuates above and below such a value. In fact, the reason for lack of convergence of the B-SD algorithm is also illustrated in Figure 8, since the iterative procedure seems to lead to a non optimal value for the capacity of process 1. 4.2. Example 2. Consider the processes network proposed by Iyer and Grossmann given in Figure 9.17 The configuration consists of 38 processes, 28 chemicals, an international market and a local market. The goal of the problem is determine the best process configuration satisfying the product demands subject to the availabilities of raw materials, which are both uncertain and represented by a continuous distribution of uniform probability. Due to the large size of this problem, the data are not shown but are available upon request. The convergence results of the SD algorithm using different sampling techniques are reported in Figure 10. As in example 1, the error upper bound was defined as 1 × 10-6. Also, using the MC sampling technique, the SD algorithm has better convergence properties than the other techniques; it needed 209 iterations. Observe that, by using MC, HSS, and LHS sampling techniques, the algorithm predicts a similar expected NPV and similar optimal configurations (shown in Table 7). Finally, in Tables 8 and 9 we present the characteristics of the problems that were solved, as well as the computational effort required to solve each one of the proposed examples, through two different methodologies: one using an analysis of scenarios (like

2820


Figure 9. Network of example 2.

the L-shaped method) and another one using incumbent SD algorithm. In these results it is possible to observe that for small size problems the methods that use an analysis of scenarios make a smaller computational effort; nevertheless, when the size of the problem increases, the incumbent SD algorithm is preferred. 5. Conclusions A strategic capacity allocation model has been reformulated as a mixed-integer stochastic program with first stage integer

decisions. Uncertainties were incorporated in the availabilities of the raw material and in the demand of the products throughout continuous probabilistic distribution functions. An SD incumbent algorithm was proposed to solve the resulting two-stage stochastic problems with recourse. A procedure to identify an incumbent solutions sequence was incorporated into our implementation of the basic SD algorithm. Hence, the incumbent solutions were used to define a stopping criterion within a specified error. This implementation substan-


2821

Figure 10. Convergence results for example 2 with different sampling techniques.

tially reduces the number of iterations needed for convergence. Also, it was shown that the SD algorithm with incumbent solutions not only solves small problems but also large scale problems. In particular, a problem with 38 binary variables, 202 continuous variables, and 33 stochastic parameters was solved. Results also show better convergence properties when using the MC sampling technique. We expect then that this extension to the SD algorithm can be particularly suitable for large supply chain optimization problems, helping to make decisions about how a process supply chain should be designed, in the face of uncertainty. Acknowledgment The authors would like to acknowledge the financial support provided by the Mexican National Council for Science and Technology (CONACYT). Literature Cited (1) Grossmann, I. E. Challenges in the New Millennium: Product Discovery and Design, Enterprise and Supply Chain Optimization, Global Life Cycle Assessment. Comput. Chem. Eng. 2004, 29, 29–39. (2) Shah, N. Pharmaceutical Supply Chains: Key Issues and Strategies for Optimization. Proceedings FOCAPO 2003; Coral Springs, FL, January 12-15, 2003; Grossmann, I. E., McDonald, C. M., Eds.; CACHE, 2003; pp 73-85. (3) Masini, G.; Petracci, N.; Bandoni, A. Supply Chain Optimization in the Fruit Industry. Proceedings FOCAPO 2003; Coral Springs, FL, January 12-15, 2003; Grossmann, I. E., McDonald, C. M., Eds.; CACHE, 2003; pp 237-240. (4) Chopra, S.; Meindl, P. Supply Chain Management. Strategy, Planning and Operation; Prentice Hall: New Jersey, 2001. (5) Shapiro, J. Challenges of Strategic Supply Chain Planning and Modeling. Proceedings FOCAPO 2003; Coral Springs, FL, January 12-15, 2003; Grossmann, I. E., McDonald, C. M., Eds.; CACHE, 2003; pp 27-34. (6) Sahinidis, N. V. Optimization under Uncertainty: State of the Art and Opportunities. Proceedings FOCAPO 2003; Coral Springs, FL, January 12-15, 2003; Grossmann, I. E., McDonald, C. M., Eds.; CACHE, 2003; pp 153-165. (7) Barbaro, A. F.; Bagajewicz, M. Managing Financial Risk in Planning under Uncertainty. AIChE J. 2004, 50 (5), 963–989. (8) Clay, R. L.; Grossmann, I. E. A Disgregation Algorithm for the Optimization of Stochastic Planning Models. Comput. Chem. Eng. 1997, 21 (7), 751–774.

(9) Birge, J. R. Louveaux F. Introduction to Stochastic Programming; Springer-Verlag: New York, 1997. (10) Dantzig, G. B. Linear Programming Under Uncertainty. Manage. Sci. 1955, 1, 197–206. (11) Beale, E. M. L. On Minimizing a Convex Function Subject to Linear Inequalities. J. R. Stat. Soc. B17: 1955, 173–184. (12) Kall, P.; Wallace, S. W. Stochastic Programming; John Wiley and Sons: Chichester, England, 1994. (13) Haneveld, W. K.; Vander Vlerk, M. H. Stochastic Integer Programming: General Models and Algorithms. Ann. Oper. Res. 1999, 85, 39–57. (14) Higle, J. L.; Sen, S. Stochastic Decomposition; Kluwer Academic Publisher: Norwell, MA, 1996. (15) Ponce-Ortega, J. M.,; Rico-Ramı´rez, V.; Hernańdez-Castro, S.; Diwekar, U. M. Improving Convergence of the Stochastic Decomposition Algorithm by Using an Efficient Sampling Technique. Comput. Chem. Eng. 2004, 28, 767–773. (16) Rico-Ramı´rez, V.; Frausto-Hernańdez, S.; Diwekar, U. M.; Hernańdez-Castro, S. Water Networks Security: A Two-Stage Mixed-Integer Stochastic Program for Sensor Placement Under Uncertainty. Comput. Chem. Eng. 2007, 31, 565–573. (17) Iyer, R. R.; Grossmann, I. E. A Bilevel Decomposition Algorithm for Long-Range Planning of Process Networks. Ind. Eng. Chem. Res. 1998, 37 (2), 474. (18) Simchi-Levi, D.; Kaminsky, P.; Simchi-Levi, E. Designing and Managing the Supply Chain; Irwin/McGraw-Hill: New York, 1998. (19) Dantzig, G. B. Planning Under Certainty Using Parallel Computing. Ann. Oper. Res. (Switzerland) 1988, 14, 1–4 (A01) 1-16. (20) Bienstock, D.; Shapiro, J. F. Stochastic Programming Models for Strategic Planning: An Application to Electric Utilities; MIT Technical Report OR 128-84, 1985; revised. (21) Painton, L. A.; Diwekar, U. M. Stochastic Annealing under Uncertainty. Eur. J. Oper. Res. 1995, 83, 489. (22) Elshorbagy, W. E.; Yakowitz, D. S.; Lansey, K. E. Design of Engineering Systems Using Stochastic Decomposition: Water Supply Planning Application. Eng. Optim. 1997, 27 (4), 279–302. (23) Brooke, A.; Kendrick, D.; Meeraus, A.; Raman, R. GAMS a User’s Guide.4; GAMS Development Corporation, 1998. (24) Sahinidis, N. V.; Grossmann, I. E.; Fonari, R. E.; Chathrati, M. Optimization Model for Long Range Planning in the Chemical Industry. Comput. Chem. Eng. 1989, 13 (9), 1049.

ReceiVed for reView June 23, 2009 ReVised manuscript receiVed January 19, 2010 Accepted January 21, 2010 IE901021B

Strategic Capacity Allocation under Uncertainty by Using a Two-Stage

Recommend Documents