Multiobjective Optimization for Design and Operation of the Chilling

May 20, 2010 - The chilling train is a very important section for ethylene plants, as its design ... (1) Because of the rapid growth of the worldwide ...
0 downloads 0 Views 4MB Size
5786

Ind. Eng. Chem. Res. 2010, 49, 5786–5799

Multiobjective Optimization for Design and Operation of the Chilling Train System in Ethylene Plants Jian Zhang, Yanqin Wen, and Qiang Xu* Dan F. Smith Department of Chemical Engineering, Lamar UniVersity, Beaumont, Texas 77710

The chilling train is a very important section for ethylene plants, as its design and operation significantly influence energy consumption and product loss rates. In this paper, a systematic methodology for the optimal design and operation of a chilling train system has been presented. A multiobjective mixed-integer nonlinear programming (MINLP) model has been developed, which aims at minimizing the ethylene loss rate and the exergy-accounted energy consumption, meanwhile maximizing the hydrogen recovery rate. Also presented are the technologies for obtaining reliable thermodynamic properties based on rigorous simulation and the methods for removing redundant variables for MINLP model simplification. On the basis of the methodology, the optimal design and operating conditions of a chilling train system can be simultaneously achieved. A real case study has demonstrated the efficacy of the developed methodology. The results show that the optimized chilling train system has better economic performance than the current system. 1. Introduction Ethylene is the most important organic product with the largest bulk productivity in the petrochemical industry. Its derivatives, like ethylene oxide, vinyl acetate, linear alcohols, ethylene dichloride, and high/low density polyethylene, are extremely important in daily life.1 Because of the rapid growth of the worldwide economy and population, global ethylene production in 2007 reached about 115 million metric tons and is expected to increase by 4.4% per year from 2007 to 2012.2 The chilling train is an important part of ethylene plants. Within this section, the charge gas from the upper stream process is refrigerated to a very low temperature (lower than -160 °C) to separate hydrogen and methane. Meanwhile, the left charge gas is liquidized for cryogenic separation. The chilling train system consumes over 40% of the refrigerating capacity of an ethylene plant. Different operational conditions with different temperature and pressure settings will affect both product separation efficiencies and energy consumption. Therefore, the operational optimization for the chilling train system has great economic potential and needs sufficient study. Figure 1 gives a typical chilling train flowsheet. The system input is the charge gas stream from the Depropanizer Tower, which normally includes 13.95 mol % of H2, 23.83 mol % of CH4, 43.43 mol % of C2H4, 5.96 mol % of C2H6, 10.62 mol % of C3H6, and 2.21 mol % of C3H8. Its input temperature and pressure are -19.6 °C and 1.75 MPa, respectively. The charge gas is refrigerated to about -70 °C and separated at the firststage flash drum (FLS1). The liquid stream is sent to the Demethanizer Tower. The vapor stream is refrigerated to about -105 °C and separated at the second-stage flash drum (FLS2). The liquid stream is also sent to the Demethanizer Tower, while the vapor stream is heated to the ambient temperature and is then compressed to 3.2 MPa. The compressed stream is refrigerated to -120 °C and separated by the third-stage flash drum (FLS3). Again, the liquid stream is directed to the Demethanizer Tower. The vapor stream, which mainly consists of H2 and CH4, is refrigerated to -140 °C and separated at the fourth-stage flash drum (FLS4). Since FLS4 works at such a low temperature that no other refrigerants can be used to further * To whom correspondence should be addressed. Phone: 409-8807818. Fax: 409-880-2197. E-mail: [email protected].

refrigerate its inlet stream, the liquid outlet stream of FLS4 itself has to be utilized as a refrigerant. Thus, the liquid stream containing mainly methane is sent to a throttling valve to reduce its pressure and temperature and is then reused to refrigerate the inlet stream from FLS3, as well as other upper streams throughout chilling train cold boxes. The vapor stream of FLS4 is further refrigerated to -160 °C and sent to the fifth-stage flash drum (FLS5). In FLS5, hydrogen and methane are separated as the vapor and liquid streams, both of which are reused to refrigerate the inlet stream of FLS5 and the upper streams throughout the chilling train. Note that the liquid stream is sent to a throttling valve to reduce its pressure and temperature before its reuse. Finally, the liquid stream of FLS5 containing a high purity of methane is heated to the ambient temperature after reversibly passing through the chilling train. It will be sent to the fuel gas system combined with the other two methane streams: one is from the Demethanizer Tower and the other is from the aforementioned liquid stream of FLS4. The vapor stream of FLS5 containing about 95% hydrogen is heated to the ambient temperature after reversibly passing through the chilling train and will be reused by other facilities. Note that there is a small amount of ethylene supposedly left in the inlet stream to FLS4. It is accounted for as the product loss because it cannot be recovered in the subsequent process. Also note that the hydrogen left in the liquid stream from all the flash drums cannot be recovered, which is also accounted for as a material loss. Figure 2 gives a simplified flowsheet for the chilling train system depicted in Figure 1. Compared with Figure 1, when a stream passes through a serial of heat exchangers (including the cold box) before reaching a flash drum or the compressor, these heat exchangers can be aggregated together and considered as one logic heat exchanger for this stream. With this approximation, the complexity of cold box calculations can be decoupled so that the energy consumption for each stream and the overall energy transfer for the chilling train system are easier to calculate. The reason that one can have this approximation is that each cold box has one refrigerant stream from outside of the chilling train system that can be freely manipulated to help accomplish the energy balance within the cold box. Certainly, the performance of such an approximation will be

10.1021/ie100455g  2010 American Chemical Society Published on Web 05/20/2010

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

5787

Figure 1. Typical chilling train system.

validated through rigorous simulations, after the results of the optimal design and operation based on the simplified chilling train system are obtained. From Figure 2, it is obvious that changing the temperature or pressure of the five flash drums will change vapor fraction and the compositions of vapor and liquid streams. These changes will affect the energy consumption and product separation efficiency to influence the ethylene loss rate and hydrogen recovery rate. On the other hand, the fourth and fifth stage flash drums should operate at high pressure (about 3.2 MPa), but the input charge gas stream has only 1.75 MPa. Thus, the charge gas should be sufficiently compressed before entering the fourthstage flash drum. Theoretically, the compressor can be placed before or after each of the first three flash drums. The front position of the compressor (e.g., the compressor is placed prior to the first flash drum) will reduce the heat duties of heat exchangers but increase the compressor load and work consumption, and vice versa. Thus, the compressor placement and the operation conditions (temperature and pressure settings of flash drums) of the chilling train system should be optimized to balance the overall energy consumption, ethylene loss rate, and hydrogen recovery rate. To obtain the optimal design and operation, a first-principle-based chilling train synthesis model has to be employed. Although the chilling train has been widely used on ethylene plants, the literature on its optimal design and operational optimization is still lacking. Chiu and Newton presented a cryogenic process analysis based on the Second Law of Thermodynamics.3 Their work focused on LNG (Liquefied Nature Gas) processes. Colmenares and Seider presented a nonlinear model for the synthesis of a refrigeration system integrated with a cryogenic process.4 Hurmes et al. presented a qualitative decision support system for a chilling train.5 The model implies sign algebra to get qualitative relations on process variables. The model generation is more straightforward than previous synthesis methods. However, this method cannot guarantee the solution optimality. Dhole and Linnhoff proposed an overall design method for cryogenic processes.6 In this method, a cryogenic system was divided into three parts: separation process, heat exchange system, and refrigeration system. Using a combined method of pinch analysis and exergy analysis, the model could optimize the shaft work target and column target. Castillo and Dhole optimized the operating pressure of the ethylene cold-end process, where pinch analysis

and exergy analysis were employed in this simulation-based study.7 Wu and Zhu presented an algorithm that simultaneously integrated a refrigeration system and heat exchanger network.8 It can be used at conceptual design to determine the main flowsheet and parameters. Lee et al. presented an optimal synthesis method for a mixed-refrigerant system.9 The algorithm employed some thermodynamic analysis to gain insight into the system and nonlinear programming techniques to determine the optimal operating conditions. On the basis of the literature survey, the exact study on multiobjective optimization for design and operation of the chilling train system has not yet been reported. Since available software tools are incapable of handling simultaneous optimization for the design and operation of a chilling train system, a novel methodology needs to be developed. In this paper, a systematic methodology for chilling train system synthesis in ethylene plants has been presented. A multiobjective MINLP (mixed-integer nonlinear programming) model based on rigorous process simulation has been developed. On the basis of the methodology, the optimal design and operating conditions of a chilling train system can be simultaneously achieved. A real case study has demonstrated the efficacy of the developed methodology. 2. Methodology Framework The developed methodology framework generalizes threestage work: model development, model simplification, and solution identification and validation. Figure 3 gives the big picture. In the first stage, the model development includes system analysis, superstructure development, and synthesis model development. The main purpose of the system analysis is to investigate the current chilling train system, discover the potential opportunities to optimize process operations, and collect process data for modeling. On the basis of the system analysis, a superstructure of a chilling train is developed with the help of industrial expertise, which contains all alternative designs. The third step is to develop the synthesis model. During this step, a multiobjective MINLP model will be developed. Lots of thermodynamic property calculations will be involved when solving the developed MINLP model. To make the optimization more efficient and utilize the reliable thermodynamic data from a commercial simulator, the synthesis model needs to be appropriately simplified based on a rigorous

5788

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

Figure 2. Simplified flowsheet for the current chilling train system.

Figure 3. Diagram of the methodology framework.

simulator. Thus, in the second stage, the model simplification includes three major tasks: model data preparation, thermodynamic function estimation, and redundant variable elimination. During the model data preparation, sufficient rigorous process simulations are conducted to get enough data for major equipment. Aspen Plus is used for the process simulations.10 Since many simulation scenarios are needed, some user-defined automation scripts are developed and imbedded into Aspen Plus to generate and collect the simulation data effectively and efficiently. The second step is to approximate thermodynamic property functions through data regressions. Linear functions, piecewise linear functions, and quadratic functions are used for

such regressions. To further reduce the complexity of the developed synthesis model, to improve the solving efficiency, redundant variables and equations in the model will be systematically eliminated in the third step. In the last stage, the optimal synthesis is programmed and solved with GAMS.11 When the solution is identified, it should be validated. For this purpose, rigorous process simulation is performed based on the obtained optimized results. If the optimized results are proven to be valid, the final solution for design and operation of the chilling train system is obtained. Otherwise, troubleshooting has to be conducted in the previous two-stage work.

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

5789

Figure 4. Superstructure for the chilling train system synthesis.

3. Chilling Train Synthesis Model The chilling train synthesis model contains a design superstructure and detailed unit models for the heat exchanger, compressor, flash drum, valve, mixer, as well as multiple objective functions and process specifications, which are presented in the following sections. 3.1. Chilling Train Superstructure. Since the input charge gas in Figure 2 should be compressed from 1.75 MPa to about 3.2 MPa at the fourth and fifth stage flash drums, it can be envisioned that the compressor can be placed before or after each of the first three flash drums. On the basis of this idea and the in-depth process analysis, Figure 4 gives the superstructure for a chilling train system design. Stream S1 has two alternative flow possibilities. It will be either refrigerated to -67 ∼ -77 °C (stream S7) or compressed to 3.4-3.6 MPa and then be refrigerated to -67 ∼ -77 °C (streams S2fS5fS6). The exact refrigeration pressure and temperature will be determined through optimization. After that, the streams are sent to the firststage flash drum (FLS1B or FLS1A) for separation into the liquid and vapor streams. The liquid stream is sent to the Demethanizer Tower. The vapor streams may have one or two alternative flow routes, depending on their pressures. If the pressure of the vapor stream is higher than 3.4 MPa, the vapor stream will be refrigerated (stream S9). If the pressure of the vapor stream is lower than 2.0 MPa, the vapor stream can either be refrigerated directly (stream S18) or be compressed to 3.4-3.6 MPa and get refrigerated (stream S16). Then, the streams are sent to the second-stage flash drums (FLS2A, FLS2B, or FLS2C) for separation into liquid and vapor streams. Similarly, the liquid stream from a second-stage flash drum is sent to the Demethanizer Tower, and the vapor stream will

Table 1. Temperature Range of Flash Drums unit

temperature high limit (°C)

temperature low limit (°C)

first stage flash drum second stage flash drum third stage flash drum fourth stage flash drum fifth stage flash drum

-67 -100 -120 -135 -160

-77 -110 -130 -140 -165

be refrigerated (streams S27 and S33) or compressed before refrigeration (stream S31). Then, the stream is sent to the thirdstage flash drums (FLS3A, FLS3B, or FLS3C). The liquid streams from the third-stage flash drums are directed to the Demethanizer Tower, and the vapor streams may get refrigerated or compressed then refrigerated. If the pressure of the vapor stream is still lower than 2.0 MPa, the vapor stream should be compressed and then refrigerated (stream S42) because the subsequent process should work at the high pressure of 3.4-3.6 MPa; otherwise, the vapor stream will be directly refrigerated (streams S35 and S37). All refrigerated streams should be mixed and sent to the flash drum of FLS4. Starting from FLS4, the downstream process flowsheet is similar to Figure 2. Stream S45 is sent to the throttling valve of VA1 to reduce its pressure and temperature first and then reused to refrigerate the stream S43 through the heat exchanger of HX14. After that, stream S47 is heated to atmospheric temperature and sent to the methane fuel gas system. The vapor stream of FLS4 is refrigerated and sent to FLS5. In FLS5, hydrogen and methane are separated as the vapor and liquid streams, both of which are heated and sent out of the chilling train system. Table 1 gives the temperature ranges of all the flash drums.

5790

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

The chilling train superstructure includes several alternative designs. To rigorously model the candidate designs and identify the optimal one, binary variables yi are employed, which is 1 if the associated process stream (represented by the index k) is selected; otherwise, it is 0 (represented by the stream index l). The model equations are described below. Fk + Fl ) Fj

(1)

Fk e FMaxyi

(2)

Fl e FMax(1 - yi)

(3)

Pk ) Pjyi

(4)

Pl ) Pj(1 - yi)

(5)

Tk ) Tjyi

(6)

Tl ) Tj(1 - yi)

(7)

Mk ) yiMj

(8)

Ml ) (1 - yi)Mj

(9)

Hk ) Hjyi

(10)

Hl ) Hj(1 - yi)

(11)

Sk ) Sjyi

(12)

Sl ) Sj(1 - yi)

(13)

(15)

Mk ) Mj

(16)

Hk ) ηk(Tk, Pk, Mk)

(17)

Sk ) φk(Tk, Pk, Mk)

(18)

Qu ) Fj(Hk - Hj)

(19)

Bu0 ) Fj{[Hk - Sk(T0 + 273.15)] - [Hj - Sj(T0 + 273.15)]} (20) -Mzu e Bu0 e M(1 - zu)

where (j,k,l) ∈ {(S1,S2,S3), (S11,S13,S17), (S25,S28,S32)} in this superstructure. yi ∈{y1,y2,y3}. F is the molar flow rate of streams. P and T are the pressure and temperature of streams. FMax is the upper limit of all streams’ flow rate. Mj ) (Mj,H2,Mj,CH4,Mj,C2H4,Mj,C2H6,Mj,C3H6,Mj,C3H8) is the molar fraction vector of the jth stream of a stream, containing species of H2, CH4, C2H4, C2H6, C3H6, and C3H8. Similarly, Mk ) (Mk,H2,Mk,CH4,Mk,C2H4,MkC2H6,Mk,C3H6,Mk,C3H8). Hj is the upstream j’s specific enthalpy (enthalpy per molar). Hk and Hl are specific enthalpies of the possible down streams k and l, which are calculated by eqs 10 and 11. Sj is the upstream j’s specific entropy (entropy per molar). Sk and Sl are specific entropies of the possible down streams k and l, which are calculated by eqs 12 and 13. 3.2. Heat Exchanger Model. Through heat exchangers, process streams are heated or cooled by refrigerants or other process streams. To address the energy consumption of the heat exchanger model, the concept of exergy is employed. Exergy is a measure of the departure of the state of a system from that of the environment. It can be defined as the maximum obtainable work from the combination of the system and the environment.12 If the charge gas stream in Figure 4 is cooled by a refrigeration system, its exergy change suggests the consumption of the work (the minimum required work) provided by the compressor of the refrigeration system. For a heat exchanger u, suppose its inlet and outlet streams are j and k. The heat exchanger model equations are generalized as below Fk ) Fj

Pk ) Pj - ∆Pu

(14)

(

Bu ) Bu0 euzu +

1 - zu eu

)

(21) (22)

where (j,k) ∈ {(S2,S4),(S5,S6), ..., (S50,S52)} according to the superstructure in Figure 4. u ∈ HX ) {HX1,HX2, ..., HX18}. ∆Pu is the fixed pressure drop of the heat exchanger u. Hk is the specific enthalpy (enthalpy per molar) of the kth stream. Sk is the specific entropy (entropy per molar) of the kth stream. Both Hk and Sk are functions of the stream’s temperature, pressure, and compositions (represented by ηk (•) and φk (•)). Qu is the heat duty of the heat exchanger. Bu0 is the theoretical exergy consumption of the stream in the heat exchanger. The positive value of Bu0 means the stream consumes exergy. The negative value of Bu0 means the stream generates exergy. eu is the exergy use efficiency. z is a preset binary indicating whether the flow is releasing exergy or consuming exergy. T0 is the reference temperature, which is 30 °C in this paper. Note that Tk is the outlet temperature of the heat exchanger. It will be optimized by the chilling train synthesis model. Once it has been determined, the heat duty (Qu) and the actual exergy consumption/generation (Bu) of the u-th heat exchanger can be identified by eqs 19 and 22, respectively. A heat exchanger always has exergy loss because of the temperature difference between the hot and cold streams. For instance, the heat exchanger CD3 in Figure 1 has two streams: the hot stream, changing from -19.6 to -29.8 °C, and the cold stream, changing from -34.0 to -42.0 °C. The average temperature difference is 5.6 °C. According to the calculation, the exergy generated by the cold stream is 3.80 GJ/h, while the exergy consumed by the hot stream is 3.53 GJ/h. The exergy use efficiency is eu ) (3.80 - 3.53)/3.80 ) 92.7%. Similarly, Table 2 gives the exergy efficiencies of all heat exchangers in the current chilling train system (see Figure 1). It shows that most of the exergy use efficiencies are larger than 86%, and the average exergy efficiency is 89.5%. Thus, the exergy use efficiency can be roughly estimated as 89.5% during the first round of optimizations. Note that the exergy calculation in eq 20 is based on the charge gas stream. What we are concerned with is actually the exergy gained or lost from an exchanger, i.e., the exergy change of the refrigerant stream on the other side of the exchanger, because this will be related to the refrigerant compressor work consumptions. Since a charge gas stream passing through an exchanger may release exergy or obtain exergy, eq 20 generalizes two scenarios with the help of the binary variable z. If the charge gas stream releases exergy, z will be 1, which means the actual released exergy will be the theoretical generated exergy from the cracked gas stream multiplied by the efficiency eu; otherwise, if the charge gas stream absorbs exergy, z will be 0, which means the actual consumed exergy will be the theoretical gained exergy by the charge gas stream divided by

Ind. Eng. Chem. Res., Vol. 49, No. 12, 2010 Table 2. Exergy Efficiencies of Heat Exchangers in the Current Chilling Train System heat exchanger ID CB1 CB2 CB3 CB4 CB5 CB6 CB7 CB8 CB9 CD1 CD2 CD3 CD4 CD5 average

heat duty (GJ/h)

exergy generated (GJ/h)

exergy consumed (GJ/h)

exergy loss (GJ/h)

4.207 3.395 5.565 16.806 10.451 12.649 4.587 2.285 2.325 5.716 29.616 15.995 3.192 0.627

0.486 0.963 0.938 6.958 5.163 8.410 4.287 2.919 3.926 1.029 9.235 3.803 1.431 0.524

0.380 0.834 0.786 6.442 4.924 8.049 4.145 2.625 3.389 0.699 8.765 3.526 1.376 0.507

0.106 0.129 0.152 0.516 0.239 0.361 0.142 0.294 0.537 0.330 0.470 0.277 0.054 0.018

exergy efficiency 78.1% 86.6% 83.8% 92.6% 95.4% 95.7% 96.7% 89.9% 86.3% 67.9% 94.9% 92.7% 96.2% 96.7% 89.5%

the efficiency eu. Therefore, eq 21 gives the logic constraint for z, and eq 22 gives the actual exergy consumption/generation from the exchanger. 3.3. Compressor Model. For a compressor u, suppose its inlet and outlet streams are j and k. The model equations are generalized as below Fk ) Fj

(23)

Mk ) M j

(24)

Tk ) τu(Tj, Pj, Pk, Mk)

(25)

Hk ) ηk(Tk, Pk, Mk)

(26)

Sk ) φk(Tk, Pk, Mk)

(27)

Wu ) Fj(Hk - Hj)

(28)

where (j,k) ∈ {(S4,S5),(S14,S15),(S29,S30),(S40,S41)} and u ∈ CP ) {CP1,CP2,CP3,CP4}. Equation 25 shows that the outlet temperature of the u-th compressor, Tk, is a function of Tj, Pj,Pk, and Mk (represented by τu (•)). Equations 26 and 27 are functions to calculate the specific enthalpy and entropy. Equation 28 calculates the work consumption of the uth compressor. Note that Pk will be optimally determined by the chilling train synthesis model. 3.4. Flash Drum Model. For a flash drum u, suppose its inlet stream is j; the vapor outlet stream is k; and the liquid outlet stream is l. The model equations are generalized as below Fk ) FjRu

(29)

5791

Fl ) Fj(1 - Ru)

(30)

Pk ) Pl ) Pj - ∆Pu

(31)

Tk ) Tl ) τu(Tj, Pj, ∆Pu, Mj)

(32)

Ru ) Fu(Tj, Pj, ∆Pu, Mj)

(33)

Mk ) νu(Tj, Pj, ∆Pu, Mj)

(34)

Ml ) (FjMj - FkMk)/Fl

(35)

Hk ) ηk(Tk, Pk, Mk)

(36)

Hl ) ηl(Tl, Pl, Ml)

(37)

Sk ) φk(Tk, Pk, Mk)

(38)

Sl ) φl(Tl, Pl, Ml)

(39)

where (j,k,l) ∈ {(S6,S8,S9),(S7,S10,S11), ..., (S48,S49,S50)} and u ∈ {FLS1A,FLS1B, ..., FLS5}. ∆Pu is the fixed pressure drop of the uth flash drum. Ru is the molar vapor fraction, which is a function of Tj,Pj,∆Pu and Mj (represented by Fu(•) in eq 33). The outlet temperatures (Tk and Tl) and the vapor stream composition (Mk) are functions of Tj,Pj,∆Pu and Mj (represented by τu(•) and νu(•) in eq 32 and eq 34, respectively). On the basis of the mole balance, the liquid stream’s compositions can be obtained by eq 35. Similarly as aforementioned, the specific enthalpy and entropy of the vapor and liquid streams can be calculated by eqs 36 through 39, respectively. 3.5. Valve Model. In the chilling train superstructure, valves are only located at the outlet liquid streams of FLS4 and FLS5. They are mainly used to reduce the stream’s temperature through throttling expansion. After passing through such a valve, the pressure of a stream decreases; a part of the liquid is vaporized; and the stream temperature decreases. For a valve u, suppose its inlet and outlet streams are j and k. The model equations are generalized below Fk ) Fj

(40)

Mk ) Mj

(41)

Tk ) τu(Tj, Pj, Pk, Mk)

(42)

Hk ) H j

(43)

Sk ) φk(Tj, Pj, Pk, Mk)

(44)

where (j,k) ∈ {(S45,S53),(S49,S54)} and u ∈ {VA1,VA2}. Pk is the specified valve outlet pressure. The outlet stream’s temperature is the function of the inlet temperature, inlet pressure, outlet pressure, and outlet compositions. Because of the throttle valve expansion, Hk should be equal to Hj. Sk is the entropy of the outlet flow, which is the function of the inlet temperature, inlet pressure, outlet pressure, and outlet compositions (represented by φk(•)). 3.6. Mixer Model. For a mixer u, suppose its inlet streams are k and l, and the outlet stream is j. The mixer model equations are generalized below Fj ) Fk + Fl

(45)

5792

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

Tj ) Tkyi + Tl(1 - yi)

(46)

Pj ) Pkyi + Pl(1 - yi)

(47)

Mj ) yiMk + (1 - yi)Ml

(48)

Hj ) Hkyi + Hl(1 - yi)

(49)

Sj ) Skyi + Sl(1 - yi)

(50)

where (j,k,l) ∈ {(S21,S23,S26), (S37,S42,S55), (S35,S55,S43)} and u ∈ {MX1,MX2}. Note that the equations are associated with binary variables of yi. 3.7. Process Specification Constraint. The vapor stream of FLS5 contains mainly hydrogen, which is used as a hydrogen source for other facilities. To meet the hydrogen purity specification in downstream facilities, a mole concentration constraint is included. MH2,S52 g MHSpec 2

(51)

where MHSpec is the hydrogen purity specification. 2 3.8. Objective Functions. The chilling train synthesis model is a multiobjective model. It has three objectives: to minimize total energy consumption, to minimize ethylene loss, and to maximize hydrogen recovery amounts, which are shown below min J1 )

∑B

u

u∈HX

+

∑W

u

(52)

u∈CP

min J2 ) FS43MC2H4,S43

(53)

max J3 ) FS52MH2,S52

(54)

The first objective function is to minimize the total energy consumption. For the chilling train system, the total energy consumption is characterized by the summation of exergy consumption from all the heat exchangers and work consumption from all the compressors. The second objective function is to minimize ethylene loss. In the chilling train system, the ethylene in the stream S43 is accounted for as a loss because the ethylene in this stream can no longer be recovered. The third objective function is to maximize the hydrogen recovery amount. In the chilling train system, only the hydrogen in stream S52 is used as the hydrogen source. The hydrogen contained in the other outlet streams can only be burned as the fuel gas. 4. Thermodynamic Data Preparation and Formula Estimation As shown in the synthesis model equations, there are many thermodynamic property functions that will be calculated, such as the stream enthalpy, entropy, vapor fraction, and stream composition calculations. These property functions will be frequently used during the optimization. If the rigorous formulas are employed, the optimization complexity will be significant. On the other hand, commercial simulators, such as Aspen Plus, can easily and reliably calculate these thermodynamic property data. Thus, this paper employs Aspen Plus simulations to generate thermodynamic data and then regresses the data with simplified property formulas. After that, the obtained property functions are integrated into the synthesis model. Certainly, the data regression quality should be sufficient; meanwhile, the optimization results based on the regressed property functions need to be examined again with a rigorous process simulation.

In the synthesis model, the following property calculations need simplification: the specific enthalpy and entropy of process streams, vapor fraction and vapor compositions of flash drums, and the outlet temperatures of compressors and valves. Obtaining the reliable property estimation is not an easy task. The key input variables of a property function and their ranges should first be determined. For instance, the input variables of the vapor fraction of the third-stage flash drum are the inlet temperature, pressure, and compositions of the six components. However, the molar fractions of propylene and propane can be omitted because they are less than 100 ppm in the inlet stream. The molar fraction of ethane can be obtained from the other three components’ compositions. Thus, the key input variables of vapor fraction are actually the inlet temperature, pressure, and compositions of three components. Second, rigorous Aspen Plus models are developed and executed to obtain sufficient data under different input scenarios. For the above vapor fraction example, each key input variable samples eight values. There are a total of 32 768 flash drum scenarios to be simulated. Such a large simulation task must be conducted in an automatic way. Thus, some automation scripts written by Java are embedded into the Aspen Plus simulator to conduct all the simulations automatically. To collect and process the simulation results effectively and efficiently, some Excel-based program scripts are also developed. When the needed thermodynamic data are ready, the next step is to regress the obtained data with appropriate functions. Three types of estimation functions are employed in this paper: linear functions, piecewise linear functions, and quadratic functions. 4.1. Linear Estimation Function. Linearization is the simplest and most widely used method for function approximation. In this paper, the function to calculate outlet temperature of compressor, Tk ) τu(Tj,Pj,Pk,Mk), is approximated by the flowing linear function Tk ) a1Tj + a2Pj + a3Pk + a4Mk,H2 + a5Mk,CH4 + a6Mk,C2H4 + a7Mk,C2H6 + a8Mk,C3H6 + a9Mk,C3H8 + b (55) where a1 through a9 and b are undetermined coefficients. Similarly, most of the streams’ specific enthalpy Hk ) ηk(Tk,Pk,Mk), specific entropy Sk ) φk(Tk,Pk,Mk), all flash drums’ vapor fraction Ru ) Fu(Tj,Pj,∆Pu,Mj), and all flash drums’ vapor compositions Mk ) νu(Tj,Pj,∆Pu,Mj) except CH4 and C2H4 compositions can use the linear equation format as shown in eq 55. In many cases, linear functions will give very accurate estimations. For instance, Figure 5 gives the comparison between rigorous simulation and linear estimation results of the outlet temperature of compressor CP1. The estimation function matches the rigorous simulation data very well. The coefficient of determination is 0.999784, and the residual sum of squares is only 0.2851. 4.2. Piecewise Linear Function Estimation. In the chilling train model, two specific enthalpy functions need piecewise linear estimation: the system inlet stream S1 and HX3’s outlet stream S7. Under this condition, their Hk ) ηk(Tk,Pk,Mk) can be approximated as Hk )

∑λ H m

m

k,m

+ a1Tk + a2Mk,H2 + a3Mk,CH4 +

a4Mk,C2H4 + a5Mk,C2H6 + a6Mk,C3H6 + a7Mk,C3H8 + b (56)

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

∑λ P

m k,m

) Pk

5793

(57)

m

∑λ

m

)1

(58)

m

where a1 through a7 and b are undetermined coefficients, which are obtained based on the regression of rigorous simulation data. m is the number of linearization functions. Pk,m and Hk,m are specified vertices of piecewise lines selected based upon rigorous simulation results. λm are the special order sets (SOS2) defined in GAMS, which means at most two variables within a SOS2 set can have nonzero values, and these two nonzero values have to be adjacent. Figure 6 gives piecewise linear estimation of the inlet stream S1’s specific enthalpy. The maximum relative error is only 0.26%. Figure 7 gives the piecewise linear estimation of stream S7. It shows that the specific enthalpy of stream S7 is piecewise linear to pressure and linear to temperature. 4.3. Quadratic Function Estimation. In the chilling train model, the functions of Hk ) ηk(Tk,Pk,Mk), Sk ) φk(Tk,Pk,Mk), Mk,CH4 ) νu,CH4(Tj,Pj,∆Pu,Mj), and Mk,C2H4 ) νu,C2H4(Tj,Pj,∆Pu,Mj) for the vapor streams of a flash drum can be accurately estimated through quadratic functions. For instance, the quadratic approximation for MCH4,k ) νCH4,u(Tj,Pj,Mj) can be written as Mk,CH4 ) a1Tj + a2Pj + a3Mj,H2 + a4Mj,CH4 + ... + a12TjPj + a13TjMj,H2 + a14TjMj,CH4 + ... 2 2 2 2 + a11(Tj) + a22(Pj) + a33(Mj,H2) + a44(Mj,CH4) + ...

(59) where a1, a2, a3, ... are coefficients for the linear terms. a12, a13, a14, ... are coefficients for the interactive bilinear terms. a11, a22, a33, ... are coefficients for the quadratic terms. These coefficients will be regressed based on sufficient simulation data. Similarly, Hk ) ηk(Tk,Pk,Mk), Sk ) φk(Tk,Pk,Mk), and Mk,C2H4 ) νu,C2H4(Tj,Pj,∆Pu,Mj) can be approximated using the equation format as shown in eq 59. Figure 8 gives the quadratic estimation of the CH4 molar fraction in the vapor stream of FLS1B. The maximum relative error is only 0.02%. Figure 9 gives the quadratic estimation of

Figure 5. Linear estimation of the temperature of stream S5.

Figure 6. Piecewise linear estimation of the enthalpy of stream S1.

specific enthalpy in vapor stream of FLS1A. The maximum absolute error is 5.9 kJ/kmol, and the maximum relative error is 2.2 × 10-6. 5. Redundant Variable Elimination The developed chilling train model has some simple equations such as Mk ) Mj. These simple equations will generate many redundant variables. To reduce the model size and complexity so as to improve the solving efficiency, some methods are used to remove simple equations and redundant variables. 5.1. Combining Flow Rate and Composition Variables. In the chilling train process, the stream flow rate and compositions will not change when passing through a heat exchanger, a compressor, or a valve. Thus, the downstream flow rate and compositions can directly use those of the very upper stream. In the superstructure shown in Figure 4, for instance, stream S13 goes through HX5, CP2, and HX6. The flow rate and compositions of streams S14, S15, and S16 are equal to S13. Thus, only variables of stream S13’s flow rate and compositions are needed. The flow rate and composition variables of streams S14, S15, and S16 are replaced. The related mass composition equations of unit HX5, CP2, and HX6 are also removed. 5.2. Separating Output Variables from the Model. For the synthesis model, since the composition of the inlet stream S1 is fixed, the composition of vapor streams at the first stage flash

5794

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

Figure 7. Piecewise linear estimation of the enthalpy of stream S7.

Figure 8. Quadratic estimation of CH4 molar fraction of stream S11.

drum (S9 and S11) is only a function of the temperature and pressure of the inlet streams. Thus, when calculating the specific enthalpy and entropy of all the streams between the first-stage and second-stage flash drums, they are just functions of temperature and pressure of FLS1A or FLS1B plus the fixed composition of stream S1. The composition variables of S9 and S11 and associated equations can be removed during the optimization. Meanwhile, since the flow rate, temperature, pressure, enthalpy, entropy, and composition of all the flash drums’ liquid streams are directly sent out of the system and have nothing to do with the optimization results, these variables and their associated equations will not be included in the synthesis model. Through these simplifications, the updated synthesis model will be easier to converge. Note that although the unnecessary data and equations are removed during the optimization they will still be calculated when the optimization results are obtained. That means they

are not included in the synthesis model, but their values will be calculated after optimization for analysis. 6. Case Study and Analysis The chilling train synthesis model is formulated with GAMS. It is solved by DICOPT, where CPLEX and CONOPT are used for solving an MILP master problem and an NLP sub problem, respectively.13-15 The model has 292 constraints, 297 continuous variables, and 21 binary variables. The model has three objective functions: minimizing energy consumption, minimizing ethylene loss, or maximizing hydrogen recovery rate. At the first step, the model is solved three times by using three objective functions individually. Table 3 gives the results of three cases. For comparison, Table 3 also gives the results of the current chilling trains system, which are obtained through the rigorous process simulation with Aspen Plus.

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

5795

Figure 9. Quadratic estimation of the enthalpy of stream S9. Table 3. Comparison Results Based on a Single Objective Function

heat exchanger exergy (GJ/h) compressor work (kW) total energy (kW) ethylene loss rate (kmol/h) hydrogen recovery rate (kmol/h) hydrogen concentration (%)

current system

minimize energy consumption only

minimize ethylene loss only

maximize hydrogen recovery rate only

33.51 1614 10923 3.82 695.2 94.4

30.28 524 9212 1.145 540.9 95.0

31.58 519 9290 0.55 536.8 95.0

37.39 4016 14403 3.953 948.9 95.0

Table 3 shows that when minimizing the energy consumption the results of total energy consumption, ethylene loss, and hydrogen recovery rate will be improved compared with those of the current system. When minimizing the energy consumption only, the total energy consumption will be reduced to 9212 kW from the current value of 10 923 kW. The ethylene loss will be reduced to 1.145 kmol/h from 3.82 kmol/h. However, the hydrogen recovery rate will also be reduced to 540.9 kmol/h from 695.2 kmol/h. When minimizing the ethylene loss only, the ethylene loss will be reduced to 0.55 kmol/h from 3.82 kmol/ h. The hydrogen recovery rate will be reduced to 536.8 kmol/h from 695.2 kmol/h. The total energy consumption will be reduced to 9028 kW from 10 923 kW. When maximizing hydrogen recovery rate, the hydrogen recovery will increase to 948.9 kmol/h from current 695.2 kmol/h, but the ethylene loss rate will also increase to 3.953 kmol/h, and the total energy consumption will increase to 14 403 kW. Since the maximizing hydrogen case consumes much more energy than the other two cases, it is not an effective objective function. Table 4 gives the heat duty and exergy consumption of each heat exchanger when the objective function of minimizing energy consumption is adopted. This table shows the difference between heat duty and exergy consumptions. For instance, HX3 has the heat duty of -68.67 GJ/h, which means this unit releases the heat rate of 68.67 GJ/h. However, its inlet and outlet stream temperatures are both lower than the ambient temperature. Thus, the unit consumes exergy of 22.64 GJ/h for the processing. On the other hand, unit HX12 has the heat duty of 4.50 GJ/h, which means the exchanger absorbs the heat of 4.50 GJ/h from the refrigerant. However, the unit generates exergy of 3.31 GJ/h. The reason for this is that the inlet stream subcools the refrigerant to a lower temperature, and the inlet stream itself is

Table 4. Exergy Changes Based on Minimizing Energy Consumption name

heat duty (GJ/h)

exergy (GJ/h)

HX3 HX7 HX11 HX12 HX13 HX14 HX15 HX16 HX17 HX18 VA1 VA2

-68.67 -11.82 -2.56 4.50 -6.49 -0.52 2.38 -1.59 2.02 3.00 0 0

22.64 7.51 2.36 -3.31 3.80 2.06 -2.03 2.55 -2.34 -1.78 -0.13 -0.14

heated from a lower temperature to the ambient temperature; thus, it will generate exergy. As indicated in the methodology framework, the chilling train synthesis model has some simplifications regarding the real chilling train process. To validate the optimization results, rigorous simulations with Aspen Plus have been conducted. Figure 10 is the optimized diagram when only the objective of minimizing energy consumption is adopted. Figure 11 gives the flowsheet and rigorous simulation results based on optimized designs and operation settings. The flowsheet shows that charge gas is refrigerated and separated by the first three flash drums at lower pressure (about 1.75 MPa). After that, the vapor flow is compressed to a higher pressure (about 3.4 MPa), and then it is refrigerated and separated by FLS4 and FLS5. The optimized flowsheet in Figure 11 has lower exergy consumption in comparison with the original flowsheet (see Figure 1) because the operational temperatures and pressures in the process are optimized to reduce the exergy consumption

5796

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

Figure 10. Optimized flowsheet when minimizing the energy consumption.

Figure 11. Flowsheet and rigorous simulation results based on the optimized result in Figure 10.

of heat exchangers. Meanwhile, the compressor is optimally allocated from originally the behind of the second flash drum to currently the behind of the third drum, which causes 68% of the compressor load reduction. Thus, the compressor work and the overall energy consumption are much smaller than the original one. In the optimized flowsheet, the temperature of the third flash drum is lower than in the original flowsheet, which causes more ethylene to be liquefied and recovered in the downstream process. Thus, the ethylene loss rate is lower than the original one. However, the operation temperatures of the fourth and fifth flash drums are lower than the original ones, which causes more hydrogen to be liquefied and discharged as the methane fuel gas. Thus, the hydrogen recovery rate is lower than the original one.

On the basis of the comparison of the original and optimal designs of the chilling train system, it can be seen that the number of major pieces of equipment does not change. Thus, the major capital costs of two designs will not change too much. Meanwhile, it is also noticed that the compressor’s load in the optimized system is less than that of the original system. Thus, the compressor size in the optimized system can be smaller than that of the original system, but the cold box of CB7 in the optimized system should be larger than that of the original system because it involves more refrigerant streams. These are obviously two contradictory aspects in terms of capital costs. Therefore, actual capital cost calculations should be conducted at the detailed design stage in the future instead of the current conceptual design stage.

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

5797

Table 5. Composition Comparisons of Simulation and Optimization Results

Vapor H2 CH4 C2H4 C2H6 C3H6 C3H8 Liquid H2 CH4 C2H4 C2H6 C3H6 C3H8

first stage flash drum (%)

second stage flash drum (%)

third stage flash drum (%)

fourth stage flash drum (%)

simulation optimization simulation optimization simulation optimization simulation optimization simulation optimization simulation optimization

41 43.63 41.4 40.29 15.9 14.57 1.2 1.10 0.4 0.29 0.1 0.17

49.8 48.29 46.7 48.08 3.4 3.51 0.1 0.12 0 0 0 0

61.9 61.72 37.9 38.13 0.2 0.15 0 0.001 0 0 0 0

80.8 81.10 19.2 18.88 0 0.02 0 0 0 0 0 0

95.9 95 4.1 5 0 0 0 0 0 0 0 0

simulation optimization simulation optimization simulation optimization simulation optimization simulation optimization simulation optimization

0.7 1.33 15.2 16.82 56.9 57.21 8.3 6.91 15.6 15.18 3.3 2.55

0.9 0.97 39.7 42.72 53.2 51.22 4.6 3.93 1.5 1.01 0.2 0.16

3.3 2.90 92.1 93.92 3.8 3.42 0.9 0.21 0 0 0 0

4.1 3.97 93.2 92.91 2.7 3.06 0 0.06 0 0 0 0

3.3 3.38 96.6 96.35 0.1 0.27 0 0.002 0 0 0 0

Table 6. Comparison of Heat Duty and Exergy on Optimization and Simulation Results heat duty (GJ/h)

HX3 HX7 HX11 HX12 HX13 HX14 HX15 HX16 HX17 HX18

exergy (GJ/h)

optimization

simulation

optimization

simulation

-68.67 -11.82 -2.56 4.50 -6.49 -0.52 2.38 -1.59 2.02 3.00

-68.12 -12.53 -2.87 4.88 -6.98 -0.68 2.11 -1.93 2.45 2.67

22.64 7.51 2.36 -3.31 3.80 2.06 -2.03 2.55 -2.34 -1.78

22.43 7.82 2.55 -3.95 4.02 2.48 -1.99 2.87 -2.81 -1.85

Tables 5 and 6 give comparisons of the rigorous simulation and optimization results. Table 5 gives the composition of vapor and liquid streams from all flash drums. Table 5 shows that the optimization results and the rigorous simulation results generally match well. Table 6 gives heat duty and exergy comparisons

fifth stage flash drum (%)

of heat exchangers. Since most of the heat exchangers are combined in the optimization model, the heat exchanger data in the simulation results are also combined. From Table 6, most of the heat duties and exergy data match well between optimization and rigorous simulation results. Figure 12 shows the optimized relation between the energy consumption and hydrogen recovery rate. In this figure, the energy consumption rate is the sole objective function. Ethylene loss rate is specified at different values as shown in the horizontal axis. Hydrogen recovery rate is not constrained. Under this situation, the energy consumption rate changes in the opposite direction from that of the hydrogen recovery rate. As shown, when the ethylene loss rate increases from 1.1 kmol/h to 3.95 kmol/h, the total energy consumption will increase from 9200 kW to 14 400 kW; meanwhile, the hydrogen recovery rate will increase from 540 kmol/h to 948 kmol/h. It should be noted that the proposed chilling train synthesis model is a nonlinear MINLP model. The solver of DICOPT is

Figure 12. Optimized relation between the energy consumption and hydrogen recovery rate.

5798

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

Figure 13. Pareto frontier of the triple-objective synthesis model.

employed, where CPLEX and CONOPT are used for solving the MILP master problem and the NLP sub problem. Thus, the possibilities of finding local optima with the proposed MINLP model do exist. However, if a reliable global optimization solver or solving strategy could be employed, the optimization results would most likely be improved and have more significant merits. Figure 13 gives the 3D Pareto frontier of the triple-objective optimization problem. The algorithm for generating the 3D Pareto frontier is as follows: In the optimization model, the ethylene loss rate and hydrogen recovery rate are treated as equality constraints, and the energy consumption is selected as the objective function. The ethylene loss and hydrogen recovery rates are sampled 100 points for each, and this helps generate 10 000 combinational cases. After running all these cases, the optimal energy consumption rates are generated accordingly. Then, the 3D Pareto frontier is developed based on the obtained data. The developed optimization model is suitable and easy to be implemented for this algorithm. Figure 13 shows that when the hydrogen recovery rate increases the energy consumption and ethylene loss rate increase too. When the ethylene loss rate increases, the energy consumption rate first increases and then decreases. Generally, the hydrogen recovery rate affects energy consumption more than the ethylene loss rate. Through this figure, the optimal total energy consumption can be estimated when the hydrogen recovery rate and ethylene loss rate are determined. 7. Concluding Remarks The chilling train is a very important section in ethylene plants. Different design and operational conditions, such as compressor location and flash drum temperature and pressure settings, will affect energy consumption and product loss rates. In this paper, a systematic methodology for the optimal design and operation of a chilling train system has been presented. A multiobjective MINLP model based on rigorous process simulation has been developed. Detailed technologies on superstructure representation, model data preparation, thermodynamic formula estimation, and redundant variable removal are also introduced. On the basis of a real case study, it shows that the reduction of energy consumption rate and ethylene loss rate and the increment of hydrogen recovery rate can not be simultaneously achieved based on the current chilling train system. However,

an improved conceptual design is that when the minimization of the energy consumption rate is selected as the sole objective function the optimized chilling train system can significantly reduce the total exergy-accounted energy consumption rate and ethylene loss rate with a mild sacrifice of the hydrogen recovery rate. Acknowledgment This work was supported by Texas Air Research Center (TARC). Nomenclature Sets and Indices i ) Index of binary variables j,k,l ) Indices of streams m ) Index of piecewise linearization functions u ) Index of equipments S ) Stream HX ) Set of heat exchangers CP ) Set of compressors Parameters and Coefficients FMax ) Upper bound of streams’ mole flow rate (kmol/h) ∆Pu ) Equipment’s pressure drop (MPa) T0 ) Reference temperature, 303.15 K MHSpec ) Hydrogen purity specification 2 a, b ) Undetermined coefficient eu ) Exergy efficiency of a heat exchanger Variables F ) Stream’s mole flow rate (kmol/h) H ) Specific enthalpy (kJ/kmol) S ) Specific entropy (kJ/(kmol · K)) Q ) Heat duty of heat exchanger (GJ/h) Bu ) Exergy consumption/generation of equipment (GJ/h) Bu0 ) Theoretical exergy consumption/generation of a stream W ) Compressor’s work (kW) P ) Pressure (MPa) T ) Temperature (°C) R ) Vapor fraction of flash drum y ) 0-1 variable indicating whether or not a candidate process is chosen

Ind. Eng. Chem. Res., Vol. 49, No. 12, 2010 z ) 0-1 variable indicating a flow is releasing or consuming exergy in a heat exchanger J ) Objective variable λ ) Special order set (SOS2) variable Vectors M ) Vector of components’ mole fractions Functions η ) Specific enthalpy function φ ) Specific entropy function τ ) Temperature function F ) Vapor fraction function ν ) Components’ mole fraction functions

Literature Cited (1) Liu, C.; Zhang, J.; Xu, Q.; Li, K. Cyclic scheduling for best profitability of industrial cracking furnace system. Comput. Chem. Eng. 2010, 34, 544–554. (2) SRI Consulting. World petrochemical report on ethylene; SRI Consulting: Menlo Park, CA, 2008. Available online at http://www.sriconsulting.com/WP/Public/Reports/ ethylene/. (3) Chiu, C. H.; Newton, C. L. Second law analysis in cryogenic processes. Energy 1980, 5, 899–904. (4) Colemenares, T. R.; Seider, W. D. Synthesis of cascade refrigeration systems integrated with chemical processes. Comput. Chem. Eng. 1989, 247–258.

5799

(5) Hurme, M.; Dohnal, M.; Ja¨rvela¨inen, M. Qualitative decision support system for cold box operation. Comput. Chem. Eng. 1994, S541–S545. (6) Dhole, V. R.; Linnhoff, B. Overall design of low temperature processes. Comput. Chem. Eng. 1994, S105–S111. (7) Castillo, F. J. L.; Dhole, V. R. Pressure Analysis of the Ethylene Cold-end Process. Comput. Chem. Eng. 1995, 19 (Suppl. 1), 89–94. (8) Wu, G.; Zhu, X. X. Design of integrated refrigeration systems. Ind. Eng. Chem. Res. 2002, 41, 553–571. (9) Lee, G. C.; Smith, R.; Zhu, X. X. Optimal synthesis of mixedrefrigerant systems for low-temperature processes. Ind. Eng. Chem. Res. 2002, 41, 5016–5028. (10) AspenTech, Inc. Aspen Physical Property System, Aspen plus version 2006.5; 2007. (11) GAMS; GAMS Development Corporation: Washington, DC, 1992. (12) Haseli, Y.; Dincer, I.; Naterer, G. F. Optimum temperatures in a shell and tube condenser with respect to exergy. Int. J. Heat Mass Transfer 2008, 51, 2462–2470. (13) CPLEX. Using the CPLEX callable library; CPLEX Optimization, Inc.: Incline Village, NV, 1995. (14) Drud, A. S. CONOPTs A Large Scale GRG Code. ORSA J. Comput. 1994, 6, 207. (15) 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 March 1, 2010 ReVised manuscript receiVed April 28, 2010 Accepted May 5, 2010 IE100455G