474
Ind. Eng. Chem. Res. 1998, 37, 474-481
A Bilevel Decomposition Algorithm for Long-Range Planning of Process Networks Ramaswamy R. Iyer and Ignacio E. Grossmann* Department of Chemical Engineering, Carnegie Mellon University, Pittsburgh, Pennsylvania 15213
The solution of the multiperiod MILP model for long-range planning of process networks by Sahinidis et al. (Comput. Chem. Eng. 1989, 13, 1049) is addressed in this paper. The model determines the optimal selection and expansion of processes over a long-range planning horizon, incorporating multiple scenarios for varying forecasts for demands and prices of chemicals. A rigorous bilevel decomposition algorithm is proposed to reduce the computational cost in the multiperiod MILP model. The decomposition algorithm solves a master problem in the reduced space of binary variables to determine a selection of processes and an upper bound to the net present value. A planning model is then solved for the selected processes to determine the expansion policy and a lower bound to the objective function. Numerical examples are presented to illustrate the performance of the algorithm and to compare it with a full-space branch and bound method. 1. Introduction A number of papers have been published on the capacity expansion problem for several industrial applications [see Roberts (1964), Manne (1967), Noonan and Giglio (1977), Luss (1982), Minoux (1987), Sahinidis et al. (1989)]. The industrial relevance of these problems is discussed in Florian et al. (1980) and Hirshfeld (1987). In this work, we address the specific application of expansion planning in the chemical industry. Given a network of processes and chemicals, the objective is to select processes and their capacity expansion policy for a given forecast of prices and demands over the planning horizon. The objective function is the maximization of the net present value. The problem is characterized by binary decision variables for each process in each time period. Thus, the combinatorial complexity of the problem increases with the number of processes and periods, with which large-scale problems may become computationally intractable. To solve the above problem, Sahinidis et al. (1989) investigated several strategies including branch and bound, use of integer cuts, strong cutting planes, Benders decomposition, and heuristics. Later, Sahinidis and Grossmann (1991) proposed a reformulation based on lot sizing substructures, leading to a problem with more constraints and variables, but with a tighter LP relaxation, which potentially reduces the number of nodes searched in the branch and bound tree. Planning problems of up to 38 processes, 26 chemicals, and 4 time periods were solved with reasonable computation times. However, for larger numbers of processes and time periods, the number of constraints and variables in the reformulated model becomes very large. To circumvent this problem, Liu and Sahinidis (1996a) developed projection techniques in order to reduce the size of the reformulated model. In addition, the uncertainty in demands and prices may be modeled using discrete random parameters, leading to multiscenario multiperiod problems [see Liu and Sahinidis (1996b)]. For * Author to whom all correspondence should be addressed. Telephone: (412)-268-2230. E-mail:
[email protected].
uncertainty in m parameters with n discrete values in each time period, there are a total of (nm × NT) scenarios. The use of full-space techniques is computationally expensive if not impossible for such large problems, motivating the development of efficient decomposition solution techniques. In this work, we address the multiperiod MILP planning model with multiple scenarios for each time period. Given that a network of processes is given, the decision variables involved in the problem are as follows: 1. Capacity expansion of existing processes. 2. Selection of new processes and their capacity expansion policy. 3. Production profiles and sales and purchases of chemicals in each time period and scenario. The paper is organized as follows. In the following section, the problem formulation is reproduced as in Sahinidis and Grossmann (1991)). The multiperiod MILP model is modified to include multiple scenarios for each time period. The proposed decomposition algorithm, which produces the same optimal solution as a full-space MILP method, is explained in section 3 followed by example problems in section 4. The performance of the proposed algorithm is compared to the reformulation algorithm as in Sahinidis and Grossmann (1991). 2. Problem Formulation The network superstructure of NP chemical processes and NC chemicals (including raw materials, intermediates, and products) is given. There are NK streams in the network, and the material balance for each process is expressed linearly based on the production rate of the main products from the process. The capacity of the process is determined by the flow rate of the main product. The demands and prices of chemicals may vary for each time period over the planning horizon. It is assumed that NS independent scenarios are considered for each time period, which are obtained from discrete random uncertainty in demands and prices. This is equivalent to considering a two-stage stochastic programming approach as described in Hiller (1986).
S0888-5885(97)00383-7 CCC: $15.00 © 1998 American Chemical Society Published on Web 01/14/1998
Ind. Eng. Chem. Res., Vol. 37, No. 2, 1998 475
The objective function to be maximized is the expected net present value (NPV) of the project over the planning horizon consisting of NT time periods. The NPV includes the discounted investment costs, operating costs, and revenues from sales in NM markets. The indices, parameters, and variables defined in the model are as follows. 1. Indices
i ) process (i ) 1, ..., NP)
yit
) decision variable for expansion of process i in period t ) 1 if there is an expansion, 0 otherwise
Qit ) capacity of process i in period t QEit ) capacity expansion of process i in period t Pjlts ) amount of product j purchased in period t for scenario s in market l Sjlts ) amount of product j sold in period t for scenario s in market l
j ) chemical (j ) 1, ..., NC) t ) time period (t ) 1, ..., NT)
Wkts ) amount of flow of stream k in period t for scenario s
k ) stream in network (k ) 1, ..., NK) l ) market (l ) 1, ..., NM)
4. Superscripts
s ) scenarios (s ) 1, ..., NS)
( )L ) lower bound 2. Parameters
( )U ) upper bound
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 ∈ Li)
µik
) material balance coefficient for process i and stream k
πst
) probability of scenario s in period t s.t.
∑sπst ) 1
∀t
The equations in model P are as follows: NT NS NC NM
P:
∑ ∑∑∑πstγjltsSjlts t)1s)1j)1 l)1
NT NP
∑ ∑ t)1 i)1
NT NS NP
(RitQEit + βityit) -
∑ ∑∑πstδitsWm ts t)1s)1i)1 i
NT NS NC NM
∑ ∑∑∑πstΓjltsPjlts t)1s)1j)1 l)1
(1)
subject to yitQELit e QEit e yitQEU it Qit ) Qi,t-1 + QEit
NEXP(i) ) maximum number of allowable expansions of the process (i) CI(t)
max NPV )
i ) 1, ..., NP, t ) 1, ..., NT (2)
i ) 1, ..., NP, t ) 1, ..., NT (3)
NT
yit e NEXP(i) ∑ t)1
) capital investment limit in period t
i ∈ I′ ⊆ {1, 2, ..., NP}
(4)
NP
For each scenario s and period t
Rits ) variable term for investment cost of process i ($/unit capacity)
(RitQEit + βityit) e CI(i) ∑ i)1
γjlts ) price of sales of chemical j in market l ($/unit sold) Γjlts ) price of chemical j in market l ($/unit purchased) ajlts ) bound for sales of chemical j in market l ($/unit sold) djlts ) bound for purchase of chemical j in market l ($/unit purchased) 3. Variables
(5)
For t ) 1, ..., NT, s ) 1, ..., NS
βits ) fixed term for investment cost of process i ($) δits ) unit operating cost of process i ($/unit production amount)
t ∈ T′ ⊆ {1, 2, ..., NT}
Wmits e Qit NM
i ) 1, ..., NP
(6)
NM
Pjlts + ∑ Wkts ) ∑Sjlts + ∑ Wkts ∑ l)1 l)1 k∈I(j)
k∈O(j)
j ) 1, ..., NC (7) Wkts ) µikWmits
k ∈ Li/mi
aLjlts e Pjlts e aU jlts
j ) 1, ..., NC,
dLjlts e Sjlts e dU jlts
j ) 1, ..., NC,
(8)
l ) 1, ..., NM (9) l ) 1, ..., NM (10)
QEit, Qit, Wkts, Pjlts, Sjlts g 0 yit ∈ {0, 1}
(11)
476 Ind. Eng. Chem. Res., Vol. 37, No. 2, 1998
Figure 1. Flowchart for the bilevel decomposition algorithm.
The objective function in eq 1 is the expected net present value, which includes the expected revenues from all t over all scenarios s and the investment costs for all t and sum of operating costs and cost of purchase of chemicals from all t over all scenarios s. 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 for process i over the planning horizon. Equation 5 limits the amount of capital investment in each time period t. Equation 6 imposes production limits based on existing capacity, and eq 7 is a mass balance for each chemical in the complex. Equation 8 is the mass balance of chemicals in each process i, and eqs 9 and 10 bound the amount of purchase and sales of chemicals. These bounds could arise, for instance, in the case of long-term contracts in which fixed amounts of purchases or sales are committed over several time periods. Note that eqs 6-11 are indexed for each possible scenario s for given time period t. The scenarios are treated independently for each time periods and are obtained from discrete sampled values of random parameters. Clearly, the size of problem P increases with the number of scenarios and number of time periods. The combinatorial complexity increases with the number of time periods and number of processes. As a result, fullspace techniques can fail to solve to optimality large multiperiod planning problems. Reformulation techniques (see Sahinidis and Grossmann (1991)) to tighten the relaxation gap result in an increase in the number of equations and variables. Planning problems with large numbers of time periods, scenarios, and processes become intractable even with reformulation due to the potentially large size of resulting MILP’s. Although the projection method of Liu and Sahinidis (1996) reduces the problem size, it can weaken the LP relaxation. We propose a bilevel decomposition algorithm based on the algorithm presented in Iyer and Grossmann (1998) for the design and multiperiod operation of utility plants. The algorithm allows incorporation of the reformulation technique and enables the solution of large-scale problems at reasonable computational expense. 3. Decomposition Algorithm The proposed decomposition algorithm (see Figure 1) solves a higher level design master problem (DP), which
is a specific relaxation of problem P, to obtain an upper bound for the NPV. The design master problem does not contain binary variables yit associated with capacity expansion decisions. It only contains binary variables Yi representing the selection of a process over the entire planning horizon. Problem DP is therefore combinatorially less complex and selects a subset of processes for design. The lower level planning problem (OP) is solved for the selected set of processes (process i with Yi ) 1), yielding a lower bound to the NPV for any feasible solution of OP. Since a subset of processes is selected for planning in OP, it contains fewer decision binary variables yit and is not as combinatorially complex as problem P. We note that for both problems DP and OP, the equations for all time periods and scenarios are included for all the processes that are considered. The computational expense is lowered by reducing the number of binary variables in each level. The problems are solved iteratively until the bounds converge. The higher level design problem master problem (DP) is defined as follows: NT NS NC NM
DP:
max NPV )
∑ ∑∑∑πstγjltsSjlts t)1s)1j)1 l)1
NT NP
∑ ∑ t)1i)1
NT NS NP
(RitQEit + βitQEit/QEU it ) -
∑ ∑∑πstδitsWm ts t)1s)1i)1 i
NT NS NC NM
∑ ∑∑∑πstΓjltsPjlts t)1s)1j)1 l)1
(12)
subject to QEit e YiQEU it Qit ) Qi,t-1 + QEit
i ) 1, ..., NP, t ) 1, ..., NT
(13)
i ) 1, ..., NP, t ) 1, ..., NT (3) eqs 6-10
QEit, Qit, Wkts, Pjlts, Sjlts g 0 Ω(Yi) ) 0
(14)
Yi ∈ {0, 1}
(15)
In model DP, eqs 4 and 5 are removed. In addition, the lower bound constraints on QEit in (2) are relaxed. The objective function is overestimated since
QEit/QEU it e yit
(16)
Thus, for every expansion,
RitQEit + βitQEit/QEU it e RitQEit + βityit
(17)
The investment costs are therefore underestimated. Also, since the investment costs have a negative coefficient, the objective function is overestimated. Thus, model DP is a relaxation of P and overestimates the objective function, yielding an upper bound on NPV. Therefore, the following property then applies to DP. Property 1. Problem DP yields an upper bound to the solution of problem P. In DP, eq 14 represents logic cuts derived from the network structure (see Raman and Grossmann (1993)). As an example, if the product from process A is the sole feed to process B, then a logical cut may be symbolically written as
Ind. Eng. Chem. Res., Vol. 37, No. 2, 1998 477
process B w process A
(18)
which in turn can be mathematically written as
YA - YB g 0
(19)
Since variable Yi represents the existence of the process, these logic cuts are valid constraints (see Hooker et al. (1994)). Similar logic cuts cannot be written for yit in P as they represent expansion decision variables in each time period. The following additional properties can be established. Property 2. Let the design solution of DP be Y h i. For fixed Y h i, a feasible solution to OP is a feasible solution to P and yields a lower bound to P, where OP is defined as
OP: problem P subject to yit e Y hi
(20)
Clearly, eq 20 has the effect of choosing a subset of processes for the planning problem, and therefore a feasible solution to (OP) will be feasible in P, yielding a lower bound to the NPV. Property 3. Let QEr be the optimal solution to OP for a corresponding fixed configuration Y h ri . For positive cost coefficients of QEit in the objective function, the following inequality is a valid design cut for problem DP,
Yri - ∑ Yri - (|Mr| - 1)] ∑ i∈M i∈N
QEit g QErit[
r
∀i, t (21)
r
where
h ri ) 1 for configuration in iteration r} Mr ){i|Y h ri ) 0 for configuration in iteration r} Nr ) {i|Y Proof: For given Y h ri , the solution to OP gives an r optimal design QEit which represents the best design for that choice of configuration Y h ri in the original problem P. For the selected value Y h ri , the multiplying factor for the right-hand side for QErit is 1. Therefore, the inequality QEit g QErit is enforced through (21). This inequality also corresponds to a valid cut to DP due to the positive cost coefficients. Furthermore, for h ri , QEit g 0 any other choice of configuration Yri * Y dominates eq 21, since
Yri - ∑ Yri - (|Mr| - 1)] e 0 ∑ i∈M i∈N
[
r
(QED)
r
Note that when the same configuration is obtained, then due to the positive cost coefficient of QEit the value of the NPVDP is lowered as QEit g QErit from (21), thereby tightening the difference between the lower and upper bounds. The following property establishes the basis for deriving integer cuts on supersets and subsets of the configurations predicted by the design problem DP. Property 4. Let Zr1 ) {i|Yi )1} and Zr0 ) {i|Yi ) 0} correspond to the optimal solution of DP in iteration r. For all iterations s > r, if Zr1 is feasible in OP, then any solution Zs1 ⊃ Zr1 will result in a solution of OP such
that NPVsOPe NPVrOP. Furthermore, the subsets Zs1 ⊂ Zr1 can be excluded from DP at iteration s > r. Proof: Consider Zr1, the solution of DP in iteration r. If any feasible solution Zs1 ⊇ Zr1 is not selected before iteration r, it implies that in the relaxed problem DP, the selection of any additional process u does not result in an increase in NPV in DP. Therefore,
NPVsDP e NPVrDP
(22)
Now assume that there exists set Zs1 ⊃ Zr1 such that NPVsOP > NPVrOP in OP. This implies that the investment in the additional process u ∈ Zs1 and u ∉ Zr1 results in an increase in NPVOP. Since the investment costs are underestimated in DP and since DP is a relaxation of OP, the additional investment in process u would result in an increase in NPVDP which implies NPVsDP > NPVrDP. This contradicts eq 22 and therefore the property holds for Zs1 ⊃ Zr1. Hence, adding the additional process will not result in a larger NPV in problem OP, which is a restriction of DP, and therefore NPVsOP e NPVrOP. Now consider any Zs1 ⊂ Zr1. Zs1 may either be feasible or infeasible in OP. If Zs1 is infeasible in OP, then we may exclude all these infeasible solutions for consideration in problem DP. Consider only those solutions in Zs1 that are feasible in OP. Since the optimal solution of OP will consider as candidates Zr1 and all possible subsets Zs1, it follows that Zs1 can be excluded from further consideration in DP for all iterations s > r. QED Thus, for Zs1 ⊃ Zr1, the solution obtained from DP is used to exclude configurations that would be suboptimal in OP. For Zs1 ⊂ Zr1, the best solution amongst all subsets of Zr1 will be chosen by OP and therefore, if any subset of Zr1 is chosen by DP, it will have been considered in OP in the previous iteration. The implication of this property is that an integer cut may be added to DP at each iteration s > r, which precludes all supersets and subsets of Zr1. This will help in eliminating a large number of feasible configurations from the solution of DP and therefore reduces the search space for binary variables significantly. The cut for precluding supersets may be logically written as
if ∧ r Yn w ∧ r (¬Yi) n∈Z1
i∈Z0
(23)
The logical proposition in (23) states that if in any solution all the y’s in any set Zr1 are 1, then all other y’s must be zero in order to prevent a superset of Zr1 from entering the solution of DP. Mathematically, this equation may be written as
∑ Yn + Yi e |Zr1|
∀i ∈ Zr0
(24)
n∈Zr1
Similarly, the cuts for precluding subsets of Zr1 are equivalent to the condition,
( ∨ Y )∨Y n∈Zr0
n
i
∀i ∈ Zr1
which gives rise to the equation
(25)
478 Ind. Eng. Chem. Res., Vol. 37, No. 2, 1998
∑ Yn + Yi g 1
∀i ∈ Zr1
(26)
n∈Zr0
A description on deriving cuts from logical inference clauses is given in Raman and Grossmann (1991). 3.1. Remarks. Problem OP is essentially problem P for a subset of selected processes from the superstructure. At the level of problem OP, the reformulation techniques proposed by Sahinidis and Grossmann (1991) may be used to reduce the relaxation gap. The projection approach of Liu and Sahinidis (1996a) may also be used at this level to improve the computational efficiency. The bilevel decomposition primarily aims at reducing the size of the planning subproblems at each stage by proposing design configurations based on a design master problem. It should also be noted that the proposed decomposition algorithm produces the same optimal solution as when a full-space MILP method like branch and bound is used. 3.2. Algorithmic Steps. Based on the above properties, the steps of the bilevel decomposition algorithm are presented below. (1) Set iteration count r ) 0, upper bound UB ) ∞, and lower bound LB ) -∞. (2) r ): r + 1. Set iteration number R ) r. Solve DP to get a design Y h ri with upper bound UB. Define
h ri ) 1} Mr ) {i|Y h ri ) 0} Nr ) {i|Y Yri
Y h ri ,
(3) For fixed ) solve OP [i.e., problem P with inequalities in (20)] to obtain a capacity expansion plan and a lower bound LBr. Set LB ) maxr{LBr}. (a) If OP is infeasible, add the following integer cut to DP.
Yri - ∑ Yri e |Mr| - 1 ∑ i∈M i∈N r
(27)
r
Go to step 2. (b) If OP is feasible, then define
h i ) 1} Zr1 ) {i|Y h i ) 0} Zr0 ) {i|Y Add to the master problem DP the design cut as in (21)
QEit g QErit [
∑ Yri - i∈N ∑ Yri - (|Mr| - 1)]
i∈Mr
r
r ) 1, ..., R, i ) 1, ..., NP, t ) 1, ..., NT and the integer cuts (24) and (26) as follows:
∑ Yn + Yi e |Zr1|
∀i ∈ Zr0
r ) 1, ..., R
n∈Zr1
∑ Yn + Yi g 1
∀i ∈ Zr1
r ) 1, ..., R
n∈Zr0
(4) If (UB - LB)/UB e (tolerance), stop. The solution corresponding to LB is the optimal solution. Else, go to step 2. Provided that the number of major iterations is modest, the computational requirements of the proposed
Figure 2. Example 1 network superstructure. Table 1. Input Data for Example 1: 3 Processes, 3 Chemicals, 4 Time Periods time period process
1
2
3
4
102
1 2 3
Variable Investment Costs (R) in $/ton 1.38 1.67 2.22 2.720 3.291 4.381 1.760 2.130 2.834
1 2 3
Fixed Investment Costs (β) in 105 $ 85.00 102.85 136.89 73.00 88.33 117.56 110.00 133.10 177.15
1 2 3
Operating Costs (δ) in 102 $/ton 0.40 0.48 0.64 0.60 0.72 0.96 0.50 0.60 0.80
3.58 7.055 4.565 220.46 189.34 285.31 1.03 1.55 1.29
method can be expected to be lower than a full-space method as it involves the solution of smaller subproblems. In contrast, the full-space method may require a very large amount of time to solve to optimality or even for finding a feasible solution. Furthermore, if it is not solved to optimality, the full-space method may not produce good suboptimal solutions due to the potentially large size of the search tree and the cost of solving large LP’s at each node. 4. Examples In this section, three example problems are solved to illustrate the performance of the algorithm as compared to a full-space method (with and without reformulation). The example problems are of different levels of complexity based on the number of processes in the superstructure and number of time periods. 4.1. Example 1. A 3 process, 3 chemical problem (as in Sahinidis et al. (1989)) is first considered (see Figure 2) for a planning horizon consisting of 4 time periods, with a 15 year planning horizon. The length of the time periods are 2, 3, 5, and 5 yr, respectively. The data for this problem are presented in Tables 1 and 2, which is representative of the input data for all example problems. The multiple scenarios for the problem were generated as follows. A normal distribution for the demands of the 3 chemicals was used with the nominal value as the mean and a standard deviation equal to 10% of the mean. Two discrete points were sampled randomly from this distribution for each demand, yielding 8 scenarios for each time period, or a total of 32 scenarios over the entire planning horizon. Table 2 shows nominal values of the uncertain parameters (demands and prices). All the cost parameters are
Ind. Eng. Chem. Res., Vol. 37, No. 2, 1998 479 Table 2. Input Data for Example 1 time period chemical
1
2
3
4
1 2
Purchase Upper Bounds (P) in kton/yr 6.00 7.26 9.66 20.00 24.20 32.21
15.56 51.87
3
Sales Upper Bounds (S) in kton/yr 30.00 36.3 48.31
77.81
Cost of Purchase of Chemicals (Γ) in 102 $/ton 1 4.00 4.84 6.44 10.37 2 9.60 11.61 15.46 24.90 Sales Price of Chemicals (γ) in 102 $/ton 26.20 31.7 42.19
3
67.95
Table 3. Results for Example 1 full space size variables solution LP timeb equa($ 106/yr) relaxation (CPU s) nodes tions cont. 0-1
method standard reformulated
444.7 444.7
problem OP problem DP
464.7 448.8
0.3 1.0
23 4
416 884
325 12 553 12
Decomposition Algorithm 1.4 6 0.8 4
415 515
325 12/8a 319 3
Solution of Subproblems at Each Iteration
b
subproblem
iteration 1
iteration 2
design problem planning problem
457.9 444.7
365.1 algorithm terminated
a Largest/smallest number of binary variables over all iterations. Solution time on HP 9000/700.
Table 4. Results for Example 2
Figure 3. Example 1 results.
full space solution LP timeb method ($ 106/yr) relaxation (CPU s) standard reformulated
23.9 23.9
77.6 44.6
size variables
nodes
616 44 730 10 100 42 890
equations cont. 0-1 1747 1341 60 4327 2601 60
Decomposition Algorithm problem 135 3000/12a 1556 1231 50/20a OP problem 46 38/10a 1724 1293 6 DP optimal solution of 23.9 obtained in 16 major iterations b
a Largest/smallest number of binary variables over all iterations. Solution time on HP-9000/700.
discounted at a specified interest rate and include the effect of taxes. In this example, process 2 had an existing capacity of 50 ktons/yr at the start of the first time period. The problem was modeled using GAMS (Brooke et al. (1992)) and solved in the full space using the CPLEX solver. The optimal solution of $444.7 million was obtained in 0.3 CPU s after branching at 23 nodes. The problem was solved using the lot-sizing reformulation by Sahinidis and Grossmann (1991), yielding the same solution in 1 CPU s after branching on 4 nodes. Note that the reformulated MILP was solved with a fewer number of nodes but took longer because of the larger MILP and the time required in determining bounds on the maximum flow (Sahinidis and Grossmann (1991)). The problem was also solved using the bilevel decomposition, yielding the same optimal solution in 2 major iterations in 2.2 CPU s. The results are summarized in Table 3. For the small example problem, the fullspace method is more efficient as the number of decision binary variables is very small.
Figure 4. Example 2 network superstructure.
The solution obtained for this problem is shown in Figure 3, indicating the capacity and production rate (maximum value of production rate over all scenarios) for the processes in each time period. The capacity expansions for processes 1 and 3 were chosen in the first and second period, respectively, as indicated by the values of the total capacities for each process. Thus, the solution from the model indicates that its NPV would be larger if process 3 is installed only after 2 yr as the demand for chemical 3 is completely satisfied by the existing capacity of process 2. However, as the
480 Ind. Eng. Chem. Res., Vol. 37, No. 2, 1998
Figure 5. Example 3 network superstructure.
demand increases, the smaller operating costs of process 3 offsets its investment cost, leading to a larger NPV by utilizing process 3. 4.2. Example 2. Consider the 6 process, 10 chemical problem shown in Figure 4. A 20-yr planning horizon is assumed with 10 time periods, each of a length of 2 yr. The demand of chemicals acrylonitrile and isopropyl alcohol is uncertain, and 2 discrete values are sampled for each demand for each time period, yielding 4 scenarios per time period, or a total of 40 scenarios over the entire planning horizon. The computational results are presented in Table 4. The same optimal solution of $23.9 million was obtained by both the full-space method and the proposed algorithm in 616 and 181 CPU s, respectively. The resulting solution chose processes 4, 5, and 6 with expansions in periods 1 and 8. In this example, the decomposition algorithm obtained the optimal solution in less than 30% of the time used by the full-space method. The optimal solution was obtained in 16 major iterations. A significant
Table 5. Example 3 Problems Problem Data name PLAN6 PLAN10 PLAN16 PLAN38
time uncertain scenarios processes chemicals periods params per period 6 10 16 38
10 12 16 28
10 10 10 10
2 1 2 3
4 2 4 8
length of time period ) 2 yr Problem Size name
equations
continuous variables
binary variables
PLAN6 PLAN10 PLAN16 PLAN38
1741 1301 3361 13381
1341 1061 2721 11221
60 100 160 380
number of configurations were eliminated in the design problem due to the integer cuts in (24) and (26). This results in faster convergence between the bounds and fewer number of iterations. Also, using the lot-sizing
Ind. Eng. Chem. Res., Vol. 37, No. 2, 1998 481 Table 6. Example 3 Results full space
bilevel decomposition
problem
solution (M$)
LP relaxation
timea
PLAN6 PLAN10 PLAN16 PLAN38
59.9 138.6 269.5 570.6
132.0 217.7 367.9 789.0
39 123 3 750 30 000b
a
nodes
solution (M$)
timea
% time for OP
major iteration
2602 8812 50 000b 22 000
59.9 138.6 279.8 623.0
76 110 2 275 23 000
71 89 79 79
5 7 20b 20b
Time in CPU seconds on HP9000/700. b Search stopped.
reformulation of the problem in full space resulted in a somewhat tighter relaxation gap, but at the expense of a much larger MILP. As a result, the computational expense for the reformulated problem (10 100 s vs 616 s) is much larger due to the large size of LP’s at each node of the branch and bound tree. 4.3. Example 3. Consider the 38 process, 28 chemical problem shown in Figure 5. A set of problems (see Table 5) was run using this superstructure as a basis. For example, a 10 process problem (PLAN10) uses the superstructure derived from the first 10 processes in Figure 5 and involves 2 scenarios for each of the 10 time periods. PLAN38 involves 38 processes, 28 chemicals, 10 time periods with 9 scenarios each, giving rise to a MILP problem with 380 0-1 variables, 11 221 continuous variables, and 13 381 equations. The problems were solved using both a full-space method and the proposed decomposition algorithm. The results are presented in Table 6. For the smaller problem (PLAN6), the full-space method performs better. However, for the larger problems (PLAN16 and PLAN38), the decomposition algorithm produces solutions with higher NPV as compared to the full-space method (279.8 vs 269.5 and 623.0 vs 570.6, respectively). In both cases, the search had to be terminated due to time limitations (8.3 h of full space; 6.4 h for the proposed method). Since the decomposition method evaluates solutions for different choices of configurations sequentially and has fewer open nodes at the design level, it performs a more intelligent search in the feasible space of the binary variables. Note that in the full-space method the search does not reflect the natural hierarchy between process selection and capacity expansion due to the use of the single binary variable yit. The problems were also attempted using reformulation, resulting in tighter relaxation gaps at the expense of a larger number of equations and variables. For problems PLAN10, PLAN16, and PLAN38, the amount of time taken by the reformulated problem was very large due to the quadratic increase in the size of the problem with respect to the number of time periods. 5. Conclusions The solution of the multiperiod network planning problem (as in Sahinidis et al. (1989)) was considered, incorporating uncertainty in demands and prices through the use of scenarios. The resulting model is a large multiscenario, multiperiod planning problem with a binary decision variable in each time period, making it computationally intractable for full-space methods. The bilevel decomposition method proposed by Iyer and Grossmann (1998) for utility systems was extended to solve the problem within a bounding iterative strategy. Special integer cuts were proposed for network problems to reduce the number of feasible configurations evaluated by the design master problem. For very small problems, the method is no faster than the full-space branch and bound method. For medium-sized problems
that can be solved to optimality, the proposed method proved to be faster. For large problems, solutions with 10% higher solutions were obtained when compared to the suboptimal solution obtained by the full-space method. Acknowledgment The authors acknowledge support from a grant from DOE under Grant No. DE-FG02-85ER13396. Literature Cited Brooke, A.; Kendrick, D.; Meeraus, A. GAMS: A User’s Guide; Scientific Press: Palo Alto, CA, 1992. Florian, M. K.; Lenstra, J. K.; Rinnooy Kan, A. H. G. Deterministic production planning: Algorithms and complexity. Manage. Sci. 1980, 26, 669. Hiller, R. S. Stochastic programming approximation methods with application to multistage production planning. Ph.D. Thesis Dissertation, Operations Research Center, MIT, 1986. Hirshfeld, D. S. Mathematical programming and planning, scheduling and control of process operations. Paper presented at the FOCAPO Conference, Park City, 1988. Hooker, J. N.; Yan, H.; Grossmann, I. E.; Raman, R. Logic cuts for processing networks with fixed charges. Comput. Oper. Res. 1994, 21, 265-279. Iyer, R. R.; Grossmann, I. E. Synthesis and Operational Planning of utility systems for Multiperiod Operation. Comput. Chem. Eng. 1998, in press. Liu, M. L.; Sahinidis, N. V. Optimization of process planning under uncertainty. Ind. Eng. Chem. Res. 1996a, 35, 4154-4165. Liu, M. L.; Sahinidis, N. V. Long range planning in the process industry. Comput. Oper. Res. 1996b, 23, 237-253. Luss, H. Operations research and capacity expansion problems: a survey. Oper. Res. 1982, 30, 907. Manne, A. S., Ed. Investments for capacity expansion; The MIT Press: Cambridge, MA, 1967. Minoux, M. Network Synthesis and Dynamic Network Optimization. Ann. Discr. Math. 1987, 31, 283. Nemhauser, G. L.; Wolsey, L. A. Integer and Combinatorial Optimization; Wiley: New York, 1988. Noonan, F.; Giglio, R. J. Planning electric power generation: A nonlinear mixed integer model employing Benders decomposition. Manage. Sci. 1977, 23, 946. Norton, L.; Grossmann, I. E. Strategic Planning model for complete process flexibility. Ind. Eng. Chem. Res. 1994, 33, 69. Raman, R.; Grossmann, I. E. Relation between MILP modeling and logical inference for chemical process synthesis. Comput. Chem. Eng. 1991, 15, 73. Raman, R.; Grossmann, I. E. Symbolic Integration of Logic in Mixed Integer Linear programming techniques for Process Synthesis. Comput. Chem. Eng. 1993, 17, 909. Roberts, S. M. Dynamic Programming in Chemical Engineering and Process Control; Academic Press: New York, 1964. Sahinidis, N. V.; Grossmann, I. E. Reformulation of Multiperiod Optimization model for long range planning in the chemical industry. Comput. Chem. Eng. 1991, 15, 255. Sahinidis, N. V.; Grossmann, I. E.; Fornari, R. E.; Chathrathi, M. Optimization model for Long range planning in the chemical industry. Comput. Chem. Eng. 1989, 13, 1049.
Received for review June 2, 1997 Revised manuscript received November 7, 1997 Accepted November 10, 1997 IE970383I
