A Hybrid Multiobjective Differential Evolution - ACS Publications

Ind. Eng. Chem. Res. 2009, 48, 11115–11132

11115

Optimization of Adiabatic Styrene Reactor: A Hybrid Multiobjective Differential Evolution (H-MODE) Approach Ashish M. Gujarathi* and B. V. Babu* Department of Chemical Engineering, Birla Institute of Technology and Science (BITS), Pilani - 333 031, (Rajasthan), India

In this study, a hybrid strategy of multiobjective differential evolution (hybrid-MODE) algorithm is proposed that consists of an evolutionary algorithm for global search and a deterministic algorithm for local search. To begin, the proposed algorithm is tested for its performance using a benchmark test function (KUR) as case study 1 with the nondominated sorting genetic algorithm-II (NSGA-II) (both binary- and real-coded versions), strength Pareto evolutionary algorithm (SPEA), Pareto archived evolutionary strategy (PAES), multiobjective differential evolution (MODE), and an improved strategy of MODE algorithms. Subsequently, the multiobjective optimization of an industrial adiabatic styrene reactor is carried out as case study 2, employing a prevalidated model using the hybrid-MODE algorithm and an improved strategy of MODE. Four cases (three sets of two-objective optimization, cases 1-3, and one set of three-objective optimization, case 4) are considered consisting of simultaneous maximization of styrene productivity, selectivity, and yield with four decision variables and two constraints. The proposed algorithm converges to a better set of nondominated solutions (possibly a Pareto front) as compared to the nondominated solutions obtained using NSGA and an improved strategy of MODE algorithms. The limitations of previous studies are reported in terms of key decision variables. The hybrid strategy of MODE is found to converge to the true Pareto front more rapidly (in fewer function evaluations), resulting in a well-diversified Pareto front as compared to the stand-alone evolutionary approach. 1. Introduction Styrene, C6H5CHdCH2 (a synthetic chemical), is used extensively in the manufacturing of plastics, rubber, and resins. Manufacturers around the globe use styrene-based resins to produce a wide variety of everyday goods ranging from cups and utensils to furniture; bathroom and kitchen appliances; hospital and school supplies; boats; sporting and recreational equipment; consumer electronics; automobile parts; and durable lightweight packaging.1 Industrially important and easily recognizable products obtained from styrene include polystyrene, acrylonitrile butadiene styrene (ABS), styrene-acrylonitrile (SAN) plastic, styrene-butadiene rubber (SBR), unsaturated polyester resins (commonly known as fiberglass), and two types of foams [extruded polystyrene (XPS) foam and expanded polystyrene (EPS) foam]. With such a large demand, even a slight improvement in the multiple process design decisions such as yield, selectivity, and product flow rate of the process would improve the profit of an organization significantly. Multiobjective optimization (MOO) refers to solving problems with more than one objective simultaneously. MOO using evolutionary algorithms (EAs) has gained popularity in the recent past because of the ability of EAs to yield a number of optimum solutions in a single run. Whether it is a single objective optimization or MOO, the decision maker is always interested in obtaining a solution that is suitable for the specific design requirements (i.e., a single solution). However, because of the multiobjective nature of the problem and the associated tradeoffs, the desired solution might vary according to the decision maker’s needs and choice. Thus, * To whom correspondence should be addressed. Tel.: +91-1596245073. Ext. 259/212. Fax: +91-1596-244183. E-mail: bvbabu@ bits-pilani.ac.in. Home page: http://discovery.bits-pilani.ac.in/discipline/ chemical/BVb/.

providing multiple solutions rather than a single optimum solution would be an advantage to the decision maker, so that one can have a choice of selecting one of several optimal solutions from the set of solutions. The specialty of such solutions is that, as one moves from one solution to another on the front (set of solutions), one objective is improved at the cost of a loss in another objective involved in the study. Such a set of solutions is referred as the Pareto optimal set, and the solutions in this set are nondominated with respect to each other.2 In the recent past, popular EAs such as the nondominated sorting genetic algorithm-II (NSGA-II) and its variants,2-5 strength Pareto evolutionary algorithm (SPEA),6 Pareto archived evolutionary strategy (PAES),7 and multiobjective differential evolution (MODE)8 and its different strategies have been proposed by various researchers. Several engineering problems have been successfully solved using these EAs. The complexity involved in engineering problems creates a demand for the development of new and efficient strategies for EAs. The multiobjective differential evolution algorithm is an extension of differential evolution (DE)9-13 for solving multiobjective optimization problems. Some successful applications of MODE are reported in the literature.8,14-17 MODE involves a simple nondominated sorting procedure coupled with the differential evolution algorithm. The nondominated sorting procedure is called after each generation to remove the dominated population points, thus refining the population. The key parameters of DE/MODE are NP, population size; CR, crossover constant; F, scaling factor; and MAXGEN, maximum number of generations.8,10 The crossover constant widens the points on the Pareto front, whereas the scaling factor maintains the diversity among the generated population points by local perturbation of the points. The trial vectors generated by application of these operators then compete with the target vector, and the better

10.1021/ie901074k CCC: $40.75  2009 American Chemical Society Published on Web 11/04/2009

11116

Ind. Eng. Chem. Res., Vol. 48, No. 24, 2009

one replaces the target vector in the next generation. Although MODE is a powerful algorithm that captures a global search space and generates the trial vector by its efficient recombination operation (crossover and mutation), it needs a sufficient minimum number of population points in order to converge to the global Pareto front. However, the nondominated sorting technique reduces the number of population points, thus limiting the number of nondominated solutions obtained in MODE algorithm. This problem of reducing the number of nondominated solutions was identified, and two additional variations of the MODE algorithm, namely, MODE II and MODE III, were proposed by Babu et al.16 Each method of optimization (deterministic and evolutionary), in isolation, has its own limitations and advantages over the other. It is the hybridization of the two in order to improve the efficiency that has been attempted in this work. For example, the evolutionary optimization algorithms start with multiple population points, and all of the population points usually converge to a single point (in the case of single-objective optimization) or a nondominated optimal set (in the case of multiobjective optimization) after the specified number of generations are met. However, the deterministic methods often start with a single initial guess (sequential simplex method starts with n + 1 points as initial simplex, where n is the dimension of the problem). The new point is created either by the method of gradients or by a certain perturbation law (in the case of direct search methods), and the new point is compared with the existing point. If the new point is found to be better than the current point, then it replaces the current point. The outcome of these methods is often dependent on the initial guess and the method of perturbation or the step size of the gradient. It might be possible for deterministic methods to converge to local minima, if the initial guess is selected wrongly, but at the same time, they have the advantage of faster convergence. Evolutionary-based optimization methods, because of multiple function evaluations in a single run, are more accurate at the cost of slower convergence. In this work, we propose a hybrid strategy of the multiobjective differential evolution algorithm taking advantage of both the deterministic local search method and the evolutionary-approach-based optimization method. The deterministic sequential simplex method is used for local search, whereas the evolutionary multiobjective differential evolution strategy is used for global search. The deterministic method is used as an accelerator that finds new superior points to converge to the Pareto front at a higher rate. Chiou and Wang18 applied a hybrid strategy of differential evolutionary algorithms to solve static and dynamic single-objective problems with application to the fed-batch fermentation process. The local search method was used to accelerate the algorithm, whereas the evolutionary algorithm was used for migration to a wider domain. The multiobjective genetic algorithm (MOGA) was hybridized with a neural network (NN) and a gradient-based optimizer to solve a complex design problem in fluid dynamics.19 A detailed working principle and application of the random-weighted genetic algorithm (RWGA) and multiobjective genetic local search (MOGLS) methods is available in the literature.1 Thus, although the general concept of hybrid algorithms is not new,20 this is the first attempt to hybridize the MODE algorithm with a local search method for solving MOO problems. The proposed algorithm is successfully tested on a benchmark test problem (KUR) as case study 1, and the performance measures such as convergence and divergence metrics are calculated. Subsequently, the multiobjective optimization of an industrial styrene reactor is carried out as case study 2 using the hybrid-MODE

(H-MODE) algorithm and an improved strategy of the MODE algorithm. Other algorithm-specific issues, such as the effects of (1) initial population size, (2) step size to generate a new point in the local search method, and (3) penalty parameter, on the Pareto front are also discussed. The following sections present a description of the working principle of H-MODE, followed by an analysis of the test function and the results for optimization of a styrene reactor. 1.1. Hybrid Multiobjective Differential Evolution. The sequential simplex method21-23 is used as a local search algorithm in the hybrid multiobjective differential evolution (H-MODE) strategy. Figure 1 shows the working principle of the H-MODE algorithm. The population set is initialized randomly within the specified bounds of the variables, and the corresponding costs (objective function value) are evaluated. Entire population points are preserved for the recombination operation. Three vectors from the initial population are selected at random to create a noisy random vector. Then, the trial (child) vector is generated by crossover between the target and noisy random vectors (parents). The cost of the trial vector is compared with that of the target vector for dominance. The winning vector is then used for the local search to obtain a better population point in the neighborhood using the sequential simplex algorithm. The local search deterministic methods cannot handle multiple objectives directly. They require a single objective function to evaluate. Therefore, to handle multiple objectives using a deterministic sequential simplex method, the overall objective (D) is calculated. The overall objective (D) of a population point for a min-min type of problem is defined as

( )

overall objective (D) ) abs

1

1/n

n

(1)

∑df

i i

I)1

n di ) 1 (n ) number of objective where 0 e di e 1 and ∑i)1 functions). For a two-objective min-min type of problem, the overall objective is given by

(

overall objective (D) ) abs

1 d1f1 + d2f2

)

1/2

(1a)

Here, to avoid priority-based weighting, the values of the two weights are assumed to be the same, that is, d1 ) d2 ) 0.5. This formulation is applicable provided that the values of both functions f1 and f2 are in the same range. If the range of function values is different with different orders of magnitude, say, f1 ∈ (0,1) and f2 ∈ (1 × 104,1 × 106), then the normalized function values are to be used, as given by the equation fi )

fi - fi,min fi,max - fi,min

(2)

where fi, fi,min, and fi,max are the actual, minimum, and maximum values, respectively, of the ith function. In the case of an industrial problem, where the minimum and maximum values of individual functions are not known, the current population is sorted with respect to individual functions to obtain minimum and maximum values of an individual function. In the present study, the actual function values are used, as the two function values are in the range of same order of magnitude. The values of all functions are evaluated corresponding to the best point obtained in the sequential


11117

Figure 1. Working principle of the hybrid-MODE algorithm.

simplex method. Even though a single-objective-based overall objective function is used in the sequential simplex method, the best point obtained in the sequential simplex method is checked for nondominance before it is accepted in the current population. If the point obtained using the sequential simplex method dominates any point in the current population, the dominated point gets replaced by the new point, and otherwise, a new selection procedure takes place. The local search is performed using the simplex method for all population points in a given generation, and hence, the same number of points is taken to the next generation. Because he local search method is applied on each point, not only does the convergence become faster, but the diversity of the solutions is also maintained. Therefore, it is not that all points are moving toward a single point but that a set of points (one from the vicinity of each previous point) are moving toward a global Pareto front. As the selection criterion is based on the concept of nondominance of the new point and not on the value of the overall objective function, the convergence of the algorithm is faster, thus ensuring that the solutions obtained are not trapped at a local Pareto front. The value of the overall objective function acts as a guideline

for determining new simplex value. n + 1 vertices (where n is the number of dimensions of the problem) are used in the local search sequential simplex algorithm. The overall objective function value (D) is evaluated at each of these vertices. In each iteration of the sequential simplex method, the simplex is changed to enhance the function values. Several possible transformations are considered, namely, reflection, contraction, and expansion. The reflection is the simplest one, where the vertex at which the function value is worst is reflected in the hyperplane formed by the other n vertices. The function value at the new vertex is then evaluated, and the process is repeated. However, it might be possible that the function value might be worst at the new vertex, in which case, if it is reflected, we will return to the old simplex and no progress can be made. It is also possible that the reflection process leads to a cyclic process. In this case, the simplex is extended or contracted in some direction. The iteration proceeds in such a manner that the volume of the simplex continues to decrease once appropriate conditions are satisfied, and ultimately, the volume might shrink to the required accuracy. The process is repeated until the termina-

11118


along the reactor length. According to the Arrhenius principle, the high-pressure steam acts as a source of heat to enhance the rate of the endothermic reaction. The steam is mixed with the reactants at the inlet of the reactor, so it also dilutes the mixture of reactant and products, shifting the position of chemical equilibrium toward the products. According to Le Chatelier’s principle, in the case of reversible endothermic reactions, a low pressure and a high temperature are preferred so as to obtain a high conversion and a high selectivity. The high-pressure steam also removes the coke that gets formed on the iron oxide catalyst through the water-gas shift reaction. Figure 2 shows a singlebed adiabatic reactor in which the entire amount of steam is mixed with ethyl benzene at the reactor inlet. The off-gases from ethyl benzene conversion reactions include light gases (such as hydrogen, carbon dioxide, carbon monoxide, methane, and ethylene) and two major aromatic byproducts (i.e., toluene and benzene). The six main reactions taking place in a catalytic tubular reactor are given by eqs 3-8 Figure 2. Schematic of single-bed adiabatic styrene reactor.

C6H5CH2CH3 T C6H5CHCH2 + H2

(3)

tion criterion is met. The following selection strategy is used in the hybrid-MODE algorithm:

C6H5CH2CH3 f C6H6 + C2H4

(4)

If the new point obtained using the hybrid-MODE algorithm dominates any point in the population { If the newly obtained point further improves the divergence metric or convergence metric { Then replace the dominated point and its cost with the new point and cost } } The newly obtained point replaces the dominated point in the population; otherwise, a new selection procedure takes place The stopping criteria can be of two kinds: (1) No new solution is added to the nondominated front for a specified number of generations, or (2) an upper bound is assigned to the number of generations. The stopping criteria can be a combination of the two as well. In this study, the second criterion is applied. The proposed hybrid-MODE algorithm is applied to two case studies (multiobjective optimization of the KUR test function and an adiabatic styrene reactor) in this study. The process description of styrene production using an adiabatic reactor is given in section 2, and subsequently, the problem formulations of both case studies chosen as applications are given in section 3. 2. Styrene Process Description and Reactor Kinetics Industrially, styrene monomer is produced in two ways, viz., through catalytic dehydrogenation of ethyl benzene and through ethyl benzene hydroperoxide. The former is the popular method of styrene production, where the reaction takes place in a catalytic tubular reactor. The reactants (ethyl benzene and a highpressure steam) are mixed at the inlet of the reactor. This configuration is known as adiabatic operation. Steam is also injected at an intermediate position, and the resulting configuration is called a steam-injected reactor or pseudoisothermal operation.24,25 The supply of steam plays an important role in controlling the kinetics of the reversible endothermic dehydrogenation reaction. Because the main reaction producing styrene is endothermic in nature, the temperature of the reactor decreases

C6H5CH2CH3 + H2 f C6H5CH3 + CH4

(5)

2H2O + C2H4 f 2CO + 4H2

(6)

H2O + CH4 f CO + 3H2

(7)

H2O + CO f CO2 + H2

(8)

At equilibrium, the reversible reaction results in 80-85% conversion of ethyl benzene in the given temperature range. Both kinetically and thermodynamically, high temperature favors the rate of dehydrogenation of ethyl benzene. At the same time, however, the byproducts (i.e., benzene and toluene; see eqs 4 and 5) are produced by thermal cracking at high temperatures, thus reducing the yield of the desired product. The high reaction temperature favors the productivity of styrene, but at the same time, unwanted byproducts are formed at high temperatures, thus decreasing the selectivity value. Steam over reactant (SOR) also affects the temperature, the pressure, and the initial concentration of reactants. It is necessary to feed an optimum value of SOR so as to obtain a compromise between the desired objectives. In accordance with Le Chatelier’s principle (as stated above), the pressure in the reactor plays an important role in the kinetics of the reactions involved. Thus, the temperature of ethyl benzene (TEB), the pressure (P), the SOR, and the initial 0 ) are considered as the decision ethyl benzene flow rate (FEB variables. Plug flow is considered for simulation of an adiabatic reactor. Axial dispersion is neglected in the pseudohomogeneous model, and any limitations of mass or heat transfer to the catalyst pellet or diffusion within the pellet are lumped into the rate constants. The differential equations employed in the current study are listed in Appendix A. A detailed description of the mathematical model, its validation, and kinetic data are reported in the literature.8,25 3. Problem Formulations for Case Studies on Multiobjective Optimization 3.1. Case Study 1: KUR Test Function. A difficult test problem was selected from the literature2 for evaluating the performance of the proposed H-MODE algorithm of this study. The KUR test problem (see Table 1) is an unconstrained, threevariable, two-objective test problem with a nonconvex and disconnected (discontinuous) Pareto front.


11119

Table 1. Multiobjective Optimization Problems, Bounds, and Constraints objective functions

constraints

decision variables

Case Study 1 n-1

f1(x) )

∑ [-10 exp(0.2√x

2

i

+ xi+12

-5 e x1, x2, x3 e 5

none

)]

i)1

n-1

f2(x) )

∑ (|x |

0.8

i

+ 5 sin xi3)

i)1

Case Study 2, Case 1 max J1 ≡ SST ) max J2 ≡ YST )

FST -

0 FST

FSTEAM < 454 kmol/h

0 FEB - FEB 0 FST - FST

850 K < Tmix1 < 925 K

0 FEB

550 K < TEB < 800 K 1 bar < P < 2.63 bar 7 < SOR < 20 0 27.56 kmol/h < FEB < 40.56 kmol/h

1000 K < TSTEAM < 1050 K Case Study 2, Case 2 max J1 ≡ SST )

FST -

0 FST

same as case 1

same as case 1

0 FEB - FEB

max J3 ≡ FST Case Study 2, Case 3 max J2 ≡ YST )

0 FST - FST

same as case 1

same as case 1

0 FEB

max J3 ≡ FST Case Study 2, Case 4 max J1 ≡ SST ) max J2 ≡ YST )

0 FST - FST

same as case 1

same as case 1

0 FEB - FEB 0 FST - FST 0 FEB

max J3 ≡ FST

3.1.1. Performance Measures. Two widely used metrics are used to compare the performance of the hybrid-MODE algorithm with those of other algorithms such as MODE III, NSGAII (both binary- and real-coded), SPEA, and PAES. The first metric, γ, measures the extent of convergence to a known Pareto set of solutions. The true Pareto front data obtained from Huang et al.26 is used for comparison. For each solution obtained with the proposed algorithm, the minimum Euclidean distance from it to the 500 uniformly distributed chosen solutions is computed. The average of these distances is used as the convergence metric, γ. The smaller the value of the convergence metric, the better the convergence to the true Pareto optimal front. The second metric, ∆, measures the extent of the spread achieved among the Pareto optimal solutions. The metric ∆ is calculated using the equation2 |N-1|

dfe + d1e + ∆)

∑ |d -d¯| i

i)1

dfe + d1e + |N - 1|d¯

(9)

where def + d1e represnts the sum of the Euclidean distances between the currently obtained extreme solutions and the

extreme solutions of the Pareto set and di is the Euclidean distance between the consecutive nondominated solutions. The parameter dj is the average of all Euclidean distances di, i ) 1, 2, ..., (N - 1), assuming that there are N solutions on the final nondominated front or the obtained Pareto optimal front. A higher value of the metric ∆ shows a worse distribution of solutions within the extreme Pareto optimal solutions. 3.2. Case Study 2: Multiobjective Optimization of Industrial Styrene Reactor. Simultaneous maximization of three objectives, namely, the productivity (FST), selectivity (SST), and yield (YST) of styrene are considered in this study. As higher productivity is a primary interest in every manufacturing industry, the maximization of productivity is considered as one of the objectives. Minimization of unwanted byproducts, such as toluene and benzene, is also important, as it results in higher costs for separating these impurities from the final product. To minimize the production of toluene and benzene (byproducts) in the reactor, the selectivity of styrene should be increased, and hence, the maximization of selectivity is considered as the second objective in this study. With an increase in the selectivity,

11120


the yield decreases; therefore, maximization of the yield of the reaction is considered as the third objective. Four decision variables, namely, feed temperature (TEB), pressure (P), steam-to-reactant ratio (also called steam over 0 ) are reactant, SOR), and initial ethyl benzene flow rate (FEB considered for optimization. The effect of this extra decision variable on the Pareto front is discussed in the Results and Discussion. These decision variables are among the operating variables of existing plants.27,28 The bounds for the decision 0 were chosen on the basis of variables TEB, P, SOR, and FEB information available in the literature.8,25 In this case study, simultaneous maximization of three objectives (namely, maximization of productivity, yield, and selectivity) is considered. Three cases of two-objective optimization (eqs 10-12) and one case of three-objective optimization (eq 13) are also considered in this study. Case 1: Simultaneous maximization of SST and YST

(10)

Case 2: Simultaneous maximization of FST and SST

(11)

Case 3: Simultaneous maximization of FST and YST

(12)

Case 4: Simultaneous maximization of FST, SST, and YST

(13)

The mathematical formulations of MOO problems considered in this study (KUR test function and adiabatic styrene reactor), along with the bounds on the decision variables and the associated constraints, are listed in Table 1. 4. Results and Discussion 4.1. Case Study 1: KUR Test Function. 4.1.1. Experimental Runs and Parameter Setting. For comparison of the algorithm (H-MODE) proposed in this study with those reported in the literature, the algorithm parameters were set to be the same as those reported in the literature.2 The population size was fixed at 100. The value of the crossover constant was fixed at 0.9, whereas the scaling factor was obtained randomly during each call of the recombination operation. The scaling factor varied between 0-1. However, the maximum number of generations was fixed at 30 (in case study 1) instead of 250 (as reported in the literature2,29) because hybrid-MODE algorithm converges to the Pareto front within less than 30 generations, as determined in our preliminary experimental runs. Figure 3a shows the convergence metric profiles, as a function of number of generations, of the MODE III algorithm and the proposed hybrid-MODE algorithm. In the H-MODE algorithm, the sequential simplex method is applied to the trial vector obtained using MODE (see Figure 1), and this hybridization results in a faster convergence to the true Pareto front, as shown in Figure 3a. The values of the control parameters used in the present study are listed in Table 2. The smaller value of the convergence metric indicates that the solutions are converged toward the true Pareto front. To calculate the convergence metric, we considered 500 uniformly diverse population points in the true Pareto front for the KUR test problem. Figure 3a shows that the hybrid-MODE algorithm converged to the true Pareto front within less than 25 generations. The algorithm further continued so as to give comparative performances of the hybrid-MODE and MODE III algorithms. The MODE III algorithm converged to the true Pareto front to the closest possible value, but computationally, it is slightly more expensive than the hybrid-MODE algorithm. Table 3 shows an analysis of the convergence metric (γ), the divergence metric

Figure 3. Comparison of convergence metric as a function of (a) generation numbers and (b) number of function evaluations (NFE). Table 2. Values of Control Parameters Used in the Reference Run of the Present Study parameter population size (NP) number of generations (MAXGEN) scaling factor (F) crossover rate (CR) reflection factor (R) contraction factor (β) expansion factor (γ)

MODE III

hybrid MODE

100 250

100 250

randomly obtained between 0 and 1 0.9 -

randomly obtained between 0 and 1 0.9 1 0.5 2

(∆), and the number of population points (NPP) converged to the true Pareto front in the specified numbers of generations (as given in the footnote of Table 4). These metrics were calculated for 10 continuous runs (to keep the same basis for comparison with other algorithms reported in the literature, as shown in Table 4 in the later part of this study), and the results are presented in Table 3. During all 10 runs, both the MODE III and hybrid-MODE algorithms converged to the true Pareto front (as it shows the smallest and most consistent value of the convergence metric). The number of function evaluations (NFE) is another important parameter that needs to be considered during the


11121

Table 3. Hybrid-MODE and MODE III Analyses of the Convergence Metric (γ) and the Divergence Metric (∆) and Number of Points Converged to the Pareto (NPP) Front in 10 Continuous Runs run no. 1

2

3

4

5

6

7

8

9

10

0.0039 0.6547 66

0.0044 0.6318 71

0.0033 0.7638 80

0.0038 0.6056 83

0.0035 0.6718 65

0.0028 0.6436 78

0.0028 0.6108 80

0.0024 0.6665 81

a

Hybrid-MODE Algorithm (Generation 300) γ ∆ NPP

0.0034 0.6207 76

0.0039 0.6540 75

0.0032 0.6935 83

γ ∆ NPP

0.0024 0.6438 66

0.0022 0.6205 78

0.0035 0.6718 78

0.0031 0.6823 82

0.0041 0.7354 79

0.0037 0.7165 80

MODE III Algorithm (Generation 300)

a

0.0029 0.6202 76

0.0030 0.6864 81

0.0029 0.5857 75

Results obtained in 30 generations remained consistent thereafter until 300 generations.

Table 4. Mean and Variance Values of the Convergence (γ) and Divergence (∆) Metrics Obtained Using Various Algorithms for the KUR Test Problema algorithm

no. of generations

mean

variance

0.028964 0.028951 0.049077 0.057323 0.03837 0.003853 0.005321

0.000018 0.000016 0.000081 0.011989 0.00057 5.45 × 10-7 7.2244 × 10-7

Convergence Metric (γ) NSGA II real-codedb NSGA II binary-codedb SPEAb PAESb MODEc MODE IIIc hybrid MODEd

250 250 250 250 300 300 300

Divergence Metric (∆) NSGA II real-codedb NSGA II binary-codedb SPEAb PAESb MODEc MODE IIIc hybrid MODEd

250 250 250 250 300 300 300

0.411477 0.442195 0.880424 1.079838 0.82097 0.642162 0.603399

0.000992 0.001498 0.009066 0.013772 0.0053 0.001049 0.00172178

a Best result in bold font; second-best in italic font. b Results reported in ref 2. c Results obtained in the present study. d Results obtained in 30 generations remained consistent thereafter until 300 generations.

design of any algorithm. Because the hybrid method involves the combination of evolutionary and local search approaches in series, the overall number of function evaluations is high. However, it is important to note that the hybrid strategy approaches the true Pareto front in fewer function evaluations (with a lower NFE), as shown in Figure 3b. The hybrid strategy took less than 5000 function evaluations to converge to the true Pareto front, whereas MODE III approached the true Pareto front in 7500 function evaluations. Table 4 shows a comparison of the MODE, MODE III, and hybrid-MODE algorithms with other well-known algorithms. The results for the other algorithms, namely, NSGA-II (both binary- and real-coded versions), SPEA, and PAES, were taken from the literature.2 The simulation runs of the existing study were carried out under identical conditions of parameters (NP ) 100 and CR ) 0.9) as given in the literature.2 However, the value of the number of generations (shown in Table 4) was varied based on observations made with reference to the number of function evaluations, as well as on the output observed in Figure 3a,b (rate of convergence). Table 4 shows that both the strategies of MODE converged to the true Pareto front with the smallest values of convergence metric as compared to the other algorithms. The results of the hybrid-MODE and MODE III algorithms in terms of the divergence metric are shown in Figure 4. During the initial generations, the divergence metric was found to increase. This is due to the nonconvergence of the algorithm to the true Pareto front. Once the algorithm converged to the true Pareto front, the divergence metric also decreased to its lowest

Figure 4. Comparison of the curve of the divergence metric with 300 generations for the KUR test problem using the hybrid-MODE and MODE III algorithms.

Figure 5. Comparison of true Pareto optimal solutions with those obtained using the hybrid-MODE and MODE III algorithms for the KUR test problem.

possible value within the next few generations. Figure 5 shows the Pareto fronts obtained for the KUR test problem using the two strategies of the MODE algorithm. Both the MODE III and hybrid-MODE algorithms converged to the true Pareto front and covered the entire discontinuous region on the true Pareto front.

11122


Table 4 indicates that the real- and binary-coded algorithms of NSGA-II resulted in convergence metric values of 0.028964 and 0.028951, respectively, whereas the MODE III and hybridMODE algortithms resulted in convergence metric values of 0.003853 and 0.005321, respectively. Thus, there is an orderof-magnitude difference in the convergence metric values obtained using NSGA-II variants and MODE strategies. In the case of MOO studies, although both convergence and diversity issues are equally important, it is also a fact that having a better diversity of solutions is insignificant if the algorithm does not converge to the true Pareto front. As per the definition (eq 9), the divergence metric would approach a value of zero (the best possible value) if the sum of all Euclidian distances matched the average of all Euclidian distances and if dei ) 0. If an algorithm does not converge to the true Pareto front and if the average of Euclidean distance is high, then the algorithm could also give a very low value of divergence metric. Thus, there exists a potential tradeoff between the stated goals of any MOO study (i.e., optimizing the convergence and divergence metrics). Hence, it is important for any algorithm to aim first for good convergence and then for satisfactory divergence rather than looking for good divergence and then satisfactory convergence. Even searching for equally good or compromised solutions (considering convergence and divergence metrics as goals) is not important, as it could keep the obtained nondominated solutions away from the true Pareto front. In all strategies of the MODE algorithm, there is a competition between the child (trial) and parent (target) vectors in each generation, and only the winner between the two moves to the next generation (unlike in other MOOs such as NSGA and its variants); therefore, MODE strategies result in very good and faster convergence to the true Pareto front. With the successful results obtained from application of the hybrid-MODE algorithm to the test problem (KUR), we applied this algorithm to an industrial application, namely, multiobjective optimization of an adiabatic styrene reactor, the results of which are discussed in the next section. 4.2. Case Study 2: Multiobjective Optimization of Adiabatic Styrene Reactor. A set of ordinary differential equations (ODEs) (as given in eqs A8-A12 of Appendix A) describing the reaction scheme were integrated and simulated using the ODE45 subroutine of MATLAB (7.0) library. ODE45 uses a fourth-order Runge-Kutta method to integrate the ODEs. The design and operating conditions, as well as the thermodynamic and kinetic data, related to the reaction scheme (as given by eqs 3-8) of the styrene reactor considered in this study were taken from the literature.8,25 The simulated model output was compared with the industrial data, and the results are tabulated in Appendix A. 4.2.1. Case 1: Simultaneous Maximization of SST and YST. The first case considers simultaneous maximization of selectivity and yield. Four decision variables (namely, TEB, P, SOR, and F0EB) are used. As discussed in section 4.1, the hybridMODE algorithm is much faster in terms of convergence for the benchmark test problem. A similar trend is observed in the case of the industrial styrene problem as well. In Figure 6, the conflicting objectives (selectivity and yield) are plotted at various generations, as indicated. The penalty function method was used to handle the constraints with a very high value of weights (104) in order to ensure that the constraint-dominated solutions were removed from the Pareto solutions. The effect of the penalty function parameter on the Pareto front was studied and is discussed in section 4.2.2.1. The numbers of points violating the constraints at several generations as obtained using the two strategies of MODE algorithms are reported in Table 5. In both algorithms, all of

Figure 6. Convergence of the hybrid-MODE algorithm toward the Pareto front for case 1. Table 5. Number of Points Violating the Constraints after Several Generations generation

hybrid MODE

MODE III

1 3 5 300

15 0 0 0

33 3 0 0

the initial population solutions are evaluated first and then sent to the generation loop. The values given in Table 5 correspond to the output obtained after the specified numbers of generations. As the initial population was generated randomly, the number of solutions violating the constraint in the initial population can vary depending on the random population generated. The hybridMODE algorithm, because of its acceleration phase (sequential simplex method coupled with evolutionary MODE algorithm), was able to remove the constraint-dominated solutions more quickly, thus leading toward the true Pareto front. The MODE III algorithm had 0 constraint-violated points in generation 5, whereas the hybrid-MODE algorithm achieved this stage during generation 3. It might be possible that, during the search process, a new point might be encountered that has a violation of constraints. The selection strategy of both MODE algorithms does not allow such solution (point) to enter into the current population list. However, if both the current point and the newly obtained point are constraint-violated points, then the better of the two points gets a place in the population. Thus, the selection strategy of both the algorithms is designed in such a way that even a not-very-efficient constraint handling technique such as the penalty approach works very well, ensuring that the constraintdominated solutions are removed from the initial generations and are not allowed to enter during later generations, unless they are better than the current constraint-dominant solutions. Table 5 and Figures 6 and 7 clearly show that the hybridMODE algorithm is able to converge toward the Pareto front more quickly than the MODE III algorithm. The inference of Figures 6 and 7 can also be correlated with the help of the set 0 ). The decision of decision variables (i.e., TEB, P, SOR, and FEB variables corresponding to the Pareto solutions in Figures 6 and 7 are shown against one of the objective functions in Figure 8a-h. Parts a and b of Figure 8 show the effect of the temperature of ethyl benzene on the objective functions (namely, selectivity and yield, respectively). Because the main reaction is reversible and endothermic in nature, a high temperature favors the rate of the forward reaction. This is apparent from the plot of yield versus temperature (Figure 8b). However, at higher temperature, side products such as toluene and benzene are also formed, thus reducing the selectivity value (Figure 8a).


11123

population sizes of 100 and 200, where the convergence is equally good for both cases. The Pareto front is able to cover a wider range (with scattered points) when the initial population size is kept at 200. However, with an initial population size of 100, the convergence of the algorithm to the Pareto front is good with nearly uniform diversity. Thus, the population size was chosen as 100 for the rest of the experimental runs in this study. 4.2.1.2. Effect of Step Size Used for Obtaining a New Neighborhood Solution. The hybrid algorithm uses a local search of solutions by creating a neighborhood point near the current point. The following mapping rule is used to generate a new neighborhood point new neighborhood solution ) current solution ( rand(0, 1) × desired step size Figure 7. Convergence of the MODE III algorithm toward the Pareto front for case 1.

Thus, a clear conflicting behavior is observed in the two objectives, which is attributed to the dominance of temperature on the objectives, and hence, it is termed a dominant variable. To select a particular value of yield, the user has to sacrifice in terms of selectivity and vice versa. If an overly high value of yield is selected, then the corresponding value of selectivity would be on the lower side, thus increasing the costs incurred for the separation of side products from the main product (i.e., styrene). The operating pressure also affects the desired values of the objective functions. According to Le Chatelier’s principle, low pressure favors the formation of the main product. By lowering the value of the operating pressure, the selectivity is increased, whereas the yield is decreased (Figure 8c,d). The feed flow rate of ethyl benzene approaches its lower bound because a lower flow rate is also responsible for generating a relatively low pressure. If the initial ethyl benzene flow rate is maintained at a lower value, then the mixture of ethyl benzene and steam would produce a relatively higher temperature (as per the energy balance of mixing streams). A high temperature and low pressure are favored at lower initial flow rate of ethyl benzene, and therefore, all of the points corresponding to the Pareto optimal solutions correspond to a lower initial ethyl benzene flow rate (Figure 8e,f). The steam-to-reactant ratio (SOR) also controls the desired objectives. A higher SOR value favors a high value of yield and a low value of selectivity. Figure 8a-h also shows comparisons of the decision variables corresponding to the Pareto optimal solutions obtained using the MODE III and hybrid-MODE algorithms. A comparatively better trend of the decision variables is observed in the case of the hybridMODE algorithm, whereas the decision variables are slightly scattered in the case of the MODE III algorithm. Figure 9 shows the Pareto optimal solutions obtained using both the hybrid-MODE and MODE III algorithms after 300 generations. The solutions obtained using the hybrid-MODE algorithm converged completely to the Pareto front, whereas those obtained using the MODE III algorithm are scattered at both the ends. The percentages of initial population points converged to the Pareto front after 300 generations as obtained using the hybrid-MODE and MODE III algorithms were 97% and 72%, respectively. 4.2.1.1. Effect of Initial Number of Population Points. The initial number of population points can play an important role in determining the distribution of solutions on the Pareto front. However, as the number of initial population points increases, the total number of function evaluations also increases. Thus, an optimum initial population size is necessary in any MOO algorithm. Figure 10 shows the Pareto fronts obtained using initial

(14)

The step size used in eq 14 also affects the quality of the Pareto front. We studied this aspect by considering various values of step size. Three different step sizes were considered for each of the variables depending on the number of digits in the variable in the current solution. For example, the upper bounds of the four decision variables considered in this study, namely, TEB, 0 , are 800 K, 2.63 bar, 20, and 40.56 kmol/h, P, SOR, and FEB respectively. Therefore, for these variables, the three sets of step sizes were 100, 1, 10, and 10 (step size 1); 10, 0.1, 1, and 1 (step size 2); and 1, 0.01, 0.1, and 0.1 (step size 3). The effect of the step size on the Pareto front was studied and is shown in Figure 11a. For step size sets 1, 2, and 3, the percentages of the initial population points converging to the final Pareto front were 86%, 96%, and 93%, respectively (Figure 11a). Because all of the results with these three sets of step sizes converged to the same front and because it is difficult to distinguish the overlapping points on the Pareto front, these results are replotted in Figure 11b, using vertical displacements in the values by +2 for the results obtained with step size 2 and by +3 for the results obtained with step size 3. Also, the diversity of the Pareto front was better when step size 2 was used. With step size 1, the location of neighborhood solution can be far from the current solution (as per eq 14), which might not give a better solution when the local search method is used. However, when step size 3 is used, the solutions are crowded on the Pareto front, as shown in Figure 11a,b. With step size 2, a well-diversified Pareto front is obtained with a maximum number of solutions on the Pareto front. Therefore, in the remaining simulation runs, step size 2 was used to create a new neighborhood point. Similar results were obtained when the hybrid-MODE algorithm was tested on several other test problems. 4.2.2. Case 2: Simultaneous Maximization of FST and SST. Figure 12 shows the Pareto optimal solutions obtained after 300 generations using the hybrid-MODE and MODE III algorithms. Twenty points were selected at random from Figure 12 and are listed along with their respective decision variables in Table 6. Table 6 shows that the ranges covered by the MODE III and hybrid algorithms in terms of objective 1 (i.e., FST) are 6.60-12.36 and 7.56-15.84 kmol/h, respectively, and in terms of objective 2 (i.e., SST), they are 87.18-95.5% and 87.22-95.41%, respectively. The hybrid-MODE algorithm produced a wider range in terms of both objectives when compared to the range of objectives obtained using MODE III and the industrial operating point (taken from the literature).27,28 The MODE III algorithm resulted in lower values of FST as compared to the hybrid-MODE algorithm because MODE III approached the lower bound of the decision variable 0 0 . If the FEB decision variable values are compared for the FEB MODE III and hybrid-MODE algorithms, then it is observed that the hybrid-MODE algorithm approached the upper bound (i.e.,

11124


0 Figure 8. (a-h) Decision variables (T, P, FEB , and SOR) plotted against one of the objective functions.

40.56 kmol/h) whereas the MODE III algorithm approached the lower bound (i.e., 27.56 kmol/h). In earlier studies by Yee et al.25 and Babu et al.,8 it was observed that the multiobjective optimization study of maximization of FST and SST was independent of this decision 0 0 variable (FEB ) because the lower bound of FEB was captured in

these studies. In contrast, in this study, it is seen that only a robust algorithm can search the entire search space and provide a possibly true Pareto optimal front for an industrial problem. A compromise (corresponding to the set of optimum nondominated solutions) with other conflicting variables (such as T, P, and SOR) is required to obtain the suitable Pareto front with


11125

Figure 9. Pareto fronts obtained using the hybrid-MODE and MODE III algorithms after 300 generations.

Figure 10. Effect of initial population size (NP) on convergence of the hybrid-MODE algorithm after 300 generations.

higher values of the initial ethyl benzene flow rate. If such a compromise is not met along with the initial ethyl benzene flow rate (one of the decision variables), then a nonoptimal solution or a set of nondominated solutions might be produced away from the true Pareto optimal front. Hence, the algorithm might result in a local Pareto front. When a higher initial ethyl benzene flow rate is obtained (in the case of the hybrid algorithm), the values of the corresponding decision variables such as T and SOR are also high. This is because the higher temperature of the mixed stream (steam and initial ethyl benzene) is responsible for the higher productivity. It is important for any efficient algorithm to visit each and every corner of the multidimensional search space so that a true Pareto front is obtained. The MODE III algorithm in this study (see Figure 12) and the MODE and NSG algorithms in earlier studies8,25 converged to the lower bound of the initial ethyl benzene flow rate and, hence, resulted in local Pareto fronts. Table 6 also shows the conflicting scenario, wherein the value of SST decreases while the calculated value of YST increases as the value of FST increases. The ranges covered by the MODE III and hybrid-MODE algorithms in terms of the calculated value of YST are 21.38-42.02% and 20.41-37.51%, respectively. A profit function can be used as a further guideline for choosing the appropriate nondominated point from the Pareto front. The profit function considers the costs of major feed materials (costs of ethyl benzene and steam) and final valued products (such as styrene, benzene, and toluene). For this purpose, a simplified profit function8,25 is defined as profit ) revenue generated by styrene and byproducts raw material cost (15)

Figure 11. (a) Effect of the step size used in obtaining a neighborhood solution on the Pareto optimal solutions obtained using the hybrid-MODE algorithm. (b) Results of part a replotted with vertical shifts in the value of the ordinate by +2 for the step size 2 data points and by +3 for the step size 3 data points.

Equation 15 in mathematical form can be written as profit ) FSTHST + FBZHBZ + FTOLHTOL (F0EB - FEB)HEB - FH2OHH2O (16) The costs of styrene, ethyl benzene, benzene, and toluene were based on the prices published in an online purchasing magazine.30 However, we considered the older cost values (as of April 2001) in this study for comparison with the previously obtained results of Yee et al.25 It was observed that the profit does not include the cost of separating byproducts from styrene (eq 16). Hence, the results related to the profit function analysis as per eq 16 with respect to the objective function values and key decision variables are not provided here. The results of the current study are compared with those obtained with the NSGA published in the literature.25 Figure 13 shows a comparison of the Pareto fronts obtained using the MODE III, hybrid-MODE, and NSG algorithms, along with an industrial operating point. From the results obtained with each of the algorithms, some data points were selected and are marked with boxes (0) in Figure 13. These selected data points were analyzed further and are listed in Table 7. Table 7 shows the selected data points from the Pareto front of Figure 13 (points A, B, C, and G using the MODE III algorithm; points A′, B′, C′, and G′ using the hybrid-MODE algorithm; and points D, E, and F using the NSGA). Figure 13 shows that the hybrid-MODE algorithm resulted in better nondominated solutions than the MODE III and NSG algorithms when the same sets of decision variables were used. A better Pareto front (as mentioned above) would

11126


Figure 12. Pareto optimal solutions obtained after 300 generations using the hybrid-MODE and MODE III algorithms. Table 6. Values of Objective Functions, Decision Variables, and Profit Functions for Selected Data Points of Figures 12 and 13 MODE III a ST

improved hybrid MODE a ST

no.

FST (kmol/h)

SST (%)

Y (%)

T (K)

P (bar)

(kmol/h)

SOR

profit ($/h)

FST (kmol/h)

SST (%)

Y (%)

T (K)

P (bar)

0 FEB (kmol/h)

SOR

profit ($/h)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20

6.60 7.52 7.80 8.39 8.50 8.72 9.10 9.29 9.37 9.57 9.65 9.71 9.93 10.10 10.17 10.31 10.61 11.17 11.67 12.36

95.50 94.78 94.49 94.03 93.70 93.38 92.81 92.74 92.47 92.24 92.04 92.03 91.75 91.40 91.33 91.01 90.36 89.79 88.67 87.18

21.38 24.38 25.49 26.66 28.08 29.08 30.54 30.66 31.37 32.04 32.60 32.60 33.29 34.12 34.30 34.95 35.99 37.30 39.33 42.02

702.13 743.52 597.88 662.48 590.35 596.11 676.66 580.08 635.9 645.45 625.49 638.46 557.29 635.89 584.44 711.99 742.34 697.68 721.25 686.97

1.025 1.199 1.696 1.687 1.955 2.004 1.806 2.154 2.061 1.99 2.137 2.009 2.353 2.068 2.374 1.823 1.599 2.071 2.133 2.476

27.77 28.12 27.97 28.98 27.88 27.71 27.62 28.13 27.75 27.77 27.56 27.75 27.82 27.63 27.71 27.60 27.63 28.17 27.99 27.82

9.470 7.416 13.79 11.38 14.62 14.81 12.44 16.09 14.21 14.53 15.33 15.29 18.28 16.44 17.78 13.43 13.62 15.40 15.29 19.21

344.6 415.9 367.3 420.8 398.2 409.7 454.3 425.7 450.9 458.2 456.1 458.7 440.1 468.5 459.1 510.2 523.6 533.9 562.1 559.1

7.72 8.14 9.09 9.33 9.57 10.41 10.80 11.18 11.50 11.76 12.22 12.47 12.96 13.30 13.52 14.17 14.75 15.12 15.68 15.81

95.41 95.20 94.73 94.67 94.50 93.76 93.29 93.05 92.68 92.42 92.07 91.72 91.14 90.76 90.50 89.71 88.81 88.17 87.56 87.22

20.41 20.96 21.88 21.87 22.08 24.38 25.59 26.08 26.86 27.75 28.39 29.27 30.48 31.33 31.60 33.18 34.93 35.85 36.90 37.55

692.90 613.54 694.78 646.1 596.47 650.35 629.02 695.94 735.88 646.92 639.64 675.58 592.57 601.05 639.60 746.45 744.74 742.20 739.63 705.58

1.201 1.571 1.415 1.665 1.821 1.860 1.921 1.793 1.651 2.064 2.157 2.099 2.449 2.565 2.573 1.979 2.035 2.270 2.278 2.556

34.56 35.64 38.52 39.6 40.32 39.96 39.6 40.32 40.32 39.96 40.68 40.32 40.32 40.32 40.68 40.68 40.32 40.32 40.68 40.32

8.788 11.3904 8.7819 10.3246 12.1912 10.9558 12.6827 9.9042 8.5383 12.7333 13.3325 12.1561 16.23 15.9961 14.2198 11.2624 12.8457 12.6923 14.207 16.1799

393.3 380.0 458.1 446.0 429.5 495.5 493.8 551.9 589.0 544.0 557.1 589.6 556.7 578.8 615.0 692.3 701.3 723.0 728.4 708.2

a

0 FEB

Calculated value of third objective function during two-objective optimization study.

mean that at least one member of the current set of solutions has a better value of at least one objective than all solutions of the other set and that none of the members of the current Pareto front are worse than the other set of solutions. Pareto solutions obtained using the hybrid-MODE algorithm covered a wider

range on the Pareto front as compared to the nondominated solutions obtained when the NSGA was used. Figure 13 also shows that both the MODE III and NSG algorithms converged to a local Pareto front compared to the Pareto front obtained using the hybrid-MODE algorithm. The Pareto front obtained

Figure 13. Comparison of Pareto fronts obtained using the MODE III, hybrid-MODE, and NSG algorithms.


11127

Table 7. Comparison of MODE III, Hybrid-MODE, and NSG Algorithms for Selected Data Points of Figure 13 MODE III (present study)

hybrid MODE (present study)

NSGA25

parameter

industrial value

A

B

C

G

A′

B′

C′

G′

D

E

F

FST (kmol/h) SST (%) YST (%) TEB P (bar) SOR 0 FEB FBZ FTOL(kmol/h) profit ($/h)

14.9 85.15 40.30 800.00 2.4 12.29 36.87 1.37 1.20 651

6.608 95.50 21.38 702.1 1.02 9.47 27.77 0.30 0.96 344

9.93 91.75 33.29 557.29 2.35 18.28 27.82 0.60 1.21 440

12.36 87.18 42.02 686.9 2.47 19.21 27.82 1.18 1.52 559

11.91 87.88 40.60 780.7 1.69 14.35 27.7 1.18 1.36 584

7.72 95.41 20.41 692.9 1.20 8.78 34.56 0.32 1.00 393

12.47 91.72 29.27 675.5 2.09 12.15 40.32 0.70 1.34 589

15.81 87.22 37.55 705.5 2.55 16.17 40.32 1.47 1.73 708

15.68 87.56 37.55 739.5 2.55 14.20 40.68 1.43 1.68 728

8.45 94.86 21.88 651.5 1.84 10.74 38.63 0.24 0.22 330

11.04 92.69 27.38 713.5 1.84 10.69 40.31 0.50 0.37 467

15.01 87.51 35.26 799.9 2.03 10.94 40.47 1.25 0.90 673

using the hybrid-MODE algorithm is far superior to the local Pareto fronts obtained using the NSG and MODE III algorithms. The hybrid-MODE algorithm is able to give better results in terms of both objectives. Table 7 also provides a comparison of the industrial operating point with the selected points obtained using different algorithms with various parameters. Moving from point A to point C (and similarly from point A′ to point C′ and from point D to point F), the values of both FST and profit increases. The values of FST and profit increase in the case of the MODE III algorithm output primarily because of the increase in SOR (from 9.47 to 19.21). The MODE III algorithm results in a local Pareto solution because a point (henceforth referred to as a chromosome, in accordance with evolutionary terminology) has a higher value of TEB (chromosome A) and a lower value of SOR and vice versa (chromosome C). As per the energy balance, this combination of decision variables (TEB and SOR) is responsible for generating the temperature of the mixed stream entering the reactor (Tmix1). In the case of the hybrid-MODE algorithm, a high value of either TEB or SOR results in a higher value of Tmix1 (e.g., chromosome A′ versus D gives TEB ) 692.9 K and SOR ) 8.78 by the hybrid-MODE algorithm versus TEB ) 651.5 K and SOR ) 10.74 by the NSGA, and chromosome C′ versus F gives TEB ) 705.5 K and SOR ) 16.71 by the hybrid-MODE algorithm against TEB ) 799.9 K and SOR ) 10.94 by the NSGA). Tmix1 has a strong opposing effect on FST and SST. A high value of Tmix1 maximizes the flow rate of styrene (point C′), whereas a low Tmix1 maximizes selectivity (point A′). This is due to the reversible endothermic nature of the main reaction (eq 3). The NSGA results in a local Pareto front because SOR approachs a lower bound. Thus, it is necessary for an efficient algorithm to have a widespread search of population points in all dimensions. Because the hybridMODE algorithm approaches the Pareto optimum values of

all decision variables, it results in a wider and better Pareto front as compared to the NSGA and the MODE III algorithm. Point G (profit ) $584/h) and G′ (profit ) $728/h) (in Figure 13 and Table 7) are the points having maximum profit values obtained using the MODE III and hybrid-MODE algorithms respectively. Point G′ has the maximum value of profit among the nondominated solutions obtained using different algorithms (i.e., the hybrid-MODE, MODE III, and NSG algorithms) and an industrial operating point as shown in Table 7. 4.2.2.1. Effect of a Penalty Parameter on the Pareto Front. The penalty function approach is used to handle the constraints. Three different weights of penalty parameters (i.e., 104, 108, and 1012) were used in the present work to study the effect of the penalty parameter on the MOO of adiabatic styrene reactor, and the Pareto fronts obtained are shown in Figure 14. The algorithm converged to the same Pareto front, but the diversity of solutions differed for the different penalty parameter values. However, the diversity of solutions was found to be independent of the penalty parameter for most of the range, except that no specific trend is observed toward the extreme ends of the Pareto front. Considering the same convergence rate, the penalty parameter value of 104 was used in this study. It should be noted that a population size of 100 was used for these simulations consistently throughout this study. However, with a population size of 200 and with a penalty parameter value of 104, the Pareto front covered the same range as with a population size of 100 and with a penalty parameter value of 1012. The diversity of the obtained solutions is one area of the hybrid-MODE algorithm that needs to be further improved. Considering the same convergence rate and computational complexity (with high value

Figure 14. Effect of the penalty parameter and population size on the Pareto optimal solutions obtained using the hybrid-MODE algorithm.

11128



11129

Figure 15. (a) Pareto optimal solutions obtained after 300 generations using the hybrid-MODE and MODE III algorithms and (b-i) corresponding decision variables. Table 8. Summary of Optimization Case Studies Analyzed Using MODE III, Hybrid-MODE, and NSG Algorithms MODE III (present study)

hybrid MODE (present study)

NSGA25

parameter

industrial value

FST vs SST

SST vs YST

FST vs YST

three-obj optim study

FST vs SST

SST vs YST

FST vs YST


FST vs SST

SST vs YST

FST vs YST


FST (kmol/h) SST (%) YST (%) TEB P (bar) SOR 0 FEB FBZ FTOL(kmol/h) profit ($/h)

14.9 85.15 40.30 800.00 2.4 12.29 36.87 1.37 1.20 651

11.91 87.88 40.60 780.7 1.69 14.35 27.7 1.18 1.36 584

12.38 86.80 42.23 762.4 2.13 15.07 27.74 1.18 1.58 601

16.73 84.75 39.84 771.51 2.62 13.49 40.31 1.67 2.20 792

14.47 88.24 36.52 791.61 1.75 11.35 37.81 1.29 1.53 716

15.68 87.56 37.55 739.5 2.55 14.20 40.68 1.43 1.68 728

12.26 87.35 41.83 667.53 2.55 19.91 27.72 1.15 1.51 547

16.64 84.89 39.97 767.15 2.62 13.86 39.96 1.67 2.16 783

16.41 86.14 39.10 732.60 2.54 16.03 40.25 1.67 1.84 741

15.01 87.51 35.26 799.98 2.03 10.94 40.47 1.25 0.90 673

10.89 88.83 39.47 787.25 1.76 12.54 27.58 0.85 0.82 482

14.48 84.70 40.80 800.00 2.45 12.47 35.49 1.36 1.25 637

13.09 88.18 38.23 776.13 1.99 13.16 34.24 1.10 0.66 565

of NP), the penalty parameter value of 104 and population size of 100 were used in this study. 4.2.3. Case 3: Simultaneous Maximization of FST and YST. The Pareto front obtained for case 3 (simultaneous maximization of FST and YST) and the corresponding decision variables are shown in Figure 15a-i. Unlike the results obtained 0 approached for cases 1 and 2 (where the decision variable FEB either the lower or upper bound), for case 3, the decision variable F0EB was also equally important in producing the Pareto solutions. While maximizing FST and YST simultaneously, the decision variable P remained practically constant and acquired the upper bound. The conflicting variables observed in this study (which are responsible for producing the Pareto solutions) are T, SOR, and 0 FEB . These results show the ability of the existing algorithm to produce more valuable and practical results that are important to the plant engineer. Figure 15a shows that the industrial point lies below the Pareto front and the nondominated solutions are present on either side of the industrial point, thus offering a wide range and choice to the decision maker. However, in the study of Yee et al.,25 the industrial point (for maximization of FST and YST) lay at one extreme (left) end of the obtained nondominated solutions (Pareto front in their work), giving a limited choice for the decision maker. Parts b, c, and f-i of Figure 15 show that the temperature of the industrial operating point is higher than the optimum Pareto temperatures and the SOR of the industrial point is lower than the optimum Pareto SOR obtained in this study. However, even though the value of F0EB (decision variable) is the same as the Pareto values, the Pareto front is away from the industrial point because of the conflicting effects of temperature and SOR on the said objectives. Therefore, an optimum and judicious compromise among all of the conflicting decision variables is important in such industrial multiobjective optimization studies.

4.2.4. Case 4: Simultaneous Maximization of FST, SST, and YST. A three-objective optimization study was carried out by simultaneously maximizing FST, SST, and YST. The point corresponding to the best profit was taken from each of the case studies. The decision variables, objectives, concentrations of unwanted byproducts, and profits corresponding to these selected points along with those for the industrial point are listed in Table 8. The profit obtained for the SST versus YST objectives (case 4) (using the MODE III, hybrid-MODE, and NSG algorithms) is lower than the profit of the industrial operating point. The profit strongly depends on the amount of styrene produced. Because the productivity of styrene is not one of the objectives in SST versus YST studies, the obtained best profit is lower than that of the industrial operating point. However, for the remaining cases (cases 1, 2, and 4) (both twoobjective and three-objective optimization studies), the hybridMODE algorithm resulted in a higher profit than was obtained with the industrial operating point or with MODE III and NSG algorithms. Unlike the hybrid-MODE and MODE III algorithms results, the NSGA results show that the profit obtained during a three-objective optimization study is less than the industrial operating profit. An evolutionary multiobjective optimization algorithm usually results in a set of solutions (probably Pareto front) in a single simulation run. The sets of decision variables corresponding to these sets of solutions are also known. The decision maker can select any one particular solution from the Pareto front suitable to his/her requirements. Once the decision maker selects a point suitable to his/her needs from the objective space, the corresponding decision variables can be obtained from the decision variable space. It is expected that, if the mathematical model output matches the industrial data (along with the satisfaction of constraints), the selected set of decision variables, if adopted, can result in the desired values of objectives as selected by the decision maker.

11130


5. Conclusions Multiobjective optimization of the KUR test function (case study 1) and an adiabatic styrene reactor (case study 2) for the simultaneous maximization of three objectives (i.e., productivity, selectivity, and yield) was studied using the proposed hybridMODE and MODE III (improved strategy of MODE) algorithms. The following conclusions can be drawn from the results of the present study: The proposed hybrid-MODE algorithm is able to converge to the true Pareto front (for case study 1) and a better Pareto front than reported in the literature (for case study 2) at a much higher rate (in fewer function evaluations), thus acquiring both accuracy and speed. The performance metrics obtained in this study for case study 1 were found to be better in terms of convergence and comparable in terms of divergence compared to those obtained using other well-known MOO algorithms such as NSGA-II (both binary- and real-coded versions), SPEA, PAES, MODE, and MODE III. The diversity aspect of the hybrid-MODE algorithm needs further improvement. The hybrid-MODE algorithm is able to reach the global Pareto front (in case study 1) in less than 30 generations (NFE of less than 5000), whereas the evolutionary approach took 80 generations (NFE of 7500) to attain the same convergence. Thus, the major expected rationale of speeding up an algorithm is achieved by using a hybrid strategy of evolutionary algorithm. Pareto optimal fronts were obtained for four cases of case study 2. The results were discussed both qualitatively and quantitatively with respect to the conflicting effects of different decision variables (i.e., T, P, SOR, and F0EB). Sets of nondominated solutions obtained using several algorithms [NSGA, one of the strategies of MODE (MODE III), and hybrid MODE] were compared with respect to the key decision variables. The hybrid-MODE algorithm was found to be able to give a widespread Pareto front that was better than the nondominated solutions obtained using MODE III and NSG algorithms. The convergence of the Pareto front obtained in this study using the hybrid-MODE algorithm was found to be independent of several parameters such as the initial population points, the step size used for obtaining a new neighborhood solution in sequential simplex method, and the penalty parameter for constraint handling. However, even though the diversity of solutions was found to be independent of the penalty parameter for most of the range, no specific trend was observed toward the extreme ends of the Pareto front. Only the hybrid-MODE algorithm (from among the algorithms considered in this study) was able to capture a point having higher profit than the industrial operating profit in the case of three-objective optimization (case 4 of case study 2). The results of this study are important for plant engineers, who can use the analysis carried out in this study to improve the performance of working plants. Nomenclature A ) frequency factor; cross-sectional area of reactor (m2) C ) molar concentration (kmol/m3) Cp ) molar heat capacity [kJ/(kmol K)] CR ) crossover constant D ) diameter of the reactor (m) E ) activation energy (kJ/kmol) F ) molar flow rate (kmol/h); scaling factor G ) mass velocity of the gas mixture [kg/(m2 h)] K ) equilibrium rate constant (bar)

k ) reaction rate constant [kmol/(kg h barn)] l ) length of the reactor (m) N ) number NP ) number of populations used in differential evolution NPP ) number of population points converged to Pareto front NFE ) number of function evaluations p ) partial pressure (bar) P ) total pressure (bar) r ) reaction rate [kmol/(kg h)] R ) universal gas constant [8.314 kJ/(kmol K)] S ) selectivity (%) SOR ) steam-to-reactant (ethyl benzene) molar ratio T ) temperature (K) x ) conversion (%) Y ) yield (%) Greek Symbols γ ) convergence metric ∆ ) divergence metric ∆H ) heat of reaction (kJ/kmol) ε ) void fraction µ ) viscosity [kg/(m s)] F ) density (kg/m3) Subscripts/Superscripts o ) initial b ) bulk BZ ) benzene c ) crossover C ) catalyst CO ) carbon monoxide CO2 ) carbon dioxide EB ) ethyl benzene ETH ) ethylene G ) gas H ) heat-transfer area (m2) H2 ) hydrogen H2O ) steam i ) ith reaction j ) jth component MET ) methane mix1 ) mixture of fresh steam and feed ethyl benzene fed at reactor inlet p ) particle ST ) styrene t ) total TOL ) toluene

Appendix A The rates of reaction are as follows: r1 ) k1(pEB - pSTpH2 /KEB)

(A1)

r2 ) k2(pEB)

(A2)

r3 ) k3(pEBpH2)

(A3)

r4 ) k4(pH2OpETH0.5)

(A4)

r5 ) k5(pH2OpMET)

(A5)

r6 ) k6(pH2OpCO)(PT /T3)

(A6)

The rate constant ki of reaction i is expressed as


ki ) exp[Ai - (Ei /RT)]

(A7)

where Ai and Ei are the apparent frequency factor and the activation energy, respectively, of reaction i.27 The six material balance equations are given by FbAtri dxi ) 0 dl FEB

for i ) 1, 2, 3

(A8)

where xi is the fractional conversion of ethyl benzene in each of the three reactions, and FbAtri dxi ) 0 dl FH2O

for i ) 4, 5, 6

(A9)

where xi is the fractional conversion of steam in each of these three reactions. The Ergun equation31 is used to compute the pressure profiles

[

]

(1 - ε)G0 150(1 - ε)µG dP ) -1 × 10-5 + 1.75G0 dl Dp Dpε3FG (A10) The energy balance differential equation can be derived from the equation 10

∑FC i)1

i pi

6

dT +

∑ ∆H A F r dl ) UA j t b j

H

dl(T - Te)

j)1

(A11) Because the reactor is adiabatic, U A H dl (T - Te) ) 0, and the energy balance differential equation is given by 6

dT ) dl

∑ ∆H A F r

j t b j

j)1

(A12)

10

∑FC i)1

i pi

The heats of reactions are computed as functions of temperature from the relation27 ∆Hj ) aj + bjT

(A13)

The set of values of aj and bj in eq A13 were taken from Sheel and Crowe.27 The molar heat capacities (Cpi and Cpj) of the components are given as functions of temperature by the equations Cpi ) Ri + βiT + γiT2 (organic components) (A14) Cpj ) aj + bjT + cj /T2 (inorganic components) (A15) The set of the values of the constants of molar heat capacities Cpi for organic and inorganic components were taken from Smith and Van Ness32 and are available in our earlier publication.8 The additional value for the equilibrium constant of ethyl benzene was found from the literature.28 The equations for conversion and molar flow rate were taken from Elnashaie and Elshishini.28 Finally, Table A1 provides a comparison of the simulated model output with the industrial data.

11131

Table A1. Comparison of Simulation Runs from Present Study with Industrial Data27,28 parameter at the exit of reactor temperature (K) pressure (bar) conversion of ethyl benzene per pass (%) styrene flow rate (kmol/h) flow rate of ethyl benzene (kmol/h) benzene flow rate (kmol/h) toluene flow rate (kmol/h)

simulation results (present study)

industrial data

850.08 2.33 46.78

850 2.32 47.25

15.37 19.61 1.46 2.07

15.57 19.45 1.5 2.03

Note Added after ASAP Publication: The version of this paper that was published online November 4, 2009 had errors in the titles of Tables 6 and 7. The corrected version was reposted to the web December 9, 2009. Literature Cited (1) The Styrene Forum. http://www.styreneforum.org (accessed May 23, 2009). (2) Deb, K. Multi-ObjectiVe Optimization Using EVolutionary Algorithms; John Wiley & Sons: New York, 2001. (3) Deb, K.; Mitra, K.; Dewri, R.; Majumdar, S. Towards a better understanding of the epoxy polymerization process using multi-objective evolutionary computation. Chem. Eng. Sci. 2004, 59, 4261. (4) Gao, X.; Chen, B.; He, B.; Qiu, T.; Li, J.; Wang, C.; Zhang, L. Multi-objective optimization for the periodic operation of the naphtha pyrolysis process using a new parallel hybrid algorithm combining NSGAII with SQP. Comput. Chem. Eng. 2008, 32, 2801. (5) Agarwal, A.; Gupta, S. K. Jumping gene adaptations of NSGA-II and their use in the multi-objective optimal design of shell and tube heat exchangers. Chem. Eng. Res. Des. 2007, 86, 123. (6) Zitzler, E.; Laumanns, M.; Thiele, L. SPEA2: ImproVing the Strength Pareto EVolutionary Algorithm; Technical Report 103; Computer Engineering and Networks Laboratory (TIK), Swiss Federal Institute of Technology (ETH): Zurich, Switzerland, 2001. Available online at http://www. tik.ee.ethz.ch/sop/publicationListFiles/zlt2001a.pdf. (7) Knowles, J. D.; Corne, D. W. The Pareto Archived Evolution Strategy: A New Baseline Algorithm for Pareto Multiobjective Optimisation. In Proceedings of the 1999 Congress on EVolutionary Computation (CEC 1999); IEEE Press: Piscataway, NJ, 1999; pp 98-105. (8) Babu, B. V.; Chakole, P. G.; Mubeen, J. H. S. Multi-objective Differential Evolution (MODE) for Optimization of Adiabatic Styrene Reactor. Chem. Eng. Sci. 2005, 60, 4822. (9) Corne, D.; Dorigo, M.; Glover, F. New Ideas in Optimization; McGraw-Hill Publications: London, 1999. (10) Onwubolu, G. C.; Babu, B. V. New Optimization Techniques in Engineering; Springer-Verlag: Heidelberg, Germany, 2004. (11) Price, K. V.; Storn, R. Differential evolutionsA simple evolution strategy for fast optimization. Dr. Dobb’s J. 1997, 22, 18. (12) Babu, B. V.; Munawar, S. A. Differential Evolution Strategies for Optimal Design of Shell-and-Tube Heat Exchangers. Chem. Eng. Sci. 2007, 62, 3720. (13) Angira, R.; Babu, B. V. Optimization of Process Synthesis and Design Problems: A Modified Differential Evolution Approach. Chem. Eng. Sci. 2006, 61, 4707. (14) Babu, B. V.; Mubeen, J. H. S.; Chakole, P. G. Modeling, simulation, and optimization of wiped film poly ethylene terephthalate (PET) reactor using multi-objective differential evolution (MODE). Mater. Manuf. Processes 2007, 22, 541. (15) Reddy, J. M.; Kumar, N. D. Multiobjective Differential Evolution with Application to Reservoir System Optimization. J. Comput. CiVil Eng. 2007, 21, 136. (16) Babu, B. V.; Gujarathi, A. M.; Katla, P.; Laxmi, V. B. Strategies of Multi-Objective Differential Evolution (MODE) for Optimization of Adiabatic Styrene Reactor. In Proceedings of the International Conference on Emerging Mechanical Technology: Macro to Nano (EMTMN-2007); p 243. (17) Gujarathi, A. M.; Babu, B. V. Improved Multi-objective Differential Evolution (MODE) Approach for Purified Terephthalic Acid (PTA) Oxidation Process. Mater. Manuf. Processes 2009, 24, 303. (18) Chiou, J. P.; Wang, F. S. Hybrid method of evolutionary algorithms for static and dynamic optimization problems with application to a fedbatch fermentation process. Comput. Chem. Eng. 1999, 23, 1277.

11132


(19) Poloni, C.; Giurgevich, A.; Onesti, L.; Pediroda, V. Hybridization of a multi-objective genetic algorithm, a neural network and a classical optimizer for a complex design problem in fluid dynamics. Comput. Methods Appl. Mech. Eng. 2000, 186, 403. (20) Fonseca, C. M.; Fleming, P. J. An overview of evolutionary algorithms in multiobjective optimization. EVol. Comput. 1995, 3, 1. (21) Rao, S. S. Optimization Theory and Applications; Wiley Eastern: New Delhi, India, 1991. (22) Babu, B. V. Process Plant Simulation; Oxford University Press: New York, 2004. (23) Nelder, J. M.; Mead, R. A simplex method for function minimization. Comput. J. 1965, 7, 308. (24) Li, C. H.; Hubbell, O. S. Styrene. In Encyclopedia of Chemical Processing and Design; Mcketta, J. J., Weismantel, G. E., Eds.; Wiley: New York, 1982. (25) Yee, A. K. Y.; Ray, A. K.; Rangiah, G. P. Multi-objective optimization of industrial styrene reactor. Comput. Chem. Eng. 2003, 27, 111. (26) Huang, V. L.; Suganthan, P. N.; Qin, A. K.; Baskar, S. MultiobjectiVe Differential EVolution with External ArchiVe and Harmonic Distance-Based DiVersity Measure; Technical Report-MODE-2005; NTU:

Singapore, 2005.Available online at http://www3.ntu.edu.sg/home/ epnsugan/index_files/papers/Tech-Report-MODE-2005.pdf. (27) Sheel, J. G. P.; Crowe, C. M. Simulation and optimization of an existing ethyl benzene dehydrogenation reactor. Can. J. Chem. Eng. 1969, 47, 183. (28) Elnashaie, S. S. E. H.; Elshishini, S. S. Modeling, Simulation and Optimization of Industrial Fixed Bed Catalytic Reactors; Gordon and Breach Science Publishers: London, 1994. (29) Agrawal, S.; Dashora, Y.; Tiwari, M. K.; Son, Y. J. Interactive Particle Swarm: A Pareto-Adaptive Metaheuristic to Multi-objective Optimization. IEEE Trans. Syst., Man Cybernetics A 2008, 38, 258. (30) Purchasing.com. http://www.purchasing.com (accessed May 23, 2009). (31) Bird, R. B.; Stewart, W. E.; Lightfoot, E. N. Transport Phenomena; Wiley: New York, 1960. (32) Smith, J. M.; Van Ness, H. C. Introduction to Chemical Engineering Thermodynamics; McGraw-Hill: New York, 1975.

ReceiVed for reView July 4, 2009 ReVised manuscript receiVed September 7, 2009 Accepted October 22, 2009 IE901074K

A Hybrid Multiobjective Differential Evolution - ACS Publications

Recommend Documents