+
+
Ind. Eng. Chem. Res. 1996, 35, 2611-2623
2611
Disjunctive Programming Techniques for the Optimization of Process Systems with Discontinuous Investment Costs-Multiple Size Regions Metin Tu 1 rkay and Ignacio E. Grossmann* Department of Chemical Engineering, Carnegie Mellon University, Pittsburgh, Pennsylvania 15213
This paper addresses the optimization of process models that involve discontinuous investment cost functions with fixed charges and are defined over several regions for the sizes. These discontinuous cost functions are naturally expressed by disjunctions. Conventional modeling and solution techniques for the optimization and synthesis of process systems with these cost models include a direct NLP approach, the use of smoothing functions and MINLP models with big-M constraints. In this paper we propose the application of disjunctive programming approaches based on the convex hull formulation of disjunctions, and disjunctive branch and bound with surrogates or linear underestimators. Theoretical comparison between the convex hull formulation, the big-M model, and the model with linear underestimators is presented. It is proved that the convex hull formulation of disjunctions gives the tightest relaxation among the alternative modeling techniques. The proposed solution techniques are tested on heat exchanger networks, process flowsheet design, and the synthesis of a large petrochemical complex. It is shown that the convex hull formulation of disjunctions gives the overall best performance among the proposed techniques. 1. Introduction The synthesis and optimization of process systems require the development of mathematical programming models that incorporate the selection of units and design and operating parameters so as to optimize a cost function satisfying physical constraints and specifications. Operating costs are in most cases simple linear functions of flowrates and power and heat consumptions. Investment cost models on the other hand are considerably more complex. They usually consist of fixed costs and variable costs in terms of the design variables (e.g., area for a heat exchanger, volume of reactor). The current practice in modeling these investment costs is to use a single correlation assuming a continuous linear or nonlinear function over the entire range of the design variable, and with or without a fixed charge (Biegler, 1992). However, realistic cost models such as Guthrie’s cost correlations (Guthrie, 1974) are often complex discontinuous functions defined over multiple regions for sizes, pressures, and temperatures. Tu¨rkay and Grossmann (1996) have shown that generalized disjunctive programming models are very useful in modeling nonlinear discrete continuous optimization problems. They have shown benefits for the modeling and solution of flowsheet synthesis problems in which discrete decisions for the selection of process units have to be made among several alternatives. In this paper, it is shown that the disjunctive programming approach can be effectively applied in the modeling and solution of flowsheets with discontinuous investment costs. We consider as a first step the case when the investment cost functions are defined over multiple regions for the design variables. The case of multiple regions for sizes, pressures, and temperatures will be considered in a future paper. The major modeling challenge consists in effectively handling discontinuities that arise * To whom all correspondence should be addressed. Phone: 412-268-2230. Fax: 412-268-7139. E-mail: ig0c@ andrew.cmu.edu.
in the selection of the appropriate cost function depending on which region the design variable lies. Unless these discrete decisions are modeled explicitly, current continuous optimization techniques can produce suboptimal solutions or fail to converge to a feasible solution. In this paper we show that the discontinuous investment cost functions described above can be modeled as a set of disjunctions to rigorously handle the discontinuities. The disjunctive model allows conditional modeling of cost functions (e.g. if a design variable lies in a certain range, the corresponding cost expression for that range is enforced). Furthermore, it is shown that the cost modeling of the selection of parallel units when there is a constraint on the maximum available size for that equipment can be accomplished with a similar disjunctive formulation. Conventional solution methods are examined first which include a direct NLP approach, the use of smoothing functions, and the use of “big-M” constraints. Basic elements of disjunctive programming are presented next and used as a basis to derive the convex hull formulation of the disjunctions, and a disjunctive branch and bound method. Examples are presented on heat exchanger network design, process flowsheet design, and petrochemical complex synthesis, and their solutions are compared with the methods described in this paper. 2. Problem Statement Consider the optimization of a process flowsheet problem in which the investment cost functions are defined over a set of intervals D of the design variable x as shown in Figure 1. The cost in each region i ∈ D consists of two parts: fixed cost, βi, and a variable cost, Rixηi, as given by
c ) Rixηi + βi
(1)
B xi-1 e x e xBi
where x is the equipment size, the xiB define the upper
+
+
2612 Ind. Eng. Chem. Res., Vol. 35, No. 8, 1996
derivatives. The discontinuous cost functions in (2) can be modeled with the use of smoothing approximations, N-1 η1
c ) R1x + β1 +
fi((x - xBi ), µi)(Ri+1xη ∑ i)1
i+1
+
βi+1 - Rixηi + βi) (4) where Figure 1. Discontinuous investment cost function.
N ) |D|
bounds in each cost region i, Ri, βi, and ηi are cost parameters with 0 < ηi e 1. The investment cost in (1) can be expressed over all the regions i ∈ D with the following disjunction:
∨
i∈D
[
c ) Rixηi + βi B xi-1 e x e xBi
]
(2)
where ∨ is the OR operator that selects exactly one region and cost function for size x. The problem considered in this paper is the optimization of a fixed flowsheet or a superstructure given that the investment cost for each unit is given as in (2).
Although the discontinuities on the systems of equations can be approximated with this approach, it introduces a number of complications to the system. First, in the case of linear cost functions (ηi ) 1) the use of smoothing functions introduces nonlinearities. Second, the resulting function in (4) is highly nonconvex and may give rise to suboptimal solutions. Third, the accuracy and performance of the solution algorithm are very sensitive to the values of the smoothing parameter, µ, in (3). 3.3. Big-M Formulation. The big-M formulation is the simplest alternative to rigorously model disjunctions within an MINLP approach. Defining a binary variable, yi, for each interval i, and a sufficiently large parameter Ui, the disjunctions in (2) can be modeled as
c g Rixηi + βi - Ui(1 - yi)
3. Conventional Approaches Two simple approaches for handling the discontinuous investment cost (2) in optimization problems are its direct implementation within a nonlinear programming (NLP) algorithm and its formulation as “big-M” constraints within a mixed-integer nonlinear programming (MINLP) framework. 3.1. Direct NLP Approach. In this case the discontinuous cost functions are defined implicitly by modeling them as embedded IF statements. These statements are used to select the appropriate function depending on the value of the variable x in (2) within an NLP solution technique. While this approach is very simple, it has the drawback of introducing discontinuities and, therefore, potentially producing failures in gradient-based NLP optimization algorithms: first, because of the assumption that the functions involved in the model are continuous and differentiable over the entire region of the continuous variables and, second, because these methods might not select the appropriate cost region. 3.2. Smooth Approximation. Another NLP approach is the use of smoothing functions at the points of discontinuity. For instance, the sigmoid function (Chen and Mangasarian, 1995)
f(z,µ) )
1 1 + e-µz
(3)
can be used to approximate the discontinuous function.
g(z) )
{
1 0
zg0 z Z* or the problem is infeasible, prune that node, and go to step 8. Step 8: If there is no open node left, the optimal solution is found. Otherwise go to step 2.
It should be noted that the use of (21) is preferred because it generally yields tighter bounds. However, as shown with the following example, both (20) and (21) may be equivalent provided the values of Ui are selected as in (6). Consider the cost function
[
] [
c ) 2x + 10 c ) x + 25 ∨ 0 e x e 10 10 e x e 20
]
(22)
which is plotted in Figure 7. From (6) it follows that the best values of Ui in (3) are U1 ) 5 and U2 ) 15:
c g 2x + 10 - 5(1 - y1) c g x + 25 - 15(1 - y2) y1 + y2 ) 1
(23)
The surrogate relaxation in (20) yields the following for the constraints in (23):
c g 1.75x + 10
(24)
The tightest linear underestimator in (21) is given by the cost function that satisfies c ) 10 at x ) 0 and c ) 45 at x ) 20,
cLU g 1.75x + 10
(25)
which is identical to (24). Thus, in this particular example (20) and (21) yield identical relaxations. Another interesting question is how the relaxations compare between the convex hull formulation (19) and the linear underestimator (21). This is given for linear disjunctive costs by the following theorem. Theorem 3: Given the linear disjunctive cost function B B in (2), with Rixi-1 + βi > Ri-1xi-1 + βi-1 and Ri-1 g Ri, the linear underestimator given by the convex envelope, cLU ) Rx + β, is such that cLU e cCH for x ) xˆ , where cCH is the cost given by the relaxation of convex hull of disjunction in (19). Proof: see Appendix A. An outline of a proof for Theorem 3 can be found in Appendix A. Here it is shown that if the cost function is concave in the sense that the linear underestimator coincides with the cost function in the lower and upper bounds of the entire domain (see case a of Figure A), then the linear underestimator is identical to the convex hull formulation. If on the other hand the cost function is not concave in the sense that the cost function does not coincide with the linear underestimator at the extreme bounds (see case b of Figure A), then the convex hull is tighter than the linear underestimator. The implication of Theorem 3 is that from a conceptual point of view the convex hull formulation offers a theoretically superior approach to the use of linear underestimators. Note, however, that the linear underestimator is better suited for a branch and bound
+
+
Ind. Eng. Chem. Res., Vol. 35, No. 8, 1996 2617
Figure 9. Discontinuous investment cost function with maximum sizes. Figure 8. Disjunctive branch and bound tree for example 1.
approach on the disjunctions, while the convex hull representation is better suited for an MINLP model. 5.3. Example 1 Revisited. The convex hull formulation in (19) was applied to the cost function in (7) of example 1. The resulting MINLP involves 9 binary variables, 42 continuous variables, and 48 constraints as seen in Table 2. The problem was solved with DICOPT++ in 2.52 s CPU time converging to the global optimum of $117 781 yr-1. The major reason for improved performance and convergence to the correct solution is the very tight lower bound of the NLP relaxation ($115 553 vs $117 781 yr-1). Also note that this lower bound is significantly better than the use of the big-M constraints with the tightest Ui values ($115 553 yr-1 vs $102 514 yr-1). When the disjunctive branch and bound algorithm is applied for this problem the optimal solution of $117 781 yr-1 is found in 8 nodes as shown in Figure 8. At every node the deviations of the reported costs from the actual costs are compared for each heat exchanger at each interval, and the maximum deviation among the three heat exchanger costs is selected for branching. For example, at the root node the deviations in the cost for heat exchangers 1, 2, and 3 are 718.2, 2855.4, and 2439.2 with the areas 23.884, 20.509, and 7.917, respectively. The information on the area of each heat exchanger (corresponding to xr in the algorithm) and the deviation of the cost from the actual cost (c - cr) are also reported in Figure 8 for each node. The deviation in the cost for heat exchanger 2 is the maximum among the deviations; therefore the actual cost function at interval 2 is enforced at node 1, and the linear underestimator for intervals 1 and 3 in heat exchanger 2 is included at node 2. The first node has an objective function value of $115 124 yr-1, and two child nodes are generated by branching on the first term of the disjunction for heat exchanger 3. The second node corresponds to step 6 of the algorithm. The subproblem reports the area for heat exchanger 2 in the second region of the cost function which was enforced in node 1. Therefore, two child nodes are generated by enforcing the first and the third terms of the disjunctive cost function for heat exchanger 2, respectively. The search continues until the optimal solution is found as described in the algorithm. In the branch and bound tree two of the nodes yield infeasible solutions indicating that the heat exchanger areas restricted to the intervals do not satisfy the heating/cooling requirements. Two nodes, 4 and 8, are pruned after locating the optimal solution at node 7. Note that although the relaxation of the disjunctive branch and bound is somewhat looser ($112 329 yr-1 vs $115 553 yr-1), the method was slightly faster than the MINLP with convex hull formulation.
It is interesting to note that the solutions with the direct NLP approach and the use of smoothing functions are the only cases that converged to a Kuhn-Tucker point finding a solution with objective function values of $141 608 yr-1 (20.23% more expensive than the solution found with the other three methods) and $130 522 yr-1 (10.81% more expensive), respectively. As shown in Figure 5, since the cost function is discontinuous, the NLP algorithms (µ ) 0.001 and µ ) 100 in the smooth approximation) were trapped in interval 3 for heat exchanger 1 and found an area of 30.263 m2. However, the other proposed algorithms could examine other intervals and found the area of heat exchanger 1 as 25 m2 which is the upper bound for interval 2. Therefore the area of heat exchangers 2 and 3 are slightly increased to satisfy the heating/cooling requirements. 6. Design of Process Networks with Restricted Unit Sizes In the optimization of flowsheets, or more generally process networks, it is often the case that a maximum size Us be specified for a given unit. If the actual size required is larger than this maximum size, it is clear that two or more parallel units may have to be installed. In this section we describe a disjunctive programming approach that can handle restricted unit sizes by a simple extension of the disjunctive cost model in (2) and which does not require duplication of equation models for each potential parallel unit. The discontinuous cost function for the case of maximum sizes for units is illustrated in Figure 9. In the presence of a maximum size and the requirement of having several parallel equipments for the same unit the cost function is modeled with the following disjunction:
[
c ) Rixηj i + βi + (j - 1)(RNUηs N + βN) ∨ ∨ x ) xj + (j - 1)Us j∈J i∈D B xi-1j e x e xBij
]
(26)
where N ) |D|, J ) {j} potential set of identified parallel units, Us is the maximum size, x is the total size, and xj is the size of unit j, assuming others i e j - 1 have size Us. The disjunctive model (26) will design the process network considering the maximum available size for each unit type. When a processing unit requires a larger size than the maximum size, then the optimization will require one or more equipment of the same kind in parallel. The optimization will be such that the units in parallel will be used one by one (i.e., in the case
+
+
2618 Ind. Eng. Chem. Res., Vol. 35, No. 8, 1996 Table 3. Data for Example 2 stream
FCp (kW/K)
Tin (K)
Tout (K)
cost ($/kW yr)
hot 1 hot 2 hot 3 cold 1 cold 2 cold 3 cold 4 cooling water steam
30.0 13.5 20.0 31.0 5.0 28.0 11.0
626 620 528 525 405 353 313 293 650
586 519 353 613 576 386 545 308 650
20 80
heat exchanger
overall heat transfer coefficient (kW/m2 K)
1 2 3 4 5 6 7 8
1.5 0.2 0.06 1.6 0.04 0.3 0.6 1.7
interval (m2)
investment cost ($/yr)
0 < A e 20 20 e A e 50 50 e A e 100
670A0.83 + 2000 640A0.83 + 8000 600A0.83 + 16000
Figure 10. Network structure for example 2. Table 4. Heat Exchanger Areas for Example 2 area (m2)
of more than one equipment in parallel the first equipment will have the maximum size and the others will be used depending on the number of parallel equipments). The disjunction in (26) can be treated with the methods presented in the previous section. From Theorem 1 it follows that the convex hull representation of disjunction (26) is represented as follows:
x)
∑ ∑xji j∈J i∈D
xji )
∑zij
i∈D
c)
∑ ∑cji j∈J i∈D
cji ) Rizηiji + [βi + (j - 1)(RNUηs N + βN)]yij
(27)
xji ) zij + (j - 1)Usyij B xi-1j yij e xji e xBij yij
∑ ∑ j∈J i∈D
yij ) 1
yij ) 0, 1
i ∈ D, j ∈ J
It is interesting to note that variables x and c are disaggregated into xji and cji for each equipment in parallel j ∈ J and size interval i ∈ D. The design of process networks with restricted maximum sizes is illustrated with the following heat exchanger network design problem. 6.1. Example 2. Consider the design of a heat exchanger network in which the data are given in Table 3. The network has three hot and four cold streams, and cooling water and high-pressure steam at 650 K are available as utilities.
heat exchanger
no restriction on size
maximum size 100 m2
1 2 3 4 5 6 7 8
3.809 46.747 231.653 26.382 76.740 55.393 51.898 7.462
3.809 46.747 100.000 23.560 117.786 43.077 56.921 9.282
The network consists of eight heat exchangers as shown in Figure 10. In order to show the application of the proposed algorithm for the design of process networks with restricted unit sizes, first we consider the case when there is no restriction on the unit sizes and then compare the results for the case with the restricted unit sizes. When this problem is solved without a limit on the sizes of the heat exchanger areas, the optimal solution is given in Table 4. Assume that there is a limit on the maximum size of a heat exchanger: maximum area of 100 m2. The sizes are expected to change since heat exchanger 3 has an area of 231.653 m2 which is above this limit. The proposed model will detect the required number of parallel units when the size of any unit is above the maximum available size. From (26) the disjunction for the case when there is a limit of 100 m2 for the area of the heat exchangers is given by
[
][ ][ ][ ][ [
]
0.83 0.83 c ) 670A11 + 2000 c ) 640A12 + 8000 ∨ A ) A12 ∨ A ) A11 0 < A e 20 20 e A e 50 0.83 0.83 c ) 600A13 + 16000 c ) 670A21 + 2000 ∨ A ) A21 + 100 ∨ A ) A13 50 e A e 100 100 < A e 120 0.83 0.83 c ) 640A22 + 8000 c ) 600A23 + 16000 ∨ A ) A23 + 100 ∨ A ) A22 + 100 120 e A e 150 150 e A e 200 0.83 0.83 c ) 670A31 + 2000 c ) 640A32 + 8000 ∨ A ) A32 + 200 ∨ A ) A31 + 200 200 < A e 220 220 e A e 250 0.83 c ) 600A33 + 16000 ∨ (28) A ) A33 + 200 250 e A e 300
[ [ [
]
] ] ]
When one heat exchanger requires an area greater than 100 m2 but less than 200 m2, two heat exchangers must be placed in parallel (fourth, fifth, or sixth term of disjunction). Three parallel heat exchangers must be
+
+
Ind. Eng. Chem. Res., Vol. 35, No. 8, 1996 2619 Table 5. Summary of Results for Example 2 item
big-M method
No. of 0-1 variables No. of cont. variables No. of constraints cost ($/yr) relaxed cost ($/yr) No. of iterations CPU time (s)d
72 88 143 423 938 268 734 3a,b 21.99
convex hull disjunctive of disjunction branch & bound 72 176 239 412 390 399 467 4a 12.11
N/A 64 47 412 390 386 471 60c 20.89
a Number of major iterations with DICOPT++. b Termination when there is no improvement on NLP subproblems. c Number of nodes. d IBM RS6000/530.
Figure 12. Process flowsheet for example 3. Table 6. Problem Data in Example 3 feed
composition (%)
cost ($/kg mol)
A B D
65 30 5
0.033
stream
specification
cost ($/kg mol)
main product
g90% C e1000 tons/day
0.28
by-product
Figure 11. Heat exchanger areas for example 2.
installed for a total area larger than 200 m2 (seventh, eighth, or ninth term of the disjunction). The summary of results with the proposed algorithms is shown in Table 5. As it is seen from the results, the convex hull formulation of the disjunction finds the optimal solution in four major iterations with DICOPT++ (Viswanathan and Grossmann, 1990). A suboptimal solution is obtained in three major iterations with the big-M formulation. This is due to the heuristic termination of DICOPT++ which stopped after no improvement was found in the NLP subproblem. The disjunctive branch and bound algorithm requires the evaluation of 60 nodes before terminating the search. In this example it is clear that despite its increased size the convex hull formulation is substantially more efficient than the big-M and disjunctive branch and bound methods. The results given in Table 4 show that in the optimal solution with a limit of 100 m2 the area of the heat exchanger 5 is increased from 76.740 to 117.786 m2. Thus, it is required to have two heat exchangers for unit 5: one with a maximum available size of 100 m2 and the other with an area of 17.786 m2. For exchanger 3 on the other hand, the maximum size of 100 m2 is selected. The resulting design for the heat exchanger network is shown in Figure 11. The proposed model was able to detect the requirement of parallel equipments for heat exchanger 5. 6.2. Example 3. In this example it will be shown that the disjunctive models and solution techniques are useful in realistic process design problems. It is desired to produce chemical C with a minimum purity of 90% and an amount of 1000 tons/day. The flowsheet contains the following units; two-stage compression with intercooling in the feed stream, a plug flow reactor, four coolers, three heaters, a single-stage recycle stream compressor, and a flash separator as shown in Figure 12. The purge stream has a market value as
0.021
utility
cost
electricity steam water
$0.03/kW h $8.0/106 kJ $0.7/106 kJ
fuel. The key optimization variables in the reactor are the operating pressure and the per pass conversion. It is expected that the reactor operates at high pressures; therefore, reactor feed stream and the recycle stream are compressed. The problem data for the example are given in Table 6. In order to satisfy operational safety and construction material specifications, the operating pressure in the reactors is restricted to at most 15 MPa and the outlet temperature from the reactor is limited to at most 873 K. The temperature in the flash unit is specified between 300 and 500 K. These restrictions are also included in the model. Each unit has three different cost functions depending on the sizes. Heat exchangers are limited to a size of 100 m2. In other words this illustrative example combines the proposed approaches for handling discontinuous investment costs and the presence of a limit on the maximum available size of a particular unit. The investment cost is annualized by considering a project life of 3 years. The cost function for each unit is given in Table 7. The cost of compressors is a function of the compressor load, W, where unit load is 1000 kW. The investment cost for reactor is composed of two parts: reactor vessel cost and the catalyst cost. The vessel cost is a function of volume of the reactor and also includes fixed pressure and construction material factors as shown in Table 7. The catalyst cost is determined as $3616 m-3 with a 20% replacement every year. Investment cost for heat exchangers is a function of actual area, pressure factor, and construction material for shell/tube. It is considered that maximum area of a single heat exchangers is 100 m2, and if more heat exchange area is required more heat exchangers in parallel will be used. The flash separator cost is a function of feed flowrate and the residence time, τf, in the separator. Note that in order to update the costs Marshall and Swift equipment cost index for the fourth quarter of 1994 (M&S ) 1004.4) is used in the model. The optimal solution with the convex hull formulation which involves 51 0-1 variables, 307 continuous variables, and 304 constraints is found in three major
+
+
2620 Ind. Eng. Chem. Res., Vol. 35, No. 8, 1996 Table 7. Investment Cost Functions for the Units unit
interval
investment cost ($1000 yr-1)a
compressor W ) (1000 kW)
0.0 e W e 1.5 1.5 e W e 2.5 2.5 e W e 4.0 0 e V e 30 30 e V e 60 60 e V e 100 0eAe1 1eAe2 2eAe3 0 e F e 10 10 e F e 20 20 e F e 30
M&S/280(50.5W0.706 + 1) M&S/280(55.9W0.706 + 5) M&S/280(60.5W0.706 + 9) M&S/280FpFm(0.125V0.738 + 1) M&S/280FpFm(0.150V0.738 + 1.5) M&S/280FpFm(0.175V0.738 + 2) M&S/280FpFe[(3.17A0.641) + 0.75] M&S/280FpFe[(3.17A0.641) + 4.67] M&S/280FpFe[(3.17A0.641) + 8.59] M&S/280FpFm[(0.240τf0.738F0.738) + 1] M&S/280FpFm[(0.289τf0.738F0.738) + 1.5] M&S/280FpFm[(0.337τf0.738F0.738) + 2]
reactor vessel V ) (m3) heat exchanger A ) (100 m2) flash F ) (kg mol/s)
a F ) 1.20 pressure factor; F ) 2.25 construction material factor (SS); F ) 2.81 construction material factor (CS/SS); τ ) 5 min p m e f residence time.
Table 8. Summary of Results for Example 3 item
big-M method
No. of 0-1 variables No. of cont. variables No. of constraints profit ($/yr) relaxed profit ($/yr) No. of iterations CPU time (s)b
51 273 273 1883.86 2100.81 4a 16.83
convex hull of disjunction 51 307 304 1883.86 1885.77 3a 8.49
a
Number of major iterations with DICOPT++. b IBM RS6000/ 530. Table 9. Optimal Sizes and Investment Costs for the Units unit
size
2: compressor 1 4: compressor 2 7: reactor
1739 kW 1739 kW 38.3 m3
9: cooler 10: flash 13: compressor 2
145.3 m2 17.9 kg mol/s 1731 kW
investment cost ($/yr) 314 331 314 331 vessel: 63 781 catalyst: 138 637 79 611 91 852 313 287
iterations with DICOPT++ (Viswanathan and Grossmann, 1990). The problem requires 8.49 s of CPU time on an IBM RS6000/530 workstation. The big-M formulation which involves 51 0-1 variables, 273 continuous variables, and 273 constraints finds the same optimal solution in four major iterations requiring 16.83 CPU s on the same workstation. The progress of the iterations with these methods are shown in Table 8. When the discontinuities in the investment costs are approximated with smoothing functions due to large scale of this problem the convergence to any solution could not be achieved with different values of smoothing parameters (µ ) 0.0001, µ ) 0.01, µ ) 1, µ ) 5, µ ) 100). The optimal sizes of the units in the flowsheet are shown in Table 9. The two-stage feed compressor has a load of 1739 kW on each stage and costs a total of $628 662 yr-1. The optimal volume of the reactor is found as 38.3 m3. The stainless steel reactor vessel costs $63 781 yr-1, and the catalyst cost is $138 637 yr-1. The flash separator cooler has an area of 145.3 m2. Since the maximum available size is 100 m2, two heat exchangers in parallel are required. The flash separator has a feed flowrate of 17.9 kg mol/s and costs $91 852 yr-1. The recycle stream compressor has a load of 1731 kW and costs $313 287 yr-1. The optimal design of the flowsheet makes a profit of $1 883 860 yr-1. 7. Synthesis of Process Networks with Discontinuous Investment Costs Using the generalized disjunctive programming formulation for process networks by Raman and Gross-
mann (1994) as a basis, the synthesis of linear process networks with linear discontinuous investment costs as in (2), and considering a set of K units, can be expressed by the following model:
min Z ) γTx + s.t.
[[
Yk Akxk e bk ck ) Rkixk + βki ∨ B B i∈Dk xi-1 e xk e xi
]
]
∑k ck
Ex e e
[ ]
¬ Yk ∨ xk ) 0 ck ) 0
k∈K
(29)
Ω(Yk) ) true x ∈ X, Y ∈ {true, false}
where xk is the subset of variables x belonging to unit k, and ck is its corresponding fixed charge. The disjunctive model (29) involves two types of decisions: (1) whether a unit exists in the optimal structure and (2) the investment cost of the existing units depending on the different size intervals. From Theorem 2 it follows that the convex hull for every disjunction k ∈ K in (29) is given by the following constraints:
xk ) x1k + x2k x1k )
∑i xki1
ck ) c1k + c2k c1k )
∑i cki1
Akx1k e bky1k 1 1 cki ) Rkixki + βkizki B 1 xi-1 zki e xki e xBi zki
(30)
x2k ) 0y2k c2k ) 0y2k yk ) y1k + y2k yk )
∑i zki
y1k, y2k ) 0, 1, zki ) 0, 1
i ∈ Dk
+
+
Ind. Eng. Chem. Res., Vol. 35, No. 8, 1996 2621
Figure 13. Process network superstructure in example 4.
From (30) the convex hull formulation reduces to
xk )
∑i xki
ck )
∑i cki
Akxk e bkyk
(31)
xLi yk e xk e xU i yk cki ) Rkixki + βkizki B xi-1 zki e xki e xBi zki
yk )
∑i zki
yk ) 0, 1, zki ) 0, 1
i ∈ Dk
The incorporation of constraints (31) in problem (29) is illustrated with the following petrochemical complex synthesis problem. 7.1. Example 4. Consider the superstructure of a petrochemical complex in Figure 13 consisting of 38 different processes and 28 chemicals; 17 raw materials, 9 final products, and 2 intermediate products. The network superstructure is described in Sahinidis et al. (1989). Each process in the superstructure has three different operating capacities: low, medium, and high depending on the amount of the main product from each process. The investment cost function of the processes are discontinuous, since different investment cost expressions are valid for each design capacity. The big-M model for this problem involves 190 binary variables, 1175 continuous variables, and 1461 constraints as shown in Table 10. The model with the convex hull formulation of disjunctions is larger than the big-M model consisting of 190 binary variables, 1403 continuous variables, and 1689 constraints. The relax-
+
+
2622 Ind. Eng. Chem. Res., Vol. 35, No. 8, 1996 Table 10. Summary of Solutions for Example 4 item
big-M method
No. of 0-1 variables No. of cont. variables No. of constraints profit ($/yr) relaxed profit ($/yr) No. of iterations No. of nodes CPU time (s)a
190 1175 1461 no solution 409 388 900 >1 000 000 >150 000 N/A
convex hull of disjunction 190 1403 1689 127 574 400 220 663 200 2779 1182 161.32
ation gap for the convex hull formulation is 72.97%, and that for the big-M formulation is 220.86%, which is consistent with Theorem 2 in which it was proved that the relaxation gap with the convex hull formulation is tighter than that with the big-M model. This problem shows that the difference in the relaxation gap plays an important role in finding the optimal solution in a reasonable amount of time. As it is seen in Table 10, the optimal solution with the convex hull formulation is found in about 3 min requiring the examination of 1182 nodes and 2779 simplex iterations with GAMS/ OSL (Brooke et al., 1992). The big-M formulation failed to even find a feasible solution after one million iterations! The optimal flowsheet configuration that is obtained includes 11 processes (processes 4, 5, 6, 8, 14, 16, 17, 28, 29, 32) with 6 raw materials and 12 products and a predicted profit of $127 574 400 yr-1.
8. Conclusions The treatment of discontinuous cost functions that are defined over several regions and with fixed charges has been addressed for the optimization and synthesis of process networks. It was shown that direct NLP techniques and the smoothing functions used to solve the optimization models with discontinuous investment cost functions may converge to suboptimal solutions, or may fail to converge to a feasible solution. It was shown that big-M models may produce suboptimal solutions in some cases, or may fail to converge to a solution in a reasonable amount of time for large scale problems due to their loose relaxation. Two alternative modeling techniques based on disjunctive programming have been proposed: convex hull formulation and disjunctive branch and bound, based mainly on linear underestimators. The big-M model requires the selection of a big-M parameter and the use of a binary variable for each size interval of a unit. The convex hull formulation eliminates the artificial big-M parameters and is a natural representation of discontinuities in realistic cost functions. The linear underestimators are relaxations of the discontinuous investment cost functions and can be used within a disjunctive branch and bound algorithm. Several theoretical results have been presented, the major one being that the convex hull formulation yields the tightest relaxation amongst the alternative modeling techniques examined in this paper. The numerical performance of the proposed solution techniques was examined with examples on heat exchanger networks, process flowsheets, and process networks. The disjunctive branch and bound algorithm performs reasonably well for small problems in which the solution of subproblems at each node of the branch
Figure A. Cost functions in Appendix A.
and bound tree is not expensive. Since the solution of subproblems for process network optimization problems is often expensive, the convex hull formulation tends to be computationally efficient and produces reliable and robust results as was shown in the examples of this paper.
Acknowledgment The authors would like to acknowledge financial support from the Engineering Design Research Center at Carnegie Mellon University and from the Eastman Chemical Company.
Appendix A: Relaxations of Convex Hull and Linear Underestimator Theorem 3: Given the linear disjunctive cost function B B + βi > Ri-1xi-1 + βi-1 and Ri-1 g Ri, in (2), with Rixi-1 the linear underestimator given by the convex envelope, cLU ) R jx + β h , is such that cLU e cCH for x ) xˆ , where cCH is the cost given by the relaxation of convex hull of disjunction in (19). Outline Proof: Rather that providing a general proof we will restrict ourselves to the case of two regions to show the above theorem. The proof can be generalized to the case of multiple intervals although the detailed proof is very cumbersome. Proof for case of two intervals: consider
0 e x e xB1 , c ) R1x1 + β1 xB1 e x e xB2 , c ) R2x2 + β2
(A1)
Case a: Concave Cost Function. First assume a concave cost function for which the two following conditions hold (see case a of Figure A):
I. II.
R2xB1 + β2 > R1xB1 + β1 R1 > R2, and R1 >
[
]
R2xB2 + β2 - β1 xB2
The linear underestimator for this case is then given by
cLU )
[
]
R2xB2 + β2 - β1 xB2
x + β1
(A2)
From (19) the convex hull representation of the disjunc-
+
+
Ind. Eng. Chem. Res., Vol. 35, No. 8, 1996 2623
tive cost function is
cCH ) c1 + c2
(A3)
x ) x1 + x2
(A4)
c1 ) R1x1 + β1y1
(A5)
c2 ) R2x2 + β2y2
(A6)
0 e x1 e xB1 y1
(A7)
xB1 y2 e x2 e xB2 y2
(A8)
y1 + y2 ) 1
(A9)
Rearranging (A12) by using (A9), we can determine the expressions for y1 and y2:
y1 )
y2 )
cCH )
[
]
R2xB2 + β2 - β1 xB2
x + β1
(A10)
Since (A2) and (A10) are equivalent to each other, the cost given by the linear underestimators is equal to the one given by the convex hull representation when assumptions I and II are true. Case b: Nonconcave Cost Function. In order to prove that the relaxation with the convex hull representation is tighter than linear underestimators, consider that assumption II of case a is not true, that is R1 g R2 (consider R1 ) R2, see case b of Figure A). From assumption I we can deduce that β2 has to be greater than β1. The linear underestimator in this case is given by
cLU ) R1x + β1
(A11)
The convex hull formulation of the disjunctive cost function is given with (A3)-(A9). If the size of unit lies in the first size interval (i.e., 0 e x e xB1 ), then the linear underestimator represents the cost function exactly. The convex hull representation becomes cCH ) R1x + β1 by setting y1 ) 1 and y2 ) 0, x1 ) x and x2 ) 0. Therefore, in the first size interval the costs predicted with two approaches are equal. If the size of unit lies in the second size interval (i.e., xB1 e x e xB2 ), we want to find the largest feasible y1 to make xB1 y1 e x in (A7) to investigate the worst cost predicted with the convex hull formulation. Substituting active constraints from (A7) and (A8) into (A4),
xB1 y1 + xB2 y2 ) x
(A12)
(A13)
xB2 - xB1 x - xB1
(A14)
xB2 - xB1
Substituting (A5), (A6), (A13), and (A14) into (A3) and using (A9) and (A12), the relaxation of the cost function with convex hull formulation is given by
[ ] [ ]
cCH ) R1x + β1
xB2 y2,
Let x1 ) 0 and x2 ) then y2 becomes equal to x2/xB2 . From (A9), y1 is equal to 1 - x2/xB2 , and since x1 ) 0, and using (A4) x ) x2. Substituting (A5), (A6) and the above expressions for y1, y2, and x into (A3), and rearranging the cost function becomes
xB2 - x
xB2 - x
xB2 - xB1
+ β2
x - xB1
xB2 - xB1
(A15)
Since 0 e y1 e 1 and 0 e y2 e 1, the cost given by (A15) is always greater than the cost given by (A11). Therefore, the relaxation of the convex hull formulation is stronger than the linear/convex underestimators. Literature Cited Balas, E. Disjunctive programming and a hierarchy of relaxations for discrete optimization problems. SIAM J. Alg. Disc. Meth. 1985, 6, 466-486. Biegler, L. T. Optimization strategies for complex process models. Adv. Chem. Eng. 1992, 18, 197-256. Brooke, A.; Kendrick, D.; Meeraus, A. GAMS a user’s guide, release 2.25; Boyd & Fraser Publishing Co.: Danvers, MA, 1992. Beaumont N. An algorithm for disjunctive programs. Eur. J. Operational Res. 1990, 48, 362-271. Chen, C.; Mangasarian, O. L. A class of smoothing functions for nonlinear and mixed complementarity problems. Mathematical Programming Technical Report 94-11; University of Wisconsin: Madison, WI, 1995. Guthrie, K. M. Process plant estimating, evaluation, and control; Craftsman Book Co.: Los Angeles, 1974. Raman, R.; Grossmann, I. E. Modeling and computational techniques for logic based integer programming. Comput. Chem. Eng. 1994, 18, 563-578. 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-1063. Turkay, M.; Grossmann, I. E. Logic-based MINLP algorithms for the optimal synthesis of process networks. Comput. Chem. Eng. 1996, 20, 959-978. Viswanathan, J.; Grossmann, I. E. A combined penalty function and outer-approximation method for MINLP optimization. Comput. Chem. Eng. 1990, 14, 769-782.
Received for review February 20, 1996 Revised manuscript received April 18, 1996 Accepted April 29, 1996X IE9600856
X Abstract published in Advance ACS Abstracts, June 15, 1996.
