Multi-objective Optimization of the Operation of an Industrial Low

3182

Ind. Eng. Chem. Res. 2006, 45, 3182-3199

Multi-objective Optimization of the Operation of an Industrial Low-Density Polyethylene Tubular Reactor Using Genetic Algorithm and Its Jumping Gene Adaptations Naveen Agrawal, G. P. Rangaiah,* Ajay K. Ray,† and Santosh K. Gupta‡ Department of Chemical & Biomolecular Engineering, National UniVersity of Singapore, Engineering DriVe 4, Singapore 117576, Republic of Singapore

In this study, a comprehensive model for an industrial low-density polyethylene (LDPE) tubular reactor is presented. The model parameters are tuned using industrial data on the temperature profile, the monomer conversion and the number-average molecular weight at the end of the reactor, and estimates of the several side products from the reactor. Complete details of the model are provided. Thereafter, a two-objective optimization of this LDPE reactor is performed; the monomer conversion is maximized while the sum of the normalized concentrations of the three important side products (methyl, vinyl, and vinylidene groups) is minimized. Three variants of the binary-coded non-dominated sorting genetic algorithmsnamely, NSGA-II, NSGA-II-JG, and NSGA-II-aJGsare used to solve the optimization problem. The decision variables used for optimization include the following: the feed flow rates of the three initiators and of the transfer agent, the inlet temperature, the inlet pressure, and the average temperatures of the fluids in the five jackets. Also, the temperature of the reaction mass is constrained to lie below a safe value. An equality constraint is used for the number-average molecular weight (Mn,f) of the product, to ensure product quality. Pareto-optimal solutions are obtained. It is observed that the algorithms converge to erroneous local optimal solutions when hard equality constraints such as Mn,f ) desired number-average molecular weight (Mn,d) are used. Correct global optimal Pareto sets are obtained by assembling appropriate solutions from several problems involving softer constraints of the type Mn,f ) Mn,d ( an arbitrary number. Furthermore, the binary-coded NSGA-II-aJG and NSGA-II-JG perform better than NSGA-II near the hard end-point constraints. The solution of a four-objective problem (with each of the three normalized side product concentrations taken individually as objective functions) is comparable to that of the two-objective problem, and the former (more) computationally intensive problem does not need to be solved. Introduction Ethylene, along with initiators and telogens, is used to produce low-density polyethylene (LDPE) via free radical polymerization in a tubular reactor at extreme conditions, namely, 150-250 MPa and 325-625 K. A typical commercial reactor has several reactions, and heating and cooling zones, with intermediate additions of initiators, monomer, and solvent, so that the conditions of polymerization differ significantly in each zone. The single-pass conversion of ethylene in this reactor is reported to be 20%-35%. The LDPE produced in these reactors contains several short-chain branches (mainly ethyl and butyl groups), which are responsible for its lower crystallinity, density, melting point, tensile strength, etc.1 The minimization of these groups would improve the quality and strength of the polymer. Some vinyl and vinylidene groups (unsaturated) are also present on the polymer chains. These are undesirable, because they make the final product susceptible to cracking due to oxide formation. Optimum operation of reactors should attempt to minimize these side products too. Yet another important requirement for optimal reactor operation is the maximization of the production (maximization of the monomer conversion for a given feed flow rate), while producing a product that has a desired value of the * To whom correspondence should be addressed. Tel.: 65-65162187. Fax: 65-6779-1936. E-mail: [email protected]. † Present address: Department of Chemical and Biochemical Engineering, University of Western Ontario, London, Ontario, N6A 5B9, Canada. ‡ Present address: Department of Chemical Engineering, Indian Institute of Technology, Kanpur, 208 016, India.

number-average molecular weight (Mn) (or the weight-average molecular weight, Mw), so that the product has the desired physical properties. Several detailed studies have been reported on the modeling of LDPE reactors in the literature. Agrawal and Han2 studied the effects of various operating parameters on the performance of the reactor. They incorporated axial mixing to simulate the effect of pressure pulsing. However, Chen et al.3 showed that axial mixing can be neglected because the Peclet number is quite high, because of turbulence. Donati et al.4 and Yoon and Rhee5 also observed that axial dispersion has a negligible effect on the reactor performance. Hence, an ideal plug flow model can be used to describe LDPE reactors. A pressure drop of 10%30% along the reactor length is observed in industrial practice.6 This also affects the molecular weight of the polymer, because the propagation rate constant is dependent on the pressure.7 In fact, the molecular weight distribution of the product is broader because of this condition. Some workers2,5,8-12 have used the quasi steady-state approximation to avoid integrating stiff ordinary differential equations (ODEs). This currently does not need to be assumed, because of the availability of powerful algorithms and computers. The most interesting observation that can be made from a survey of these modeling studies is the significant discrepancies in the values of the rate constants. These have been alluded to by Gupta et al.9 In more-recent studies,6,12 the rate constants are estimated using industrial data. For example, Brandolin et al.6 simulated an industrial LDPE reactor and calculated the kinetic parameters using the temperature profiles and the properties of the product. Because industrial data are involved, these workers do not provide

10.1021/ie050977i CCC: $33.50 © 2006 American Chemical Society Published on Web 03/29/2006

Ind. Eng. Chem. Res., Vol. 45, No. 9, 2006 3183

complete information, for proprietary reasons. The same is true for the study of Zabisky et al.,12 who also did not provide exact values of their tuned kinetic parameters. Asteasuain et al.13 provided the design features, steady-state operating conditions, measured temperature profiles, and the monomer conversions and number-average molecular weights at the reactor exit for yet another industrial reactor. Again, a few important details are not provided. The dynamic model of Asteasuain et al.13 has been modified to the steady-state model for use in this study. The model incorporates the axial variation of physical properties and pressure, in addition to temperature and concentration, as well as several main and side reactions, e.g., intramolecular chain transfer, chain transfer to polymer, β-scission of secondary and tertiary radicals, etc. (the latter give the extent of long- and shortchain branching and the amount of unsaturation). We assume (and provide) reasonable values for all the missing information. Therefore, our model description is quite complete and useful for researchers. Best-fit values of several model parameters are obtained using the reported industrial data.13 A similar approach has been used earlier for LDPE,6,12,14,15 poly(ethylene terephthalate) (PET),16 and nylon 6.17 This model is then used to optimize the LDPE reactor operation. Several studies on the optimization of LDPE reactors have been reported. A variety of objective functions, decision variables, and constraints have been used.18 Yoon and Rhee5 used the maximum principle to obtain the optimum temperature policy required to maximize the monomer conversion at the exit. Polymer quality (through Mn, Mw, or the side products) is not a concern in this study. Brandolin et al.19 obtained the optimal temperature and initiator concentration profiles that maximized the monomer conversion, while using several constraints on the polymer properties (through Mn, the polydispersity index (PDI), etc.). Mavridis and Kiparissides11 maximized the productivity of ethylene while controlling Mn and PDI in the product, using the reactor wall temperature and the concentrations of the initiator and the chain transfer agent as decision variables. Kiparissides et al.20 tuned an industrial reactor and used almost the same objective function for the online optimization of an LDPE tubular reactor. Yao et al.21 used a genetic algorithm (GA) to obtain the optimal jacket temperature profile required to maximize the polymer production. No requirements on the product properties are considered. A few optimization studies have also been reported under unsteady (dynamic) operation of LDPE reactors. Asteasuain et al.13 first presented a dynamic model of an LDPE reactor and then obtained the optimal start-up policies. They maximized the outlet conversion and minimized the time required for it to stabilize, while forcing the properties at some desired values during start-up. The feed flow rates of the two initiators and of telogen are used as decision variables. Cervantes et al.22 presented a dynamic model for an entire LDPE plant with a feed mixture of ethylene, methane, butane, and impurities. They minimized the switching time between two steady states, corresponding to two different grades of polymer. In all the studies involving more than one objective, a weighted sum of the multiple objectives is used as a single, scalar objective function. This allows the use of simpler optimization algorithms, but the solution is dependent on the values selected for the weighting factors; thus, some degree of arbitrariness is involved. A more important disadvantage of the combination of the many objectives into a scalar quantity (“scalarization”) is that the algorithm may miss some optimal solutions.23 Recently, several multi-objective adaptations of GA

Figure 1. Schematic diagram of an industrial LDPE reactor.

that can solve such problems have become available. These have been used in the present study to obtain solutions of a few meaningful multi-objective optimization problems for an industrial LDPE reactor. GA24-26 is an extremely robust technique that mimics natural genetics, using operators such as reproduction, crossover, and mutation to guide the search in the feasible domain. It requires only the values of objectives and does not require any initial guesses. One of its recent adaptations, the elitist non-dominated sorting genetic algorithm (NSGA-II),26 can be used to solve multi-objective optimization problems. This algorithm and its earlier versions (NSGA/NSGA-I) have been applied successfully to optimize several industrially important systems.27 The performance of GA including NSGA-II has been further enhanced by incorporating one of several jumping gene (JG) adaptations.28-31 Very recently, multi-objective differential evolution was applied to optimizing an adiabatic styrene reactor49 that was solved earlier50 by NSGA for multiple objectives. In this study, the operation of an industrial LDPE tubular reactor is optimized using two objectives: maximization of the monomer conversion and minimization of the (weighted average value of the) concentration of the undesirable side products (methyl, vinyl, and vinylidene groups). The binary-coded NSGA-II26 and its JG adaptations, NSGA-II-JG28 and NSGAII-aJG,31 are used. Pareto-optimal solutions (sets of non-dominated or equally good solutions, in which, on moving from any one point to any other, one objective function improves while the other worsens) are obtained. A decision maker can be provided with these and he/she can use his/her industrial intuition to select any one of these points as the “preferred” solution. To the best of our knowledge, this is the first study using multiple objectives (without scalarization) for an industrial LDPE reactor. Reactor Modeling and Simulation Formulation. Figure 1 shows the industrial LDPE tubular reactor13 that has been simulated and optimized in this study. It is typical of several industrial reactors. Equations T1-1-T1-10 in Table 1 give the fairly general kinetic scheme, which captures all the molecular developments in the LDPE product. All the reactions are considered to be elementary except reaction T1-1, which is of the order of 1.1, with respect to oxygen.15 The elementary reaction mechanism includes initiation, propagation, termination by combination, thermal degradation, transfer to polymer, transfer to solvent, β-scission of secondary and tertiary radicals, and intramolecular chain transfer. The reactor model includes mass, energy, and momentum balances for a tubular reactor. Ideal plug flow conditions are assumed in the reactor and jacket sides (i.e., there are no radial temperature or concentration gradients in the tubular reactor and jackets, and no axial mixing). These underlying assumptions are valid for high Reynolds (Re) numbers3 and a very high Lt/Di ratio.12 (Here, Lt is the total reactor length and Di is the inside diameter

3184

Ind. Eng. Chem. Res., Vol. 45, No. 9, 2006

Table 1. Kinetic Scheme and Model Equations for the LDPE Reactora activity oxygen initiation peroxide initiation propagation termination by combination

reaction/equation number

reaction/equation ko Kinetic Scheme (x, y ) 0, 1, 2,...; i, j ) 1, 2,...) O2 + M 982R1(0) fmkdm

Im982R1(0)

(T1-1)

(for m ) 1, 2)

(T1-2)

kp

Ri(x) + M98Ri(x + 1)

(T1-3)

ktc

Ri(x) + Rj(y)98Pi+j-1(x + y)

(T1-4)

ktdt

thermal degradation

Ri(x + 1)98Pi(x) + R1(0)

(T1-5)

chain transfer to telogen or solvent chain transfer to polymer

Ri(x) + S98Pi(x) + R1(0)

(T1-6)

intramolecular chain transfer (short-chain branching) β-scission of secondary radical (vinyl group formation) β-scission of tertiary radical (vinylidene group formation)

ktrs

ktrp

Ri(x) + Pj(y)98Pi(x) + Rj+1(y)

(T1-7)

kbb

Ri(x)98Ri(x)

(T1-8)

kb1

Ri(x + 1)98Pi(x) + R1(0)

(T1-9)

kb

Ri(x + 1)98Pi(x) + R1(0)

(T1-10)

Model Equations dV V dF )dz F dz dCIm dV V (for m ) 1, 2) ) -kdmCIm - CIm dz dz dCO2 dV V ) -koCMCO21.1 - CO2 dz dz dCS dV V ) -ktrsCSλ00 - CS dz dz dCM dV 1.1 V ) -koCMCO2 - kpCMλ00 - CM dz dz

( ) ( ( ( (

V

dCMe dz dCVi

(T1-11)

)

(T1-12)

)

(T1-13)

)

(T1-14)

)

(T1-15)

dV ) kbbλ00 - CMe dz

(

(T1-16)

)

dV ) kb1(λ00 - CR1(0)) - CVi dz dz dCVid dV V ) kb(λ00 - CR1(0)) - CVid dz dz 2 dCR1(0) V ) 2koCMCO21.1 + 2 fjkdjCIj - kpCMCR1(0) dz j)1 V

{

(T1-17) (T1-18)

∑

}

dV ktcλ00CR1(0) + ktdt(λ00 - CR1(0)) + ktrsCS(λ00 - CR1(0)) - CR1(0) dz 4U(T T ) J dT FCPV ) + kpCMλ00(∆H) (from Asteasuain et al.14) dz 1000Di

{

(

)

}

2frFV2 dP dV + FV ) - 10-6 dz Di dz ∞

λnp )

(T1-20) (T1-21)

Moment Equations

∞

∑∑ i x C

(T1-19)

n p

Ri(x)

(for n ) 0, 1, . . . ; p ) 0, 1, 2, ...; from Katz and Saidel32)

(T1-22)

(for n ) 0, 1, ... ; p ) 0, 1, 2, ...; from Katz and Saidel32)

(T1-23)

i)1 x)0 ∞ ∞

µnp )

∑∑ i x C n p

Pi(x)

i)1 x)0

{

2 dλnp V ) 2koCMCO21.1δp0 + 2 fjkdjCIjδp0 + pkpCMλn,p-1 - ktcλ00λnp + dz j)1

∑

n

ktdt(λ00δp0 - λnp) + ktrsCS(λ00δp0 - λnp) + ktrp(λ00

}

dµnp dz

)

2

j

j,p+1

- λnpµ01) -

(T1-24)

(for n ) 0, 1, ... ; p ) 0, 1, 2, ...)

{ ∑∑ ∑ ktc

n

j)0

dV λnp dz V

∑( C )µ

n

j

p

(nCj)(jCk)(pCi)(-1)n-jλkiλj-k,p-i + ktdtλnp + ktrsCSλnp

j)0 k)0 i)0

+ ktrp(λnpµ01 - λ00µn,p+1) - µnp

}

dV dz

(T1-25) (for n ) 0, 1, ... ; p ) 0, 1, 2, ...)

Ind. Eng. Chem. Res., Vol. 45, No. 9, 2006 3185 Table 1 (Continued) Moment Closures

( ) ( )

µ03 ) µ00

µ02 µ01

3

µ13 ) µ10

µ12 µ11

3

( (

Mw ) 28 PDI )

(T1-27) Definitions

( )

V CM XM ) 1 Vin CMin Mn ) 28

(T1-26)

(T1-28)

) )

µ01 + λ01 µ00 + λ00

(T1-29)

µ02 + λ02 µ01 + λ01

(T1-30)

Mw Mn

SCB per 1000 CH2 )

(T1-31) 500CMe

Vinyl per 1000 CH2 )

500CVi

(T1-33)

µ01 + λ01

Vinylidene per 1000 CH2 ) a

(T1-32)

µ01 + λ01

500CVid

(T1-34)

µ01 + λ01

Data taken from ref 6, unless noted otherwise.

of the reactor.) The reaction mixture is assumed to be homogeneous (single-phase), i.e., an ethylene-polyethylene mixture behaves as a supercritical fluid in the range of the given operating conditions. The model is written using the axial length (z) as the independent variable, and the differential equations are integrated along the reactor length. Equations T1-11-T1-19 constitute the mass balance on each species in the reactor, whereas eq T1-20 represents the heat transfer from the reaction mixture to the coolant through jacket walls, followed by a momentum balance (eq T1-21). The characteristic equations (eqs T1-24 and T1-25) describe the growing and dead polymer concentrations in terms of bivariate moments (see eqs T1-22 and T1-23) of orders m and n. The complete set of model equations (eqs T1-11-T1-25)6,32 for steady-state operation are similar to those presented by Brandolin et al.6 Moment closure equations relating the third-order moments to the lower-order moments, based on the log-normal distribution,12 are also included in this table (eqs T1-26 and T1-27). The physical properties of the reaction mixture such as the density (F), viscosity (µ), and the thermal conductivity (K), vary along the axial location (eqs T2-1, T2-5, and T2-7, respectively, in Table 2). Linear additivity is assumed for these relationships.9 Table 29,11-13,33-38 gives the correlations that have been used. An average constant temperature of jacket fluid in each zone is assumed in this study. The jacket fluid normally flows countercurrently in industrial LDPE reactors. The “correct” modeling of such systems requires the solution of the coupled set of several ODEs for the inner (reacting) fluid and the energy balance of the outer (counter-currently flowing) fluid. This needs iterative solutions14 and is prohibitively time-consuming, and unsuited for optimization studies (where the model must be solved several times over for each chromosome in each generation). The assumption of a constant (average) temperature of the coolant, although not exact, circumvents this problem. In addition, because data are being “tuned”, errors that are associated with this assumption will be addressed by tuning. The coupled nonlinear ODEs that describe the reactor are integrated using the D02EJF subroutine in the NAG library.

This subroutine uses Gear’s technique39 to integrate the stiff equations. A tolerance (TOL) of 10-5 is used. A decrease in the value of this parameter to 10-8 changes the results only in the fourth decimal place. Oxygen is used as the initiator, and n-butane15 is used as the inert solvent in the feed stream, while initiator I1 (tert-butyl peroxypivalate8), and initiator I2 (tertbutyl 3,5,5 trimethyl-peroxyhexaonate10) are used as intermediate feeds, as shown in Figure 1. Details of the industrial system13 are given in Table 3. This table also includes reasonable values (assumed for this study) of the missing details13 of the reactor. The rate constants6 are provided in Table 4. When subjected to integration, the model gives the profiles of several molecular properties of LDPE (Mn, Mw, PDI, short-chain branching (SCB), and the vinyl and vinylidene group concentrations, as defined in eqs T1-29-T1-34 in Table 1) and the temperature, pressure, and concentrations of the monomer, telogen, and initiators, each as a function of the axial position, z. The model parameters are tuned using three sets of industrial data:13 the temperature, Tind(zj), which is read from the plot,13 with an accuracy of (2 K at several discrete points, zj; j ) 1, 2,..., 33 (with zj ) 0, 23, 51, 56, 79, 107, 135, 166, 180, 205, 266, 308, 350, 387, 429, 467, 509, 546, 597, 635, 677, 714, 751, 803, 855, 898, 929, 971, 1022, 1101, 1162, 1241, 1321 m for the five zones (0 m e z e 60 m, 60 m e z e 160 m, 160 m e z e 340 m, 340 m e z e 850 m, and 850 m e z e 1390 m)); and the values of the monomer conversion, the numberaverage molecular weight Mn, and the side-product concentrations in the final product. The sum-of-squares of the normalized error, I, between the model-predicted values and the industrial values,

I(u) )

(

1∑ i,j

)

Sind i (zj) Smi (zj)

2

(1)

is minimized. In eq 1, Si is the value of the ith property, and the superscripts m and “ind” represent the values predicted by the model and the industrial values, respectively. The term u

3186


Table 2. Property Correlations equation number

expression ∆H (kJ/kmol) ) -21500 × 4.1868 F (kg/m3) )

1 + 28CM(Vp - VM) Vp

FM (kg/m3) ) 745.18 - 0.51T Vp (m /kg) ) -6.793 × 10 3

-3

(from Gupta et al.9)

(T2-2a)

(T given in Kelvin; from Micheles and Geldermans33) -5

+ 1.558 × 10 T - 4.828 × 10

Vp (m3/kg) ) 9.135 × 10-4 + 5.731 × 10-7T WM )

(T2-1)

(from Asteasuain et al.13)

-8

T + 5.118 × 10 2

-11

T

3

(T2-2b) (for T e 426 K)

(for T > 426 K; from Parks and Richards34)

28CM F (from Zabisky et al.12)

[(η - ηo)ξ + 1]0.25 ) 1.023 + 0.23364Fr + 0.58533Fr2 - 0.40758Fr3 + 0.09332Fr4

(

)

Tc

ηoξ ) 4.61Tr0.618 - 2.04 exp(-0.449Tr) + 1.94 exp(-4.058Tr) + 0.1

[

( )

ηs ) η exp 2.00 + 0.001av

µ01 µ00

0.556

µ01 +

)]

Ev 1 1 R T 423

(

(for 0.1 < Fr < 3; all viscosities given in Pa s; from Poling et al.35)

(all viscosities given in Pa s; from Kiparissides et al.36)

Ev (kJ/kmol) ) -4.1868(500 - 560µ01) -1

K )) 0.03 × 41868

hi (W m-2 K-1) )

K(Nu) Di

Nu ) 0.026Re0.8Pr0.33 ho (W m

-1

ho (W m-2 K-1) ) ReJm )

De (m) )

(T2-8)

(from Mavridis and Kiparissides )

[ ()] Di L

2/3

(T2-10a)

(for 1398 < Re < 10 000)

(T2-10b)

(for m ) 1; from Coulson et al. )

KJm(0.026ReJm0.8PrJm0.33) De

DeFJmVJm µJm

(T2-7a) (T2-7b)

11

(for Re > 10 000; from Zabisky et al.12)

K ) ) 10 000

(T2-6c)

(T2-9)

Nu ) 0.166(Re2/3 - 125)Pr0.33 1 +

-2

(T2-6a) (T2-6b)

4

M′ (9.869Pc)

hw (W m

(T2-5)

1/6

3

-2

(T2-3b) (T2-4)

K (W m-1 K-1) ) 418.68(5.0 × 10-4WM + 3.5 × 10-4Wp)

ξ ) 107

(T2-3a)

38

(T2-11)

(for m ) 2, ..., 5; from Lacunza et al.37)

(for m ) 2, ..., 5)

(DJi2 - Do2) Do

Di 1 1 1 ) + + U hi hw Doho

represents the vector of parameters that are tuned. Binary-coded NSGA-II26 is used to minimize I. Estimation of Model Parameters. The rate constants in Table 4 are taken mostly from Brandolin et al.,6 except for the parameters Atrs, Etrs, and ∆Vtrs, which characterize the chain transfer to the telogen. These are taken from Asteasuain et al.14 The kinetic parameters (Ad1, Ed1, ∆Vd1, Ad2, Ed2, and ∆Vd2) for the two initiators are given6 as ranges, for proprietary reasons. The activation energies of the two initiators (Ed1 and Ed2) must be tuned using the industrial data, whereas average values (of the ranges given by Brandolin et al.6) of Ad1, ∆Vd1, Ad2, and ∆Vd2 are used in the model. Similar tuning6,12,14,15 of the kinetic parameters has been used earlier, too, to obtain the rate constants for ethylene polymerization. The first four industrial values (j ) 1-4) of the temperatures, Tind(zj), are not used for tuning, because no reaction is happening in this zone, and the reactor is only acting as a heat exchanger. It is observed that the model that has been tuned with only these two parameters underestimates the values of the monomer conversion (XM,f) and the number-average molecular weight (Mn,f) of the product. Moreover, the temperature peak in the first zone is overestimated

(T2-12) (T2-13) (T2-14) (T2-15)

while the peak in the second zone is underestimated. This suggests that we use additional parameters for an accurate prediction of the industrial data. Several simulations were performed to study the effect of the individual parameters on the results. The set of tuning parameters used is expanded in stages until a satisfactory agreement is attained. The activation energy that characterizes initiation by oxygen (Eo), and the activation energy that describes the propagation reaction (Ep) are incorporated in the earlier set (to give a total of four parameters). The former should help in improving the agreement of the temperature peak in the first zone (because oxygen affects the exothermic polymerization in this zone), whereas the latter affects the polymerization in the entire reactor and should help to improve the agreement of the results elsewhere. Unfortunately, this also does not help, because the temperature peak in the second zone is still underestimated. The parameters that influence the rates of heat transfer are then incorporated. Because steam is used to heat the reaction mixture in the first zone, the use of the high value of ho (10 000 W m-2 K-1)38 for condensing steam would not help. Hence, the four volumetric flow rates, VJ2 - VJ5, of

Ind. Eng. Chem. Res., Vol. 45, No. 9, 2006 3187 Table 3. Details of the Industrial LDPE Tubular Reactor Studieda quantity

numerical value

total reactor length, Lt inside diameter of reactor, Di wall thickness of reactor, t number of zones, Nz inner diameter of outer (jacket) wall, DJi axial lengths of zones, Lzm, (m ) 1, ..., 5) flow rate of monomer, FM flow rate of oxygen, Fo flow rate of telogen, FS flow rate of inert, Finrt flow rate of initiator-1, FI,1 flow rate of initiator-2, FI,2 flow rates of jacket fluids, VJm, (m ) 2, ..., 5)

1390 m 0.05 m 0.0254 mb 5 0.2032 mb 60, 100, 180, 510, 540 m 11 kg/s 6.8 × 10-5 kg/s 7.4 × 10-2 kg/s 0.22 kg/s 1.0 × 10-3 kg/s 1.6 × 10-4 kg/s 4.03 × 10-3,b 3.94 × 10-3,b 3.32 × 10-3,b 0.26 × 10-3 m3/sb 349.15 K 227.98 MPa 441.15, 498.15, 498.15, 441.15, 441.15 K 2.42834, 2.42834, 3.1401, 3.1401, 4.01933 kJ kg-1 K-1 0.0 kmol/m3

inlet temperature, Tin inlet pressure, Pin mean jacket temperatures, TJ,m, (m ) 1, ..., 5) specific heat of reaction mixture, CPm (m ) 1, ..., 5) initial conditions for moments, λnp, µnp (n ) 0, 1; p ) 0, 1, 2) a

Data taken from Asteasuain et al.13

Table 4. Rate Constants

Table 5. Bounds, Final Tuned Values, and Reported Values of the Parameters

(

b

k ) A exp -

Values assumed in this study.

)

E + 103P∆V RT

rate constanta

A

E (kJ/kmol)

∆V (m3/kmol)

ko kd1 kd2 kp ktc ktdt ktrsd ktrp kbb kb1 kb

1.6 × 1011 m3.3 kmol-1.1 s-1 1.0 × 1014 m3 kmol-1 s-1 c 1.0 × 1012 m3 kmol-1 s-1 c 4.0 × 105 m3 kmol-1 s-1 8.7 × 108 m3 kmol-1 s-1 7.7 × 109 m3 kmol-1 s-1 7.0 × 104 m3 kmol-1 s-1 5.2 × 104 m3 kmol-1 s-1 1.2 × 1010 m3 kmol-1 s-1 1.4 × 109 m3 kmol-1 s-1 4.4 × 109 m3 kmol-1 s-1

132168b 119929b 123117b 17431b 15282 79968 18406b 36844 60537b 84747b 70205b

-12.1 × 10-3 14.0 × 10-3 c 11.6 × 10-3 c -16.8 × 10-3 9.2 × 10-3 -10.0 × 10-3 0.0 -19.0 × 10-3 0.0 -9.90 × 10-3 -9.90 × 10-3

Note: Efficiencies (f1 and f2) of the two initiators are 0.98 (fitted in the range of 0.75-1.00)6 and 1.00,6 respectively. The efficiency of oxygen is assumed to be 1.00. a

Data taken from refs 6 and 14. Pressure P given in units of MPa, and temperature T given in Kelvin; R ) 8.314 kJ kmol-1 K-1. b Values obtained in the present study. These values differ from those of Brandolin et al.6 and Asteasuain et al.14 c Average values of the parameters based on the ranges reported by Brandolin et al.6 d Data taken from Asteasuain et al.14

the jacket fluid (that influence ho in those zones) are added on to the set of parameters used for tuning. Kiparissides et al.10 cited the value of av (eq 53 in their work) as being 0.0225. However, the same group of workers36 used a different value of 0.017 (eq 68) in another paper. This suggests that this parameter, which influences the viscosity of the reaction mass and, hence, the value of hi, also must be included in the set of tuning parameters. The tuning with two parameters in the first stage also led to the underestimation of Mn,f. Because the solvent (chain transfer) controls the molecular weight of the product without affecting the temperature of the reaction mass much, the corresponding rate parameter, Etrs, is also included for tuning purposes. Asteasuain et al.13 did not provide the concentrations of the side products (methyl, vinyl, and vinylidene groups) in the final polymer. The tuned model predicts results for these products

c

bounds

final tuned values

reported values6

125604 < Eo < 138164 117230 < Ed1 < 136071 117230 < Ed2 < 133977 14653 < Ep < 18003 14653 < Etrs < 20934 56521 < Ebb < 66988 71175 < Eb1 < 87922 62802 < Eb < 87922 0.5 × 10-3 < VJ2 < 7.0 × 10-3 b 0.5 × 10-3 < VJ3 < 7.0 × 10-3 b 0.5 × 10-3 < VJ4 < 5.0 × 10-3 b 0.1 × 10-3 < VJ5 < 5.0 × 10-3 b 0.009 < av < 0.0185

132168 119929 123117 17431 18406 60537 84747 70205 4.03 × 10-3 3.94 × 10-3 3.32 × 10-3 0.22 × 10-3 0.018

135945 94621-133140 94621-132721 17626 17253a 61964 79967 79967 1.2 × 10-3 c 1.2 × 10-3 c 1.2 × 10-3 c 1.2 × 10-3 c 0.017d

a Data taken from Asteasuain et al.14 b Values given in units of m3/s. Data taken from Yao et al.21 d Data taken from Kiparissides et al.36

that agree qualitatiVely with those reported by Brandolin et al.6 However, quantitatiVe agreement is necessary for a good model. Most workers6-10,40-42 have reported methyl (Me), vinyl (Vi), and vinylidene (Vid) contents in the product to be 25-30, 0.08-0.13, and 0.4-0.8, respectively, per 1000 CH2. Values ([Me]f, [Vi]f, and [Vid]f) of SCB/1000 CH2 ) 30, vinyl/1000 CH2 ) 0.1, and vinylidene/1000 CH2 ) 0.7 (from refs 8 and 9) are used as “industrial” values for tuning the model. The parameters, Ebb, Eb1, and Eb, associated with the corresponding reactions, also need to be used for tuning. The complete set of 13 model parameters used finally for tuning all the available results are: u ≡ [Eo, Ed1, Ed2, Ep, Etrs, Ebb, Eb1, Eb, VJ2, VJ3, VJ4, VJ5, av]. The exact objective function to be minimized is given by 33

I(u) )

( ) ( ) ( ) ( ) ( ) ( )

WT,j 1 ∑ j)5

W XM 1 -

Tind(zj)

2

Tm(zj) 0.30

+ WMn 1 -

2

XM,fm

W Vi 1 -

+ WMe 1 0.10

[Vi]fm

21900

2

+

Mmn,f

30

2

+

[Me]fm

2

+ WVid 1 -

0.70

[Vid]fm

2

(2)

In eq 2, Wj is the weighting factor associated with the normalized square error of the jth quantity. The values of these weighting factors are set to unity (1), except for WMn (15) and WXM (5), to give more emphasis for better prediction of Mn,f and XM,f. The lower and upper bounds and the final tuned values of these parameters are shown in Table 5. The values used for the computational parameters (the best values) in the binary-coded NSGA-II are given in Table 6. Figure 2 shows that the agreement between the model predictions (for the temperature profile, the monomer conversion and the number-average molecular weight at the end, and estimates of the several side products8,9,13) and the industrial results is quite good. Table 7 shows a comparison of these quantities at the end of the reactor. The temperature profile in Figure 2a shows that the first two zones are acting as preheating zones to heat the reaction mixture. Reaction is not occurring in these two zones. The free radical population is generated in the third zone after feeding in the first initiator, and these radicals react with monomer molecules to form a polymer via a polymerization reaction. A sharp temperature peak is observed in this zone, because of the exothermic nature of the polymerization reaction. Monomer conversion (XM) also shows a sudden jump, corresponding to

3188


Table 6. Values of the (Best) Computational Parameters Used in Binary-Coded NSGA-II, NSGA-II-JG, and NSGA-II-aJG

parameter Ngen Npop lsubstr lchrom laJG pc pm pJG Sr

NSGA-II for parameter for estimation optimization

NSGA-II-JG

100 100 30 390

600a 200 30 330

900a 200 30 330

0.8 0.01

0.95 0.015

0.9

0.95

0.9 0.005 0.8 0.9

Table 7. Comparison of the Model-Predicted Values to the Industrial Data Quantity at the Reactor Exit

NSGA-II-aJG (for two- and four-objective optimization) 700a 200 30 330 70 0.8 0.01 0.8 0.6

a

Generations required for convergence of the Pareto-optimal set for Mn,f ) 21 900 ( 200 kg/kmol.

this temperature peak, and almost all of the initiator is exhausted at this point. After that, XM remains practically constant until another initiator is added in the fifth zone. The second initiator is added in the fifth zone to boost the monomer conversion, as shown in Figure 2b. This initiator is also depleted soon after its introduction into the reactor. Note that the reaction mixture is cooled in the fourth zone to the optimal level for the half-life

c

property

industrial data

model prediction

monomer conversiona number-average molecular weighta CH3 groups (SCB) per 1000 CH2b vinyl groups per 1000 CH2c vinylidene groups per 1000 CH2b

0.30 21900 kg/kmol 30 0.1 0.7

0.2971 21901 kg/kmol 30.13 0.1 0.7

a Data taken from Asteasuain et al.13 Data taken from Goto et al.8

b

Data taken from Gupta et al.9

of the second initiator, which governs the maximum efficiency of the initiator. The number-average molecular weight (Mn) profile shows that it suddenly decreases in the third zone and becomes steady as the reaction mixture enters the cooling zone (Figure 2c). This decrease in molecular weight corresponds to the increase in the reaction temperature. Initiator decomposition is more temperature-dependent than chain growth, whereas chain termination is hardly affected.51 Thus, the free radical population is increased and they react with monomer molecules to form polymer chains;

Figure 2. (s) Model predictions, (1) industrial data (from Asteasuain et al.13), ([) industrial estimate from Goto et al.,8 and (9) industrial estimates from Gupta et al.9 for the LDPE reactor described in Figure 1.


however, these chains have a tendency to be smaller. Hence, the average molecular weight is reduced. However, the change in Mn at the beginning of the fifth zone is less abrupt than that observed in the third zone, which might be due to the presence of an already-formed polymer in the reaction mixture.6 Figure 2d shows that SCB in LDPE is observed after polymer is formed in the third zone and it remains almost constant until the second peroxide is injected in the fifth zone. The less-abrupt change in methyl-group concentration again accounts for the already-formed polymer change. The trends for vinyl and vinylidene concentrations can similarly be explained (see Figure 2e and f). Note that the tuning of the model parameters is performed using a single set of operating conditions, even though several sets of operating conditions should be used for reliable modeling. Jacket fluid flow rates are included in the set of tuned model parameters, which should have been part of the operating data used in the model parameter estimation. Both of these could have been avoided if more data were available13 or we had industrial data of our own. We did not provide the statistical inferences (confidence intervals and correlation structure) on the parameter estimates, because of high complexity of the problem and limited industrial data. However, overparametrization of the problem is avoided by the addition of extra parameters only after careful examination. We started by tuning only 2 parameters and progressively added more parameters to improve the predictions. In fact, the model was tuned with 26 parameters (pre-exponential factors and activation volumes, along with activation energies for the respective rate parameters); however, later, we discovered that almost the same results could be obtained with only 13 model parameters. Multi-objective Optimization of the LDPE Tubular Reactor Formulation. Multi-objective optimization (MOO) of the industrial LDPE reactor previously described is now performed using the model that has been developed. The two objective functions used are maximization of the monomer conversion (XM,f) and minimization of the (weighted average value of the) undesirable side product contents ([Me]f, [Vi]f, and [Vid]f). In industry, the number-average molecular weight of the product (Mn) is normally fixed, depending on the required polymer grade, but it can vary within an acceptable range about this fixed value. On the other hand, undesired side products should be minimized to improve the polymer quality and strength, although some may place them appropriately in constraints. Although thermodynamics and safety considerations govern conversion, other factors such as economics can lead to an optimal conversion. Also, MOO provides the opportunity to consider more objectives. Hence, we have chosen conversion and side products as objectives. Solution of this problem provides a range of solutions to the decision maker, who can choose one of them, depending on other considerations (such as safety, ease of operability, market demand, economics, etc.). The binary-coded NSGA,26 as well as its two jumping gene adaptations, NSGAII-JG28 and NSGA-II-aJG31 (both binary-coded), are used. In this study, the optimization assumes that the reactor geometry (reactor length (Lt), inside diameter (Di) and jacket diameter (DJi), monomer feed rate (FM), and flow rates of the jacket fluid are fixed and only the “operating” variables are used as the decision variables for the optimization, which are the inlet temperature (Tin), the feed flow rates (Fo, FS, FI,1, and FI,2) of oxygen, solvent and the two additional initiators added between the five average temperatures, TJ,1-TJ,5, of the jacket fluids and

the inlet pressure (Pin). In total, there are 11 decision variables. The temperature of the reaction mass along the axial length at an interval of 1 m is stored, and then the maximum temperature in the reactor (Tmax(z)) is determined. A local constraint, Tmax(z) e Tmax,d () 610.15 K), is imposed on the temperature to ensure safety, whereas the number-average molecular weight of the product (Mn,f) is constrained to lie (exactly) at a desired value of Mn,d of 21 900 kg/kmol.13 These two constraints are handled by incorporating them as penalty functions26,43 with weighting factors of w1 ) 109 and w2 ) 1010, respectively, in both of the objective functions. Both of these constraints are used in the normalized forms. Note that the constraint on temperature is an inequality constraint, whereas that on Mn is an equality constraint. The bracketed inequality constraint26 is used in the penalty term for the former. Lower and upper bounds are placed on each of the decision variables. The above problem is written mathematically as

(

Max G1 t XM,f - w1 1 Max G2 t

) 〈

〉

Mn,f 2 Tmax(z) - w1 1 Mn,d Tmax,d

2

(3a)

1 1 + ([Me]f/30) + ([Vi]f/0.1) + ([Vid]f/0.7) Mn,f 2 Tmax(z) 2 - w2 1 (3b) w2 1 Mn,d Tmax,d

(

) 〈

〉

subject to the following bounds:

323.15 K e Tin e 403.15 K

(from Brandolin et al.15) (3c)

5 × 10-5 kg/s e Fo e 10 × 10-5 kg/s (from Brandolin et al.15) (3d) 5 × 10-5 kg/s e FS e 0.5 kg/s (from Asteasuain et al.13) (3e) 5 × 10-5 kg/s e FI,1 e 5 × 10-3 kg/s (from Asteasuain et al.13) (3f) 5 × 10-5 kg/s e F I,2 e 5 × 10-3 kg/s (from Asteasuain et al.13) (3g) 383.15 K e T J,m e 543.15 K

(for m ) 1, ..., 5) (3h)

192.52 MPa e Pin e 253.31 MPa (from Brandolin et al.15) (3i) model equations given in Table 1

(3j)

The bracket operator shown in eqs 3a and 3b, 〈R〉, denotes the absolute value of the operand R, if the operand is negative. Otherwise, it returns a value of zero if R is nonnegative. The bounds (eqs 3c-3g, and eq 3i) for most of the decision variables have been chosen based on information in the literature.13,15 The bounds for the five average jacket temperatures (eq 3h) have been selected to give a range around the values reported by Asteasuain et al.13 The available NSGA-II codes maximize all the objective functions. Hence, a problem involving the minimization of a function, J, is converted to a maximization problem, using the transformation, G ) 1/(1 + J) (see the first term on the right-hand side in eq 3b). While solving eq 3, it was observed that the simulation was taking an excessive amount of CPU time for some chromosomes. These chromosomes were then studied individually, in detail.

3190


It was determined that this occurred only when FS was selected below a certain value. Under these conditions, the balance equations became extremely stiff in certain ranges of z, as reflected by the very high number of function evaluations called by the NAG library subroutine, D02EJF. Hence, the lower bound of FS was increased. Because FS also affects the other decision variables, the bounds of these also were changed from those selected initially (eq 3). Rajesh et al.44 had encountered a similar problem in the MOO of steam reformers. One of the decision variables, (H/C)in, had to be constrained to lie within a certain range of values selected by GA for two other decision variables, viz., the inlet temperature and (S/C)in, to avoid getting negative values of the intrapellet mole fractions (chromosome-specific bounds). The modified bounds for the decision variables used in place of those in eq 3 are given below:

323.15 K e Tin e 403.15 K

(4a)

5 × 10-5 kg/s e Fo e 10 × 10-5 kg/s

(4b)

2 × 10-2 kg/s e FS e 0.5 kg/s

(4c)

5 × 10-5 kg/s e FI,1 e 5 × 10-3 kg/s

(4d)

5 × 10-5 kg/s e FI,2 e 5 × 10-3 kg/s

(4e)

413.15 K e TJ,1 e 543.15 K

(4f)

473.15 K e TJ,2 e 543.15 K

(4g)

473.15 K e TJ,3 e 543.15 K

(4h)

413.15 K e TJ,4 e 543.15 K

(4i)

413.15 K e TJ,5 e 543.15 K

(4j)

182.39 MPa e Pin e 248.25 MPa

(4k)

The solution of eq 3 with these bounds overcomes the problem of excessive CPU time. Results and Discussion The solution of the multi-objective optimization (MOO) problem is obtained using an empirically determined best set of values of the several computational parameters. These are given in Table 6. The CPU time for a typical (reference) run of 700 generations on a P4 computer (3.0 GHz, 1 GB RAM) is 7 h, 5 min. This computer system can do 220 MFlops (million floating point operations per second), according to the LINPACK benchmark program (available at http://www.netlib.org) for a matrix of the order of 500. Some solutions, perhaps local optimal (see Figure 3), are obtained with NSGA-II for Mn,f ) Mn,d ( 0 kg/kmol (a “hard” constraint) rather than a Pareto-optimal set of solutions. We then relax (“soften”) the end-point constraint and allow Mn,f to lie within a small range of Mn,d (well within the experimental error of (10% for molecular weights), in particular, Mn,f ) Mn,d ( 200 kg/kmol, Mn,f ) Mn,d ( 20 kg/kmol, and Mn,f ) Mn,d ( 2 kg/kmol. Interestingly, Pareto sets of optimal points are obtained with excellent spreads (see Figure 3a). The Pareto sets for the first two problems are superimposed, giving confidence on the solutions. Because it is difficult to distinguish the overlapping points, these results are re-plotted in Figure 3b, using vertical displacements of 0.2. It is observed from Figure 3a that the results for the Mn,f ) Mn,d ( 0 kg/kmol case are quite far from those for the

Figure 3. (a) Converged solutions for several end-point constraints on the number-average molecular weight (Mn,f) using NSGA-II. Numbers in parentheses refer to the number of generations. (b) The results of Figure 3a are replotted with vertical shifts of 0.2 (i.e., the values of the ordinate for Mn,f ) 21 900 ( 20 kg/kmol are displaced vertically upward by 0.2, etc.)

other two cases. This creates a suspicion that the algorithm may be converging to local optima. The results obtained with a smaller range of Mn,f () Mn,d ( 2 kg/kmol) do not converge to the common Pareto set, even for a very large number of generations (it seems to converge to the same Pareto set, however, in 10 000 generations, when different algorithms, such as NSGA-II-JG and NSGA-II-aJG, are used; this is discussed later in this work). These possibly reflect the failure of the binary-coded NSGA-II to converge to the global optimal solution when one attempts to satisfy the constraint on Mn,f exactly (hard constraint). We attempt to improve the solutions, using NSGA-II-JG and NSGA-II-aJG. The best values of the computational parameters are again obtained empirically for each of the JG adaptations. These values are also listed in Table 6. Figure 4 shows the converged Pareto sets generated using the three techniquess NSGA-II, NSGA-II-JG, and NSGA-II-aJGsfor the case when Mn,f ) Mn,d ( 200 kg/kmol. Hereafter, the problem with Mn,f ) Mn,d ( 200 kg/kmol using NSGA-II-aJG is referenced as the reference case. All three techniques clearly give almost


Figure 4. (a) Converged Pareto-optimal sets for Mn,f ) 21 900 ( 200 kg/kmol using NSGA-II and its JG adaptations. Numbers in parentheses indicate the number of generations. (b) Results of Figure 4a replotted with vertical shifts of 0.2, as in Figure 3b.

Figure 5. (a) Pareto-optimal sets for Mn,f ) 21 900 ( 200 kg/kmol (reference case) using NSGA-II-aJG for different number of generations (indicated in parentheses). (b) Results of Figure 5a replotted with vertical shifts of 0.2, as described in Figure 3b.

similar Pareto sets, with a larger range and a good distribution of points (good spread), although NSGA-II converges to the Pareto-optimal set in the lowest number of generations (600), as compared to NSGA-II-aJG (700 generations) or NSGA-II-JG (900 generations). Figure 5 shows how the solutions converge for the reference case with NSGA-II-aJG. The solutions at the 700th generation clearly can be considered to be acceptable (converged). In contrast, Figure 6 shows that, when the MOO problem is solved using Mn,f ) Mn,d ( 20 kg/kmol, the convergence is extremely slow, and the same Pareto set is obtained only after ∼6000 generations. The solutions for several cases ((1100, (20, (2 kg/kmol) by both NSGA-II-JG and NSGA-II-aJG converged to the reference Pareto (in Figure 6). (The results for NSGA-II-JG are not shown in Figure 6, for clarity, but can be provided on request.) Because NSGA-II did not converge for the Mn,f ) Mn,d ( 2 kg/kmol case (Figure 3) while NSGA-IIaJG did (Figure 6, as did NSGA-II-JG), this indicates better performance of the latter technique(s) than NSGA-II when (near) hard end-point constraints are used. NSGA-II uses the concept of elitism, borrowed from nature, in which better chromosomes are copied to the next generation. However, diversity decreases

because of elitism. To avoid this phenomenon, Kasat et al.28 introduced jumping genes (JG) into NSGA-II. It seems that the relatively poor performance of NSGA-II for problems with (near) hard end-point constraints is due to the loss in diversity of chromosomes while NSGA-II-JG and NSGA-II-aJG introduce higher exploratory capability into the algorithm. Thus, they perform better than NSGA-II to solve difficult problems similar to that studied herein. In fact, Kasat et al.28 also observed that NSGA-II could not converge to the global Pareto-optimal set for the ZDT4 problem52 but NSGA-II-JG did, indeed, converge. We could not converge to the reference Pareto set for the Mn,f ) Mn,d ( 0 kg/kmol case by any of the three techniques (NSGA-II, NSGA-II-JG, and NSGA-II-aJG). This is shown in Figure 7. Our attempts to use different computational parameters to improve these results failed. All these indicate that the solutions (in Figures 3, 6, and 7) of the Mn,f ) Mn,d ( 0 kg/kmol case are local optima or that the methods have failed. NSGA-II could have failed due to the loss in diversity of the chromosomes in this particular problem. However, jumping genes (JG) are supposed to counteract this problem, but they also failed to create a diversified pool for chromosomes. Optimization methods generally are guaranteed to converge only when the

3192


Figure 7. Solutions for Mn,f ) 21 900 ( 0 kg/kmol using NSGA-II and its JG adaptations. Numbers in parentheses indicate the generation number. Results for NSGA-II-aJG (1050) and NSGA-II (1600) are the same as those in Figures 6 and 3, respectively.

Figure 6. (a) Converged Pareto sets for problems having different endpoint constraints on Mn,f using NSGA-II-aJG. Numbers in parentheses indicate the generation numbers. (b) Vertically shifted converged Pareto sets of Figure 6a (as described in Figure 3b).

underlying assumptions (such as continuity and convexity) are satisfied. Furthermore, stochastic methods including NSGA-II have no guarantee that they will converge in a limited number of generations/iterations. Convergence to the global optimum is even more difficult for complex problems with equality constraints. An interesting observation was made by identifying solutions that had Mn,f ) Mn,d ( 0.1 kg/kmol (the range of values actually present in the solutions for the Mn,f ) Mn,d ( 0 kg/kmol case) in the several Pareto sets in Figure 6. These are shown in Figure 8. In the Pareto set corresponding to Mn,f ) 21 900 ( 20 kg/kmol, two chromosomes (XM,f ) 0.3003, with normalized side products ) 2.2597; and XM,f ) 0.3401, with normalized side products ) 2.5106) are observed to have values of Mn,f ) 21900.04 kg/kmol and 21899.98 kg/kmol, respectively. Six solutions (shown by open triangles in Figure 8; two are very close) that have a value of Mn,f ) 21 900 ( 0.1 kg/kmol are identified in the near-converged Pareto set for Mn,f ) 21 900 ( 2 kg/kmol. No such solutions are observed in the Pareto sets for Mn,f ) 21 900 ( 1100 kg/kmol and Mn,f ) 21 900 ( 200 kg/kmol. The existence of such solutions, lying on the converged Pareto set and satisfying the end-point constraint of Mn,f )

Figure 8. Points having Mn,f ) 21900 ( 0.1 kg/kmol from among the Pareto sets of Figure 6a. These points are compared to the reference case.

21 900 ( 0.1 kg/kmol, and the fact that they were not “caught” by the algorithms when used with Mn,f ) 21 900 ( 0 kg/kmol, confirms the failure of the binary-coded NSGA-II and its JG variants for problems in which the end-point constraint on Mn,f is forced exactly. We suggest that solutions of such problems should be assembled by screening the solutions of several MOO problems with softer constraints of the type Mn,f ) Mn,d ( µ kg/kmol, where µ is an arbitrary number. Our earlier study45 on the MOO of the third-stage wiped-film PET reactor also involved similar hard end-point constraints on the molecular weight, and unique solutions were obtained. However, the MOO code was operated using different values of the random seed; a computational parameter and Pareto sets were then assembled. One must be extremely careful while solving MOO problems that involve hard end-point constraints on the molecular weight (and, possibly, other properties) before inferring that the solution is unique rather than a Pareto set. Indeed, it may be worthwhile to revisit some of the earlier studies involving hard end-point constraints and explore if the correct solutions are Pareto sets.


Figure 9. Pareto-optimal points and the corresponding decision variables and constraints for the reference case (Mn,f ) 21 900 ( 200 kg/kmol; NSGAII-aJG). Industrial data (denoted by a solid inverted triangle, 1) also are shown.

Figure 9 shows the Pareto-optimal set for the reference case, as well as plots of the decision variables and constraints corresponding to the several points in the Pareto set. Plots of the methyl, vinyl, and vinylidene contents are also shown. It is

observed from the Pareto set that higher monomer conversions can be achieved only with higher side products. The actual operating point (shown by the filled inverted triangle in Figure 9) for the industrial reactor gives much higher concentrations

3194


Figure 10. Temperature (T), monomer conversion (XM), and number-average molecular weight (Mn) profiles for (- - -) chromosome A, (s) chromosome B, and (- - -) chromosome C, shown in Figure 9a.

of the side products (for the same conversion), and so this type of study offers a scope of considerable improvement of industrial LDPE reactors. The plots of the decision variables reveal that the optimal solution is dependent, to a large extent, on four decision variables: FS, FI,1, FI,2, and Pin. When the flow rates, FI,1 and FI,2, of the two initiators are increased, higher conversions (at the cost of higher side products) are obtained, as expected. The effect of increasing these two flow rates must be counteracted by a decrease in FS (to maintain the molecular weight). Other decision variables are almost constant, with some amount of scatter. Validity of the entire range of optimization variables with the industrial data is shown by generating temperature, monomer conversion, and number-average molecular weight profiles for chromosomes A, B, and C chosen from the Pareto-optimal set (see Figure 9a). These chromosomes cover the complete range of non-dominated solutions in the Pareto-optimal set. The profiles in Figure 10 for these three chromosomes are comparable to those in Figure 2 for the industrial operation. The temperature profile for chromosome A in Figure 10a shows that the initial rate of polymerization in the fifth reactor zone decreases, because of the low initiator flow rate (FI,2). It reflects lower monomer conversion as against the chromosome C, as shown in Figure 10b. This is obvious, because a higher initiator concentration increases the concentration of free radicals and, subsequently the conversion of the monomer molecules.21 At this point, the butane flow rate (FS) has reached the lower bound, thus making the chain-transfer reaction less significant;22 consequently, the product molecular weight (Mn) increases. These molecular weight profiles for the chromosomes are shown in Figure 10c. Ehrlich and Mortimer46 have mentioned that an increase in pressure helps to reduce the SCB, vinyl, and vinylidene contents significantly. Figure 9 shows that higher inlet pressures (and, therefore, higher pressures in the entire reactor) are indicated at higher monomer conversions to keep the side product concentrations in check. The decrease of the concentration of the methyl group as the pressure increases (until it attains its upper bound; see Figures 9l and o) is attributed to the fact that the pressure dependence of the propagation rate constant is more significant than that of the branching reactions, as shown by Machi and co-workers.47,48 In the Pareto-optimal set, there exist some points for which monomer conversion attains values of ∼38% (see Figure 9a). These values are higher than the usual reported values (20%-35%) in the industrial reactors, and there might be some problems of high viscosities, reactor fouling, and even clogging at that level of conversion. In fact, this is the advantage with multiple non-dominated solutions (equally good

Figure 11. Effect of FM on the Pareto-optimal set (reference case).

points) in the Pareto-optimal set; the decision maker can choose a point (based on his/her industrial experience and intuition) that has an acceptable/lower monomer conversion (∼35%) for operating the plant. Generally, “higher-level qualitative considerations” are required to decide on the preferred solution, as suggested by Deb.26 We also studied the effect of monomer flow rate (FM) on the Pareto-optimal set. As expected, Figure 11 shows that monomer conversion can be reduced with higher amounts of monomer being fed to the tubular reactor. These results provide more options from which the decision maker may choose. A four-objective optimization problem is now studied in which each of the side products is taken as an independent objective function, i.e., there are four objective functions: maximization of the monomer conversion (XM,f); and minimization of the methyl ([Me]f), vinyl ([Vi]f), and vinylidene ([Vid]f) contents in the product, respectively, per 1000 CH2. This problem is formulated to study the formation of each side product individually, instead of the weighted average of the side products, which was done in the two-objective problem. These results will identify the need and potential for selectively reducing a particular side product. Binary-coded NSGA-II-aJG is used for MOO of the reference case (i.e., Mn,f ) 21 900 ( 200 kg/kmol). The decision variable set and local constraint on reactor temperature are the same as those described in the formulation of the two-objective optimization problem. Again, the constraints are handled by the penalty function method with weighting factors, w1 ) 109 and w2 ) w3 ) w4 ) 1010,


Figure 12. Results for the four-objective optimization problem in eq 5 (NSGA-II-aJG).

respectively, in all the objective functions. The mathematical formulation of above problem is written as follows:

Max G1 t XM,f

(5a)

Max G2 t

1 1 + ([Me]f/30)

(5b)

Max G3 t

1 1 + ([Vi]f/0.1)

(5c)

Max G4 t

1 1 + ([Vid]f/0.7)

(5d)

subject to bounds that have been described by eq 4,

eq 4

(5e)

and the following constraints:

Mn,f ) 21 900 ( 200 kg/kmol

(reference)

(5f)

Tmax(z) e 610.15 K

(5g)

model equations

(5h)

Figure 12 shows a solution of the optimization problem. The best set of computational parameters is given in Table 6. Because it is not possible to plot these results as a fourdimensional Pareto set, the four objective functions are plotted as a function of the chromosome number (after rearranging the results so that the conversion increases continuously with the chromosome number). Much more scatter is observed for the methyl group concentrations in this four-objective problem than was observed for the earlier two-objective problem. This scatter could not be reduced by a change of the computational parameters. However, the vinyl and vinylidene group concentrations show the increasing trends with monomer conversion (see Figure 12c and d). Similar trends were also observed in the two-objective problem. The optimal values of the individual (normalized) side-product concentrations have been added and plotted as a function of the monomer concentration in Figure 12f. The two-objective reference Pareto set (in Figure 5a) is also plotted. The results for the two- and fourobjective optimization problems are comparable. This suggests that one can easily combine the three side products into a single objective for this problem, and it is not necessary to solve a four-objective optimization problem. Nevertheless, the Paretooptimal set for the two-objective problem is superior to the

3196


Figure 13. (a) Flowchart of NSGA-II and the JG adaptations.53 (b) Evaluation of objective functions for optimization in Figure 13a, using a reactor simulation module.

Pareto-optimal solutions obtained from the four-objective problems. The optimal solutions for the four-objective problem show some scatter, but with a somewhat increased range of solutions. The computational parameters were varied, one by one, in the range of (10% of the best values in Table 6, to determine their effect on the results. It was observed that the Pareto-optimal set is not too sensitive to these variations in the test range. Therefore, these results are not provided here (but can be supplied on request). Conclusions A comprehensive mathematical model for the production of low-density polyethylene (LDPE) in high-pressure tubular reactors is developed. Complete details are provided. Tuned values of the several model parameters are obtained to get good agreement of model predictions with industrial data on the temperature profile, the monomer conversion, and the numberaverage molecular weight of the product, as well as estimates of the concentrations of the several side products. Thereafter, a

two-objective optimization study of the operating reactor (with constraints on the molecular weight and the temperature of the reaction mass) is conducted. Considerable improvement in the reactor performance is indicated. Pareto-optimal solutions are obtained. The present study suggests that solutions of problems involving hard (equality) end-point constraints should be assembled by obtaining solutions of several MOO problems with softer constraints, rather than by solving the problem only once, lest erroneous results are obtained. Furthermore, the binarycoded NSGA-II-aJG and NSGA-II-JG perform better than NSGA-II near the hard end-point constraints. The results of a four-objective problem (with each of the three normalized sideproduct concentrations taken individually as objective functions) are found to be comparable to that of a two objective problem in which these three are added together and taken as a single objective function. Acknowledgment Authors would like to thank the National University of Singapore-Indian Institute of Technology (NUS-IIT), Kanpur Student Exchange Program under which N.A. spent about six


months at IIT Kanpur. N.A. expresses his deepest gratitude and appreciation to the NUS for providing the financial support through research scholarship. Financial support from the Department of Science and Technology, Government of India, New Delhi [through Grant No. III-5(13)/2001-ET] is gratefully acknowledged by S.K.G. Appendix. Fundamentals of Genetic Algorithm and Binary-Coded Elitist Non-dominated Sorting Genetic Algorithm with the Jumping Gene Operators Genetic algorithms (GAs)24,25 are biologically inspired computing techniques that have a tendency to mimic the basic Darwinian concepts of natural selection. They use only the Values of the objective functions (not the gradients), and do not require initial guesses in the close proximity of the optimal solution. The randomly generated initial population is genetically processed following some basic principles of natural selection: (a) reproduction operation for identifying the candidates for the next generation; (b) crossoVer operation for the probabilistic exchange of genetic information between two randomly picked parents, facilitating a global search for the optimum in the process; and (c) mutation operation inducing a small, probabilistic change in the genetic makeup, resulting in a local search. Many problems in the engineering domain have a tendency to involve a situation that is multi-objective in nature, where more than one objective function must be optimized simultaneously. The optimal solution to a multi-objective function results in the Pareto-optimal set, rather than a unique solution. Hence, the simple GA is modified to a new algorithm known as non-dominated sorting genetic algorithm (NSGA) for finding the Pareto-optimal solutions. Two versions of this technique are available; NSGA-I and NSGA-II.26 The latter uses the concept of elitism, borrowed from nature, which gives the better parents a chance to be part of the next generation. In contrast, the likelihood of this happening in NSGA-I, which did not incorporate this concept, is very small. The three variants of NSGA-II, without and with the jumping gene (JG) operator, vary in the total number of operators. The binary-coded NSGA-II-JG has a new operator, introduced by Kasat et al.,28 along with the usual operators of NSGA-II. It probabilistically selects two sites in the chromosome string and replaces the in-between portion with a new, same-sized randomly generated binary string. The motivation behind borrowing the jumping genes concept in nature is to bring genetic diversity in the population.28 While in binary-coded NSGA-II-aJG, the second site in the chromosome is selected by the specified string length of jumping genes, as described by Bhat and Gupta.31 The working procedure of binary-coded NSGA-II with the jumping gene operators (JG and aJG) is described with the help of the flowchart displayed in Figure 13a. Initially, a binarycoded parent population P of size Npop is created using a randomnumber code. Each chromosome of this population is mapped into a set of real values of the decision variables. The reactor simulation program (Figure 13b) is used to compute the values of all the objective functions (for each chromosome). These objective function values are called fitness values, which are the indices of the merit of an individual. The population is classified into various fronts, based on the non-domination;26 each solution is assigned a fitness (or rank) equal to its nondomination level (1 is the best level, 2 is the next-best level, and so on). The crowding distance, Idist, for the chromosomes in all fronts then is calculated and all chromosomes are rearranged in their respective fronts in ascending order of the values of any one of their fitness functions. The crowding

distance of the ith solution in its front is the average side length of the cuboid (rectangle for two fitness functions) enclosing i, which just touches its nearest neighbors in the F-space. The better of the Npop chromosomes of P′ are selected in a new box P′′. This selection is done by comparing chromosomes (chosen randomly, irrespective of fronts) based on either Irank or Idist (if ranks are same for two randomly chosen chromosomes). The crossover and mutation operators (using probabilities, pc and pm respectively) then are used on the better parents to create an offspring population O of size Npop. At this stage, JG/aJG operation is checked (whether needed or not) on each chromosome (for example, 1001 | 10011 | 0) sequentially, as per specified JG probability (pJG). If JG operation is needed, a random number is generated between 0 and 1 and this number is multiplied by lchrom, the total numbers of binaries in the chromosome. The resulting number is rounded off, to convert it to an integer, which defines the position of the beginning of a transposon (for instance, at the end of the fourth binary in the aforementioned chromosome). Similarly, the second location is identified by generating another random number (for example, after the ninth binary in the chromosome). Whereas, in aJG operation, this second location is determined using the specified string length, laJG (for example, laJG ) 5; therefore, a bar is placed after the 4 + 5 ) ninth binary), of the jumping gene. The set of binaries between these two locations then are replaced by a new set of randomly generated binaries of the same size. For instance, we may get 1001 | 00100 | 0. All Npop members of P′′ and all the Npop members of O are copied into box R. The population R is of size 2Npop. These 2Npop chromosomes are reclassified into fronts using only nondomination and the best Npop are taken from box R′. These are put into box P′′′. This completes one generation and the process is repeated until the stopping criterion (Ngen,max) is met. The population from box P′′′ is copied into box P, for the next generation. Nomenclature A ) frequency factor (1/s; m3 kmol-1 s-1; m3.3 kmol-1.1 s-1) Ci ) concentration of the ith component (kmol/m3) CPi ) specific heat of the reaction mixture in the ith zone (kJ kg-1 K-1) De ) equivalent diameter of the jacket (m) Di ) inside diameter of reactor (m) DJi ) inner diameter of jacket wall (m) Do ) outer diameter of the inner (reactor) pipe (m) E ) activation energy (kJ/kmol) Ev ) activation energy for viscous flow (kJ/kmol) Fi ) flow rate of the ith component (kg/s) fm ) initiator efficiency of the mth peroxide fr ) friction factor Gi ) ith objective function in multi-objective optimization problem ∆H ) heat of polymerization (kJ/kmol) hi ) inside (the reactor) film heat-transfer coefficient (W m-2 K-1) ho ) outside (jacket side of reactor) film heat-transfer coefficient (W m-2 K-1) hw ) wall (reactor) heat-transfer coefficient (W m-2 K-1) Ii ) ith initiator Ii,dist ) crowding distance of ith chromosome Ii,rank ) rank of the ith parent/offspring chromosome K ) thermal conductivity of the reaction mixture (W m-1 K-1) k ) kinetic rate constant (1/s; m3 kmol-1 s-1; m3.3 kmol-1.1 s-1)

3198


L ) reactor length (m) laJG ) length of the replacing jumping gene lchrom ) total length (number of binaries) of a chromosome lsubstr ) length (number of binaries) of a substring representing a decision variable Lt ) total reactor length (m) Lzi ) axial length of ith zone (m) M ) monomer M′ ) molecular weight of ethylene (kg/kmol) Me ) methyl end group (short-chain branches) Mn ) number-average molecular weight Ngen ) generation number Npop ) total number of chromosomes in the population Nu ) Nusselt number Nz ) number of zones P ) reactor pressure at any axial position (MPa) Pr ) Prandtl number Pc ) critical pressure (MPa) Pi(x) ) dead polymer molecule with x monomeric units and i long-chain branches pc ) crossover probability pJG ) jumping probability for the JG operator pm ) mutation probability R ) ideal gas constant (kJ kmol-1 K-1) Re ) Reynolds number Ri(x) ) growing macroradical with x monomeric units and i long-chain branches S ) solvent (telogen) Sr ) seed for random number generator t ) wall thickness of reactor (m) T ) temperature of the reaction mass (K) Tc ) critical temperature (K) TJ,i ) jacket fluid temperature in the ith jacket (K) Tr ) reduced temperature U ) overall heat-transfer coefficient (W m-2 K-1) ∆V ) activation volume (m3/kmol) Vi ) vinyl group Vid ) vinylidene group VJi ) flow rate of the jacket fluid in the ith jacket (m3/s) VM ) specific volume of the monomer (kg/m3) Vp ) specific volume of the polyethylene (kg/m3) V ) velocity of the reaction mixture (m/s) VJ ) velocity of coolant in the jacket (m/s) wi ) weighting factor in the ith objective function WM ) monomer weight fraction Wp ) polymer weight fraction XM ) monomer conversion at any axial position z ) axial distance (m) Greek Symbols δij ) Kronecker delta η ) dense gas viscosity of the monomer (Pa s) η° ) low-pressure monomer viscosity (Pa s) ηs ) viscosity of the ethylene-polyethylene solution (Pa s) λnp ) n, p order moments for the chain-length distribution of macroradicals (kmol/m3); n ) 0, 1; p ) 0, 1, 2 µJ ) viscosity of the jacket fluid (Pa s) µnp ) n, p order moments for the chain-length distribution of the dead polymer molecules (kmol/m3); n ) 0, 1; p ) 0, 1, 2 ξ ) defined in eq 6 b in Table 2 ((Pa s)-1) F ) density of the reaction mixture (kg/m3) FJ ) density of the jacket fluid (kg/m3)

FM ) monomer density (kg/m3) Fr ) reduced density Subscripts b ) β-scission of a tertiary radical bb ) backbiting (intramolecular chain transfer) b1 ) β-scission of a secondary radical d ) desired value dm ) decomposition of the mth peroxide (initiator); m ) 1, 2 f ) (final) reactor exit I,m ) mth initiator, for m ) 1, 2 in ) reactor inlet inrt ) inert J ) jacket side M ) monomer max ) maximum o ) oxygen (initiator) p ) propagation S ) solvent tc ) termination by combination tdt ) thermal degradation trp ) chain transfer to polymer trs ) chain transfer to telogen or solvent Literature Cited (1) Luft, G.; Kampf, R.; Seidl, H. Synthesis Conditions and Structure of Low-Density Polyethylene. I. Short and Long Chain Branching. Angew. Makromol. Chem. 1982, 108, 203. (2) Agrawal, S. C.; Han, C. D. Analysis of the High-Pressure Polyethylene Tubular Reactor with Axial Mixing. AIChE J. 1975, 21, 449. (3) Chen, C. H.; Vermeychuk, J. G.; Howell, J. A.; Ehrlich, P. Computer Model for Tubular High-Pressure Polyethylene Reactors. AIChE J. 1976, 22, 463. (4) Donati, G.; Marini, L.; Marziano, G.; Mazzaferri, C.; Spampinato, M.; Langianni, E. Mathematical Model of Low-Density Polyethylene Tubular Reactor. In Chemical Reaction Engineering; Wei, J., Georgakis, C., Eds.; American Chemical Society Symposium Series 196; American Chemical Society: Washington, DC, 1981; p 579. (5) Yoon, B. J.; Rhee, H. K. A Study of the High-Pressure Polyethylene Tubular Reactor. Chem. Eng. Commun. 1985, 24, 253. (6) Brandolin, A.; Lacunza, M. H.; Ugrin, P. E.; Capiati, N. J. HighPressure Polymerization of Ethylene. An Improved Mathematical Model for Industrial Tubular Reactors. Polym. React. Eng. 1996, 4, 193. (7) Anspon, H. D. In Manufacture of Plastics, Vol. I; Smith, W. M., Ed.; Reinhold: New York, 1964. (8) Goto, S.; Yamamoto, K.; Furui, S.; Sugimoto, M. Computer Model for Commercial High-Pressure Polyethylene Reactor Based on Elementary Reaction Rates Obtained Experimentally. J. Appl. Polym. Sci. 1981, 36, 21. (9) Gupta, S. K.; Kumar, A.; Krishnamurthy, M. V. G. Simulation of Tubular Low-Density Polyethylene. Polym. Eng. Sci. 1985, 25, 37. (10) Kiparissides, C.; Verros, G.; Kalfas, G.; Koutoudi, M.; Kantzia, C. A Comprehensive Mathematical Model for a Multi-zone Tubular HighPressure LDPE Reactor. Chem. Eng. Commun. 1993, 121, 193. (11) Mavridis, H.; Kiparissides, C. Optimization of a High-Pressure Polyethylene Tubular Reactor. Polym. Process Eng. 1985, 3, 263. (12) Zabisky, R. C. M.; Chan, W. M.; Gloor, P. E.; Hamielec, A. E. A Kinetic Model for Olefin Polymerization in High-Pressure Tubular Reactors: A Review and Update. Polymer 1992, 33, 2243. (13) Asteasuain, M.; Tonelli, S. M.; Brandolin, A.; Bandoni, J. A. Dynamic Simulation and Optimization of Tubular Polymerization Reactors in gPROMS. Comput. Chem. Eng. 2001, 25, 509. (14) Asteasuain, M.; Pereda, S.; Lacunza, M. H.; Ugrin, P. E.; Brandolin, A. Industrial High-Pressure Ethylene Polymerization Initiated by Peroxide Mixtures: A Reduced Mathematical Model for Parameter Adjustment. Polym. Eng. Sci. 2001, 41, 711. (15) Brandolin, A.; Capiati, N. J.; Farber, J. N.; Valles, E. M. Mathematical Model for High-Pressure Tubular Reactor for Ethylene Polymerization. Ind. Eng. Chem. Res. 1988, 27, 784. (16) Bhaskar, V.; Gupta, S. K.; Ray, A. K. Modeling of an Industrial Wiped Film Poly(ethylene terephthalate) Reactor. Polym. React. Eng. 2001, 9, 71.

Ind. Eng. Chem. Res., Vol. 45, No. 9, 2006 3199 (17) Wajge, R. M.; Rao, S. S.; Gupta, S. K. Simulation of an Industrial Semibatch Nylon 6 Reactor: Optimal Parameter Estimation. Polym. Eng. Sci. 1994, 34, 1161. (18) Lee, K. H.; Marano, J. P. Free Radical Polymerization: Sensitivity of Conversion and Molecular Weights to Reactor Conditions. In Polymerization Reactors and Processes; Henderson, J. N., Bouton, T. C., Eds.; American Chemical Society Symposium Series 104; American Chemical Society: Washington, DC, 1979; p 221. (19) Brandolin, A.; Valles, E. M.; Farber, J. N. High-Pressure Tubular Reactors for Ethylene Polymerization Optimization Aspects. Polym. Eng. Sci. 1991, 31, 381. (20) Kiparissides, C.; Verros, G.; Pertsinidis, A. On-line Optimization of a High-Pressure Low-Density Polyethylene Tubular Reactor. Chem. Eng. Sci. 1994, 49, 5011. (21) Yao, F. Z.; Lohi A.; Upreti, S. R.; Dhib, R. Modeling, Simulation and Optimal Control of Ethylene Polymerization in Non-Isothermal, HighPressure Tubular Reactors. Int. J. Chem. React. Eng. 2004, 2, 1. (22) Cervantes, A.; Tonelli, S.; Brandolin, A.; Bandoni, A.; Beigler, L. Large-scale Dynamic Optimization of a Low-Density Polyethylene Plant. Comput. Chem. Eng. 2000, 24, 983. (23) Haimes, Y. Y. Hierarchical Analysis of Water Resources Systems: Modeling and Optimization; McGraw-Hill: New York, 1977. (24) Holland, J. H. Adaptation in Natural and Artificial Systems; University of Michigan Press: Ann Arbor, MI, 1975. (25) Goldberg, D. E. Genetic Algorithms in Search, Optimization and Machine Learning; Addison-Wesley: Reading, MA, 1989. (26) Deb, K. MultiobjectiVe Optimization Using EVolutionary Algorithms; Wiley: Chichester, U.K., 2001. (27) Bhaskar, V.; Gupta, S. K.; Ray, A. K. Applications of Multiobjective Optimization in Chemical Engineering. ReV. Chem. Eng. 2000, 16, 1. (28) Kasat, R. B.; Gupta, S. K. Multi-objective Optimization of an Industrial Fluidized Bed Catalytic Cracking Unit (FCCU) Using Genetic Algorithm with the Jumping Genes Operator. Comput. Chem. Eng. 2003, 27, 1785. (29) Simoes, A.; Costas, E. Transposition: A Biologically Inspired Mechanism to Use with Genetic Algorithm. In Proceedings of the Fourth International Conference on Neural Networks and Genetic Algorithms (ICANNGA99); Dobnikar, A., Steele, N., Pearson, D., Eds.; SpringerVerlag: Portoroz, Slovenia, 1999. (30) Man, K. F.; Chan, T. M.; Tang, K. S.; Kwong, S. Jumping Genes in Evolutionary Computing. In Proceedings of the Thirtieth Annual Conference of the IEEE Industrial Electronics Society (IECON) [at Busan, Korea, Nov. 2-6, 2004]; IEEE: Piscataway, NJ, 2004; Vol. 2, p 1268. (31) Bhat, S. A.; Gupta, S. K. Multi-objective Optimization of an Industrial Semi-batch Nylon-6 Reactor Using a New Jumping Gene Adaptation of Genetic Algorithm. Manuscript in preparation, 2005. (32) Katz S.; Saidel, G. M. Moments of the Size Distribution in Radical Polymerization. AIChE J. 1967, 13, 319. (33) Micheles, A.; Geldermans, M. Isotherms of Ethylene up to 3000 Atmospheres between 0° and 150°C. Physica 1942, 9, 967. (34) Parks, W.; Richards, R. B. The Effect of Pressure on the Volume, Thermodynamic Properties and Crystallinity of Polyethylene. Trans. Faraday Soc. 1948, 45, 203. (35) Poling, B. E.; Prausnitz, J. M.; O’Connel, J. P. The Properties of Gases and Liquids; 5th Edition; McGraw-Hill: New York, 2001. (36) Kiparissides, C.; Verros, G.; MacGregor, J. F. Mathematical Modeling, Optimization, and Quality Control of High-Pressure Ethylene Polymerization Reactors. J. Macromol. Sci.sReV. Macromol. Chem. Phys. 1993, C33, 437.

(37) Lacunza, M. H.; Ugrin, P. E.; Brandolin, A.; Capiati, N. J. Heat Transfer in a High-Pressure Tubular Reactor for Ethylene Polymerization. Polym. Eng. Sci. 1998, 36, 992. (38) Coulson, J. M.; Richardson, J. F.; Backhurst, J. R.; Harker, J. H. Coulson & Richardson’s Chemical Engineering. Volume 1, Fluid Flow, Heat Transfer and Mass Transfer, 5th Edition; Butterworth-Heinemann: Oxford, U.K., 1996. (39) Ray, A. K.; Gupta, S. K. Mathematical Methods in Chemical and EnVironmental Engineering; Thomson Learning: Singapore, 2001. (40) Gaylord, H. G.; Mark, H. F. Linear and Stereoregular Addition Polymers; Wiley: New York, 1959. (41) Kalyon, D. M.; Chiou, Y. N.; Kovenklioglu, S.; Bouaffar, A. HighPressure Polymerization of Ethylene and Rheological Behavior of Polyethylene Product. Polym. Eng. Sci. 1994, 33, 804. (42) Woodbrey, J. C.; Ehrlich, P. The Free Radical High-Pressure Polymerization of Ethylene. II. The Evidence for Side Reactions from Polymer Structure and Number Average Molecular Weights. J. Am. Chem. Soc. 1963, 85, 1580. (43) Edgar, T. F.; Himmelblau, D. M.; Lasdon, L. S. Optimization of Chemical Processes; 2nd Edition; McGraw-Hill: Boston, 2001. (44) Rajesh, J. K.; Gupta, S. K.; Rangaiah, G. P.; Ray, A. K. Multiobjective Optimization of Steam Reformer Performance using Genetic Algorithm. Ind. Eng. Chem. Res. 2000, 39, 706. (45) Bhaskar, V.; Gupta, S. K.; Ray, A. K. Multi-objective Optimization of an Industrial Wiped Film Poly(ethylene terephthalate) Reactor: Some Further Insights. Comput. Chem. Eng. 2001, 25, 391. (46) Ehrlich, P.; Mortimer, G. A. Fundamentals of the Free Radical Polymerization of Ethylene. AdV. Polym. Sci. 1970, 7, 386. (47) Machi, S.; Kawakami, W.; Yamaguchi, K.; Hosaki, Y.; Hagiwara, M.; Sugo, T. Structure and Properties of Polyethylene Produced by γ-Radiation Polymerization in Flow System. J. Appl. Polym. Sci. 1968, 12, 2639. (48) Machi, S.; Tamura, T.; Hagiwara, M.; Gotoda, M.; Kagiya, T. ShortChain Branching in γ-Radiation-Induced Polymerization of Ethylene. J. Polymer Sci.: Part A-1. 1966, 4, 283. (49) Babu, B. V.; Chakole, P. G.; Mubeen, J. H. S. Multiobjective Differential Evolution (MODE) for Optimization of Adiabatic Styrene Reactor. Chem. Eng. Sci. 2005, 60, 4822. (50) Yee, A. K. Y.; Ray, A. K.; Rangaiah, G. P. Multi-objective Optimization of an Industrial Styrene Reactor. Comput. Chem. Eng. 2003, 27, 111. (51) Luft, G.; Kampf, R.; Seidl, H. Synthesis Conditions and Structure of Low-Density Polyethylene II. Average Molar Mass and Molar Mass Distribution. Angew. Makromol. Chem. 1983, 111, 133. (52) Zitzler, E.; Deb, K.; Thiele, L. Comparison of Multiobjective Evolutionary Algorithms: Empirical Results. EVol. Comput. 2000, 8, 173. (53) Guria, C.; Bhattacharya, P. K.; Gupta, S. K. Multi-objective Optimization of Reverse Osmosis Desalination Units Using Different Adaptations of the Non-dominated Sorting Genetic Algorithm (NSGA). Comp. Chem. Eng. 2005, 29, 1977.

ReceiVed for reView August 29, 2005 ReVised manuscript receiVed December 9, 2005 Accepted February 24, 2006 IE050977I

Multi-objective Optimization of the Operation of an Industrial Low

Recommend Documents