Biomimicking Altruistic Behavior of Honey Bees in Multi-objective

Ind. Eng. Chem. Res. 2009, 48, 9671–9685

9671

Biomimicking Altruistic Behavior of Honey Bees in Multi-objective Genetic Algorithm Manojkumar Ramteke+ and Santosh K. Gupta* Department of Chemical Engineering, Indian Institute of Technology Bombay, Powai, Mumbai 400 076, India

The altruistic behavior of honey bees provides an interesting contrast to natural selection in evolutionary biology. This is biomimicked in the framework of a multiobjective optimization algorithm, namely, genetic algorithm, GA, by exploiting the concept of elitism (preserving good parents). The effects of altruism and natural selection on the total fitness of the colony are compared. This basic algorithm is used for studying the evolution process. It is then modified to enhance the convergence rates of optimization problems and to simulate the carcinogenesis of cells using multiple queens, unlike in honeycombs, mimicking other species of hymenopterans, e.g., ants, wasps, etc. This algorithm provides a new approach for studying three problems, bee evolution, optimization, and cancer, and is used to understand conflicts in animal behavior, increase the speed of convergence of optimization problems, and for an improved understanding of the causes of cancer. 1. Introduction Biological phenomena are mimicked in various modern-day technologies. The field of optimization is no exception. Many recent optimization algorithms like simulated annealing (SA),1 ant colony optimization,2 etc. have been inspired by nature. A similar inspiration from evolutionary biology for solving optimization problems led to genetic algorithm (GA).3 A further improved multiobjective version of GA is the binary-coded elitist4 nondominated sorting genetic algorithm with the jumping gene5 adaptation,6 NSGA-II-aJG. In this algorithm, a population of several (Np) parent strings (called chromosomes), each comprising several binary numbers, is generated randomly. These represent values of the decision variables. The binaries associated with each decision variable are then coded/mapped into decimal values so that they all lie within their bounds. The decoded values of each variable are used in the model equations to calculate the objective functions in a multiobjective optimization problem. The chromosomes are then classified into “ranks” or “fronts” based on the concept of nondomination.4,7 Also, the Euclidean distance in the objective-function space between nearest neighbors, referred to as the “crowding distance”, is evaluated for all the chromosomes (in each front). Two chromosomes are picked randomly from the population of parent chromosomes, and the better of the two, based on the rank and the crowding distance, is copied into a mating pool. This is known as tournament selection. Np chromosomes are, thus, copied for subsequent operation. It may be mentioned that this procedure may lead to some inferior chromosomes also finding their way into the mating pool (this is a characteristic of GA, since it is possible that a poor chromosome may give a good daughter). These are the better “parents”. Two chromosomes are picked randomly from the Np in the mating pool to undergo crossover,4,7 mutation,4,7 jumping gene6 operation, etc. to produce two (initial) daughter chromosomes. The process of crossover among any two randomly selected strings is referred to as natural selection. The Np daughters so produced are mixed with the Np better parents. These 2Np chromosomes are reranked (using nondominance), and the best Np (final) daughter chro* To whom correspondence should be addressed. E-mail: sk.gupta@ che.iitb.ac.in. Tel.: 91-22-2576 7256. Fax: 91-22-2572 6895. On leave from IIT Kanpur 208016. + Department of Chemical Engineering, IIT Kanpur 208016.

mosomes are selected from these. The process of picking up Np (final) daughter chromosomes from among the (2Np) better parents and the initial daughters is called elitism since some of the better (elite) parents pass on to the next generation. Though no biological inspirations have been offered,4,7 elitism leads to enhanced speeds of convergence. The procedure continues over several generations. The best set of solutions (in terms of the fitness of the multiple objective functions) emerges after a sufficiently large number of generations (Darwin’s law of survival of the fittest). Assigning higher fitness (selfishness or altruism) to some elite chromosomes in the population has been used in adaptations like selfish genes, Dawkins memes, (benevolent agent) prisoner dilemma in game theory, etc. by several workers8-11 to improve the performance of the algorithm. However, these studies do not explain one important question, as to why altruism should be favored over the natural selection built-in genetic algorithm. In contrast, the haplo-diploid nature of the chromosomes in honey bee colonies favors altruism over Darwin’s natural selection since nature only recognizes the increase in the inclusiVe fitness. In this work altruism is adapted in the framework of NSGA-II-aJG quite differently by mimicking the altruistic behavior of honey bee colonies. In the first step, the algorithm developed is used to simulate bee colony as it actually is using values of the crossover and mutation probabilities that are encountered in nature (in absence of jumping genes), without worrying about its optimization efficiency. While doing this, the origin of altruism in the framework of NSGA-II-aJG is also explained. Thereafter, the optimization efficiency of the algorithm is improved using a new crossover scheme (with higher probability of mutation and also the jumping gene operator) and the multiqueen adaptation. This sequential presentation shows how biological phenomena can be mimicked in optimization algorithms. The new adaptation, Alt-NSGA-IIaJG, is then used to solve three benchmark test problems and a multiobjective optimization (MOO) problem of an industrial phthalic anhydride reactor. In the third part, the algorithm is modified for simulating carcinogenesis. Before we proceed, three benchmark problems and the procedure to analyze the convergence are described below.

10.1021/ie9004817 CCC: $40.75  2009 American Chemical Society Published on Web 06/11/2009

9672

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

2. Benchmark Problems4,7 The ZDT4 Problem. The Zitzler, Deb, and Thiele problem 4,4,7 referred to as the ZDT4 problem, is used as the first test problem for this new algorithm. This problem is described below: min f1 ≡ x1

f1 g(x)

s.t.: 0 e x1 e 1;

n

∑x

(3c)

i ) 1, 2, ..., n

(3d)

9 n-1

i

i)2

with n ) 30. The Pareto-optimal front corresponds to xj ) 0; where j ) 2, 3,..., 30. This is a good test problem for any multiobjective optimization algorithm since the Pareto front is discontinuous.

[ ( )]

min f2 ≡ g(x) 1 -

g(x) ≡ 1 +

1/2

subject to (s.t.)

3. Statistical Analysis of Convergence 0 e x1 e 1 -5 e xj e 1;

In order to analyze the Pareto set of solutions statistically, as they evolve over generations, the standard deviation (σ2), or the mean sum-of-square errors, with respect to the final converged results,12 is used as a measure of convergence. σ2 is defined as

j ) 2, 3, ..., n

where g(x), the Rastrigin function, is given by n

g(x) ≡ 1 + 10(n - 1) +

∑ [x

2 i

- 10cos(4πxi)]

Np

(1)

i)2

where n is (often) taken as 10. This problem has 99 Pareto fronts4 of which only one is the global optimal. The latter corresponds to the first decision variable, x1, lying between 0 and 1 (and so, 0 e f1 e 1). All the other decision variables, xj; j ) 2, 3,..., 10, corresponding to the global Pareto set have values equal to 0 (and so, 0 e f2 e 1). The binary-coded NSGA-II (without the aJG operation) has been found4,7 to converge to local Pareto sets, rather than to the global optimal set. The importance of this problem lies in the slow rate of convergence to the solution. This also occurs in nature in various biological processes, e.g., the evolution of species, cancer, etc. Thus, the above problem is well suited to study, qualitatively, several phenomena in nature and is also used as a model problem for testing new algorithms. In contrast, convergence is much more rapid for two other problems in this family, namely, ZDT2 and ZDT3, and the effects of complex phenomena like altruism, jumping genes, etc., get diluted. Drawing conclusions is quite difficult from these cases. These two problems are described as follows: ZDT2 Problem. This problem is given below min f1 ) x1

(2a)

min f2 ) g(x)[1 - [f1 /g(x)]2]

(2b)

where g(x) ≡ 1 +

9 n-1

s.t.: 0 e xj e 1;

n

∑x

i

(2c)

i)2

j ) 1, 2,..., n

(2d)

with n ) 30. The Pareto-optimal front corresponds to the first decision variable, x1, lying between 0 and 1. All the other decision variables, xj; j ) 2, 3,..., 30, corresponding to the global Pareto set have values equal to 0. The complexity of the problem lies in the fact that the Pareto front is nonconvex. ZDT3 Problem. This problem is given below min f1 ) x1

(3a)

min f2 ) g(x)[1 - {f1 /g(x)}1\2 - {f1 /g(x)}sin(10πf1)] (3b) where g(x) is given by

σ2 )

∑ i)1

(

(f2,i - f2,opt,i) range of f2,opt,i Np

)

2

(4)

In eq 4, (f1,opt,i, f2,opt,i) is the ith point in the converged Pareto set, the points being arranged in increasing order of f1. f2,i is the interpolated value of f2 corresponding to f1,opt,i, in any (earlier) generation. σ2 gives the mean square deviation of the Pareto set in any generation, with respect to the near-optimal Pareto front. σ2 e 0.1 is assumed to correspond to the converged solution, while higher values of σ2 reflect either initial unconverged solutions or final convergence to local Pareto sets. The total fitness of the population in any generation is reflected by the value of σ2; the lower the value of σ2, the better is the total fitness of the population. It may be mentioned that the computation of σ2 requires obtaining the converged solution and then calculating it for the previous generations. Its use is, therefore, somewhat limited (better parameters characterizing the degree of convergence are available4,7 to decide when convergence is attained). 4. Adaptation of Honey Bee Colony Description of Honey Bee Colonies. Each honeycomb has one queen. The honey bee colonies have a large population of daughter-worker bees and a very small number of son-drones. The queen and all the female worker bees are diploid as they have 2n (n ) 16 for this species) chromosomes. In these colonies, son-drones are haploid, having only n chromosomes, these being randomly selected from among the 2n chromosomes of their mother (queen). They, thus, have 100% chromosomal relatedness with the queen mother. Similarly, the sister worker bees (2n chromosomes) produced by the crossing of the fatherdrone (all n chromosomes) and the queen-mother (random n chromosomes from among her 2n), are related to each other by 0.5 × (0.5 + 1), i.e., 75%.13-15 In contrast, if a daughter workerbee crosses with a son-drone to produce granddaughters, the latter are related to their own mother worker-bees by 50% only [Table 1 gives the aVerage chromosomal (chromosomewise) relationships between different members of the bee family, with the queen mating with either a single or several drones]. Hence, a large population of the daughter-worker bees prefers to rear the queen’s offspring (sibling sisters; altruistic behavior) rather than produce their own offspring (selfish behavior). This maximizes the inclusiVe fitness of the population. Initial Altruistic Algorithm. Altruism (Alt) is incorporated in NSGA-II-aJG by taking a single queen (Nqueen ) 1). This can be done in two ways. In the first of these adaptations,

Ind. Eng. Chem. Res., Vol. 48, No. 21, 2009 Table 1. Relationships (Stringwise) between Different Members of the Honey Bee Familya if queen mother, M1, has mated with with several drones, Fj no.

relationship between

a M1 and any Di,j* b M1 and any Si c Di,j and full/half sisters, Di′,j* d Di,j and full/half brothers, Si′ e Di,j and their own sons (queen’s grandsons), GSi,j,k f

Di,j and their nephews, GSi*,j′,k

only from same from different one drone, F1 father, Fj* fathers 0.5 (j* ) 1) 0.5 0.75

0.5 (j* ) same) 0.5 0.75

0.5 (any j*) 0.5 0.25

(j* ) j )1) 0.25 (j* ) j )1) 0.5

(j* ) j) 0.25 (j* ) j) 0.5

(j* * j) 0.25 (j* * j) 0.5

(j )1) 0.375 (j′ ) j )1)

0.375 (j′ ) j)

0.125 (j′ * j)

a

Sample calculations: The queen mother, M1, is chromosomically related to her daughters, Di,j, by an aVerage of 50%. This is because, of her 2n strings, only n randomly selected strings are present in any daughter. The full sisters (from the same father) are chromosomically related to each other by an average of 75% since 50% of the strings from the queen are common (from the egg) and 100% of their n strings from the sperm are common. Hence, the average relationship is 0.5 × (0.5 + 1.0) and so on.

referred to as the one-good-queen-Alt-NSGA-II-aJG, a “good” queen-string is inserted in the population right at the beginning. Such a queen is generated by solving the multiobjective optimization (MOO) problem (in this case, only the ZDT44,7 problem) until near-convergence using the original NSGA-IIaJG algorithm (in absence of altruism). Any one string of the converged Pareto front is selected as the good queen. This good queen will, most probably, pass on unchanged to the next generation (elitism), and its chances of remaining as the queen are quite high. In the second approach, the queen is selected from within the population itself, right from the first generation, as the “best/fittest” (least crowded) string. The procedure is repeated in eVery generation. Hence, the queen could change over generations. Since the initial, randomly generated population is far from the converged solutions, the queen picked up by the second procedure will not be too good, particularly in the initial generations, and is referred to as a “bad” queen. This adaptation is, thus, being referred to as the one-bad-queen-AltNSGA-II-aJG. In this, the queen would, most likely, have a rank number,4,7 Irank, of 1 [see the Appendix] and, since it is the best in any generation, is likely to be transferred unchanged to the next generation through elitism, possibly as the queen. Thus, what happens over time in a honeycomb occurs over generations in these two adapted algorithms through elitism. Clearly, the first of these algorithms would be faster, but it requires a-priori input of a good queen (converged solution of the same problem using the original NSGA-II-aJG). This requires the solution of the MOO problem twice, a futile exercise, and suggests that the one-bad-queen-Alt-NSGA-II-aJG is the preferred approach. However, the former adaptation mimics a bee colony more closely. The altruistic adaptation differs from the natural selection process in NSGA-II-aJG, wherein there are no queens and both the strings are selected randomly from the population. Crossover (as well as mutation and the aJG) operations are carried out only between the one (or more) queen-string and the remaining worker-bee strings. In both the adaptations of the one-queen-Alt-NSGA-II-aJG, of the two strings selected for crossover, mate (or string) 1 is (most likely) selected as the queen-string (in box P′′; see Figure A-1 in the Appendix). The other string, mate 2, is selected randomly from the remaining

9673

Np - 1 worker-bee strings. It is likely (since the queen may change over the generations) that the mate 2 string has been produced from a mate 1 queen-string of a preVious generation. The daughter strings in the present generation are, thus, related to the mate 1 queen-string of the previous generation by 50%. Hence, the current mate 2 string is related to the current mate 1 queen-string by an average of 75% (25% comes from the previous to previous generation). The daughter strings in the current generation [strings 2 - Np] are related to each other, on an average, by 62.5% [0.5 × (0.5 + 0.75); it will actually be 0.5 × (0.5 + [∑1Ngen (1/2N)]), but neglecting the terms for Ngen > 2, the average relationship is 62.5%]. If a worker-bee string is crossed with any other worker-bee string, then the granddaughter (of the current generation queen-string) string will be related to their mother worker-bee strings by 50% only (while the sibling-strings are related by 62.5%). Hence, the current daughter-strings prefer to crossover (and die) with the mate 1 queen-string, than with the other worker-bee strings, 2 to Np, so as to maximize the total fitness of the population. This produces more sisters in the algorithm (or brothers, since both are indistinguishable in Alt-NSGA-II-aJG). Mating among the sister-strings, 2 to Np, is not preferred much in any generation, since this is equivalent to daughter-worker bees producing their own offspring. We do not have the possibility of protecting sisters in evolutionary algorithms, and generating sisters (and not offspring) is an algorithmic means of doing this. This is how altruism is implemented in the proposed algorithm. Haplo-diploidity13-15 is mainly responsible for the relatedness in the bee colony, but the same effect is mimicked in the present approach using the concept of elitism, as described above. There is yet another difference between what occurs in the algorithms and in bee colonies. In the latter, the n chromosomes in the sex cells of the queen and the n chromosomes in the sperm produce daughter-worker bees having 2n chromosomes. This “crossing” involves n crossover operations. While mimicking this operation in GA, the n chromosomes of the sex cells of the queen, as well as of the drone, are represented as a single computational binary string, and the birth of a female worker bee is represented in terms of the crossover of any two strings. In GA, crossover between a mate 1 queen-string and a mate 2 worker-bee string in any generation (to increase the fitness) kills both the starting mates and produces two daughter-strings. The best strings from among the parents and the daughters are passed on to the next generation by elitism. Elitism leads to a further improvement of the total fitness of the population, over and above the effect of altruism. Results for the Initial Algorithms. Figure 1a shows the results (circles) of the ZDT4 problem using the one-good-queenAlt-NSGA-II-aJG. The parameters used are those actually present16 in honey bee colonies, namely, Pcross ) 0.9, Pmut ) 1.6 × 10-6, PJG ) 0.0, Nqueen ) 1. The good queen used in generating these results is the following: xqueen ) [1, 0.00001, 0.00001, 0, 0.00001, 0.00001, 0, 0, 0, 0]T. Results (unfilled squares in Figure A-1) are also generated (with the one good solution inserted externally in the population, but not as a queen) using NSGA-II-aJG (no altruism). No preference is given to the externally introduced good string, that is, natural selection is used. The natural crossover operation is used in both implementations described above (see Figure A-2). It is observed that the one-good-queenAlt-NSGA-II-aJG performs better in terms of the total fitness of the population (lower σ2, reduces faster) as compared to natural selection. Figure 1a also shows the results (filled circles) of this problem obtained using the one-bad-queen-Alt-NSGA-II-aJG. It is observed that use of this algorithm is not too effective when

9674


Figure 1. (a) Standard deviation for the ZDT4 problem using the one-good-queen-Alt-NSGA-II-aJG (circles) and the one-bad-queen-Alt-NSGA-II-aJG (filled circles). (old) Crossover is used. Results (unfilled squares, filled squares) are also shown for the one-good-queen-NSGA-II-aJG and the one-badqueen-NSGA-II-aJG (natural selection), respectively. (b) Standard deviation for four different optimization algorithms for the ZDT4 problem using multiple queens. Technique number given in the text. Several intermediate points have been omitted to ensure clarity.

compared to the results from the one-good-queen-Alt-NSGA-IIaJG. In addition, this algorithm gives worse results than obtained from the one-bad-queen-NSGA-II-aJG (filled squares). This suggests that altruism helps only to propagate good solutions faster, once they have been generated or introduced. In contrast, natural selection is better when a good solution is not introduced externally. Since the one-bad-queen-Alt-NSGA-II-aJG is preferred (for reasons mentioned earlier), we need to develop further adaptations to generate better queens right at the beginning in this algorithm. The objective in this section has been primarily to mimic the altruistic behavior of honey bees in the framework of NSGA-II-aJG. 5. Improving Optimization Improved Altruistic Algorithm, Alt-NSGA-II-aJG. Now that the basic framework of altruistic NSGA-II-aJG has been developed, we now attempt to improve it (somewhat intuitively) so that it can be a more useful algorithm for solving complex, real-life MOO problems in chemical engineering. Two successful adaptations of the one-bad-queen-Alt-NSGA-II-aJG algorithm have been found. Both are incorporated in the final algorithm, named Alt-NSGA-II-aJG. The first adaptation for improving the one-bad-queen-AltNSGA-II-aJG algorithm is to use several, Nqueen (taken as about 10), internally generated parallel queens. In fact, some members of the class hymenopterans, e.g., ants, wasps, etc., have more than one queen. This is akin to studying Nqueen parallel queens in Nqueen parallel honeycombs simultaneously. This leads to more diversity and helps in producing better queens. This adaptation brings in the benefit of altruism (increasing the speed of propagation of better solutions over generations). Some amount of natural selection is still maintained. Exact and complete details are provided in the Appendix. The second adaptation is to use an additional, two-point, three-mate crossover operator (see Figure A-3b). This improves the chances of evolving good internal queens further. In the three-mate crossover procedure, Mate 1 is selected randomly from the Nqueen queen-strings while mates 2 and 3 are selected from the remaining Np - Nqueen worker-bee strings. A hybrid mate is first produced from mates 2 and 3, with its first half as the binaries from 1 to lchr/2 of mate 2, and its second half as

binaries from lchr/2 to lchr of mate 3. A two-point crossover is now carried out using mate 1 and the hybrid mate. Two locations are selected randomly in mate 1, one between 1 and lchr/2 and the other between lchr/2 and lchr. This gives four substrings of lengths (number of binaries), a, b () lchr/2 - a), c, and d () lchr/2 - c), in mate 1. The substrings are now exchanged in one of the following two ways: (a) The lengths and locations of the substrings are the same in the hybrid mate. Substrings a and c remain in mate 1. Substring b is exchanged with the corresponding second substring in the hybrid mate, while substring d is exchanged with the fourth substring in the hybrid mate, (b) The locations of the substrings are altered in the hybrid mate. The first substring in the hybrid mate has length c (and so the length of the second substring is d). Similarly, the third substring in the hybrid mate has length a (and so the length of the fourth substring is b). In crossover, substrings a and c remain in mate 1. Substring b is exchanged with the fourth substring (having the same length) in the hybrid mate, while substring d is exchanged with the second substring of equal length in the hybrid mate (see Figure A-3a). Both of these exchanges are assigned equal probability. The exact and fine details of the complete algorithm incorporating the multiqueen adaptation and the new crossover scheme, referred to as Alt-NSGA-II-aJG, are provided in the Appendix. These adaptations, in some manner, represent the improvement of the broods associated with the queen mating with a diverse set of drones.17 This new crossover operation increases the chance of evolving good solutions, though it does not follow actual genetics, unlike the simple crossover used in the onebee adaptations. The values of some of the probabilities of the common GA operators are significantly higher (best values: Pcross ) 0.9, Pmut ) 1/lchrom, PJG ) 0.4, Nqueen ) 10) than those used in the one-queen adaptations. Results Obtained for the Three Benchmark Problems using Alt-NSGA-II-aJG. The results of the ZDT4 problem for four different algorithms (with identical GA parameters, referred to as the reference values) are shown in Figure 1b. It is observed that Alt-NSGA-II-aJG is the best, with the (near-) optimal solutions obtained in about 50 generations. This is because the new crossover strategy increases diversity and helps produce


9675

Figure 2. Pareto optimal fronts for four adaptations of a multiple-queen algorithm after (a) 50, (b) 100, (c) 150, and (d) 450 generations for the ZDT4 problem (reference values: Pcross ) 0.9, Pmut ) 1/lchrom, PJG ) 0.4, Nqueen ) 10). The small windows of points shown by boxes are also enlarged.

better strings while the multiqueen altruism helps to propagate these better solutions over generations. In fact, use of the new crossover strategy leads to better results even for NSGA-IIaJG (no altruism), than those obtained with the old crossover strategy. Alt-NSGA-II-aJG with the old crossover strategy gives the worst results (worse than even the conventional NSGA-IIaJG) since it does not facilitate the generation of good solutions. The newly developed algorithm, therefore, is quite useful for solving real-life problems in science and engineering, which are computationally intensive. The Pareto optimal fronts of the ZDT4 problem using four different algorithms (based on the multiqueen adaptation), namely, (technique 1) Alt-NSGA-II-aJG (new crossover), (technique 2) Alt-NSGA-II-aJG (old crossover), (technique 3) NSGA-II-aJG (new crossover), and (technique 4) NSGA-IIaJG (old crossover), after (a) 50, (b) 100, (c) 150, and (d) 450 generations are shown in Figure 2. It is observed that nearglobal-optimal Pareto solutions are obtained for Alt-NSGA-IIaJG (new crossover) in about 50 generations, while more than 100 generations are needed for the other three algorithms. The box plots4 indicating the spread of the solutions in terms of four quartiles of the entire set of solutions are plotted for the above algorithms at the 450th generation in Figure 3. Technique 2 is observed to be good for the objective function, f1, while technique 4 is good for f2. Technique 1 covers more points in the 95% ranges (in general, evolutionary algorithms miss the end parts of the Pareto optimal solutions). Technique 1 is superior in terms of convergence and also gives

reasonably well spread-out Pareto optimal solutions as compared to techniques 2-4. It is interesting to observe that the ZDT4 problem fails to converge if the JG operator is not used (i.e., for PJG ) 0), even for Alt-NSGA-II-aJG. Figure 4a shows that σ2 is well above the required criterion of optimal convergence (σ2 e 0.1) for this case, and this is because of convergence to a local Pareto front. Figure 4c shows similar results where the asymptotic value of σ2 for NSGA-II-aJG (new crossover) approaches unity. Figure 4d shows higher values of σ2 for NSGA-II-aJG (old crossover) with PJG ) 0 than for the new crossover scheme. For PJG ) 0, optimal solutions do not emerge even when altruism is used. The importance of the JG operation is well-illustrated by this example. The speed of convergence of Alt-NSGA-II-aJG (new crossover) can be controlled by PJG and Nqueen, the number of queens (see Figure 4a and b). It is observed that varying the values of these two parameters around their reference values leads to a worsening of the rate of convergence. A comparison of Figure 4a and c shows that use of altruism speeds up or maintains the rate of convergence for all values of PJG. In contrast, the speeds of convergence for the ZDT2 and ZDT3 problems4 are quite high (near-optimal solutions are obtained in about 30 generations; see Figure 5), and thus, the benefit due to altruism is smaller. Alt-NSGA-II-aJG (new crossover) gives better results for almost all cases. However, the gains are larger if the speed of convergence is low. Thus, the proposed algorithm is very

9676


Figure 3. Box plots for the four techniques after 450 generations for the ZDT4 problem.

Figure 4. Standard deviation with respect to the near-optimal solution for the ZDT4 problem with (a) Alt-NSGA-II-aJG (new crossover) for varying values of PJG with Nqueen ) 10, (b) Alt-NSGA-II-aJG (new crossover) for varying Nqueen with PJG ) 0.4, (c) NSGA-II-aJG (new crossover) for varying PJG, and (d) NSGA-II-aJG (old crossover) for varying PJG. Reference values of Pcross ) 0.9, Pmut ) 1/lchrom, PJG ) 0.4, and Nqueen ) 10 were used where not varied.

important for solving real-life industrial problems where the speed of convergence is usually low. According to “no free lunch” theorems,18,19 if one algorithm is performing better with one class of problems then it will compensate with another (here every repeated solution should be counted only once). In the case of evolutionary algorithms, the same solution is visited several times (lack of memory) due to their stochastic

nature. Thus, there is a saving in the computational time because of the reduction of the frequency of visiting the same solution again. The altruistic algorithm developed here precisely follows the same strategy since altruism gives more preference for good solutions to crossover and thus speeds up the transfer of good characters to the next generation, which is kind of adding some intelligence to the algorithm.


9677

Figure 5. Standard deviation with respect to near-optimal solutions for the ZDT2 (a and b) and the ZDT3 (c and d) problems and their corresponding Pareto sets at the 100th generation.

MOO of an Industrial Phthalic Anhydride (PA) Reactor System. Alt-NSGA-II-aJG is now used to solve the twoobjective optimization of an industrial PA reactor system.12,20 PA is used widely in industry for the production of vinyl chloride-based plastics and for the production of polyester resins, plasticizers, dyes, insect repellants, etc. It is prepared commercially using gas-phase catalytic oxidation of o-xylene (see Figure 621). Multitubular reactors with several zones of catalyst are used with the intermediate regions between the catalyst packings being used for cooling the flowing gas and removing the exothermic heat of reaction. Such a configuration helps to avoid the temperature of the gas from becoming too high. Two “switch” condensers, used alternately, are used to treat the gas leaving the reactor. The exit gas from the reactor is sent to one of the cooled finned condensers where PA solidifies on the metal surface. The treated gas from the condenser is then scrubbed with water, or incinerated catalytically or thermally. Once a sufficient amount of PA has solidified on the condenser surface, the exit gas from the reactor is sent to the second switch

condenser. During this time, PA is removed from the previous condenser by melting. The reaction scheme suggested by Skrzypek et al.21 (see Figure 6a) is used for modeling the reactor. The LangmuirHinshelwood-Hougen-Watson (LHHW) rate expressions corresponding to this reaction scheme, and the model equations, kinetic, and adsorption parameters are given in reference 20 and are not provided here. Steady state mass and energy balance equations20 between the bulk gas and the outer surface of the imperVious catalyst are solved using the DNEQNF program in the IMSL library. The differential equations corresponding to the steady state balance20 for the bulk gas phase are solved using the code, DIVPAG, in the IMSL library. The simulation results are then used to maximize the yield, XPA, of PA and to minimize the total length of the catalyst bed. This problem involves nine catalyst zones [20 decision variables; see Figure 6b] and is described20 below: max f1(u) ) XPA

(5a)

9678

Ind. Eng. Chem. Res., Vol. 48, No. 21, 2009 9

min f2(u) ) Lcat )

∑

Li

(5b)

i)1

[

subject to Constraints

to a pure maximization problem by using4,7 max f ) 1/(1 + min f)4,7 for the second objective function in eq 5b). Equation 5c is rewritten as follows: max f1(u) ) XPA + w1 1 -

Tmax e 510 °C

(5c)

total length of the reactor tube, Ltube ) 3.5 m

(5d)

Bounds 3

65 e cin e 85 g OX/(m air at NTP)

(5e)

147 e TF,in e 287 °C

(5f)

337 e Tc,in e 447 °C

(5g)

0.001 e m ˙ e 0.005 kg coolant/s

(5h)

0.2 e Si e 0.45 m;

i ) 1, 2, 7

(5i)

max f2(u) )

(

XPA 1.2

]

2

) [

XPA 1 + w1 1 1 + Lcat 1.2

[

(5j)

0.05 e Li e 0.9 m

(5k)

0.01 e L1 e 0.2 m;

i ) 2, 3, ...,8

]

2

Lcat 3.6

]

2

+ w3 (6a)

+

w2 1 -

Lcat 3.6

]

2

+ w3

(6b)

(values of XPA ) 1.2 and Lcat ) 3.6 m in the penalty functions are somewhat arbitrary and are used to ensure that these two variables do not exceed expected or real values) subject to 8

L9 ) (3.25 -

∑

8

Li -

i)1

0.1 e S8 e 0.45 m

[

+ w2 1 -

∑S) g 0 i

(6c)

i

eqs 5e-51

(6d)

The weighting functions in the penalty terms are given by

(5l)

If XPA e 1.1, w1 ) -500;

else w1 ) 0

(7a)

The upper bound on the feed temperature, TF,in, is taken to be 287 °C since several optimal solutions for the modified problem occur in the range 247 e TF,in e 287 °C. Constraints have been put on the spacings, Si (i ) 1, 2,..., 8) between the nine catalyst zones [see Figure 6b] as well as on the lengths, Li (i ) 1, 2,..., 8 Li - ∑8i 8) of the catalyst zones [note that L9 )(3.25 - ∑i)1 Si)]. In order to ensure near-complete conversion of o-xylene and to take care of the constraints, the objective functions are reformulated using penalty functions (the problem is converted

If L9 e 0 m, w2 ) -3000;

else w2 ) 0

(7b)

If Tmax e 510 °C in bed i, i ) 1, 2, ...,9, w3 ) 0; else w3 ) -3000 + 250(i - 1) If Li g 0.01 m, i ) 1, 2,...,9, w3 ) 0;

else w3 ) -300 (7d)

The decision variables for this problem are ˙ , S1, S2, ...,S8, L1, L2, ...,L8]T u ) [cin, TF,in, Tc,in, m

Figure 6. (a) Reaction scheme21 [reaction numbers indicated on the arrows] (b) and reactor set up with nine beds for PA production.

(7c)

(8)

In eq 8, cin is the mass of OX per cubic meter air at NTP. The values of L1-L8 and S1-S8 are selected optimally by the optimization algorithm. L9 is then computed using Ltube () 3.5 m). The values of wi have been decided using trial and error. They should not be too low (so as to be ineffective in forcing the solutions toward satisfying the constraints), nor too large (so as to distort the contours of the modified objective functions too much and make it difficult to obtain the optimal solutions). The catalyst of Skrzypek et al.21 is used with the diameter of the catalyst particle as 3 mm. Other details of this reactor are the following: diameter of each reactor tube ) 25 mm and the mass flux, G ) 19 455 kg/(m2 h). Additional details are given in the work of Bhat and Gupta.20 The results for this problem are shown in Figure 7a and b. Figure 7a shows how σ2 varies with the number of generations. It is observed that near-converged solutions (σ2 e 0.1) are obtained in about 40 generations, but a better spread of the Pareto points is obtained in about 50 generations using Alt-NSGA-II-aJG. In contrast, NSGA-II does not converge to the final Pareto front. Figure 7b compares the Pareto optimal fronts obtained at the 50th generation, using both Alt-NSGA-II-aJG (converged) and NSGAII-aJG (near-converged). The point A (f1, f2 ) 1.17162, 0.795) is the last point on the Pareto optimal front obtained using Alt-NSGAII-aJG. This is superior to that corresponding to point B (f1 ) 1.17161, f2 ) 1.02) using NSGA-II-aJG. The advantage of AltNSGA-II-aJG is clearly brought out. Here, the CPU time (Pentium


9679

Figure 7. (a) Standard deviation for the industrial PA reactor optimization problem using Alt-NSGA-II-aJG and NSGA-II-aJG.20 (b) Pareto optimal front using both algorithms at 50 generations each (points A and B are identified). Best (reference) values: Pcross ) 0.95, Pmut ) 0.025, PJG ) 0.5, Nqueen ) 10.

4, 3.4 GHz, 1.0 GB RAM) for this problem is 42.0 h for 50 generations using both the algorithms. 6. Simulation of Cancer The algorithm developed here is now extended to solve an interesting problem, namely, the simulation of carcinogenesis. Groups of living forms, e.g., bees, ants, chimpanzees, lions, monkeys, etc., and even cells all exhibit altruistic behavior, where each individual fights for increasing its own fitness.15 They also exhibit the other essential ingredients of evolutionary theory like reproduction, mutation, and selection. The process of carcinogenesis also includes these ingredients.22,23 Hence, Alt-NSGA-II-aJG (or any other similar evolutionary algorithm), can be used with some modifications to simulate carcinogenesis even if these algorithms do not represent its exact nature. Komarova23 has formulated the initiation and growth of cancer in terms of an optimization problem and finds that there are some (optimal) conditions that help initiate it. It is known that, in any organism, stem cells form first. They divide asymmetrically to produce one stem cell and one progeny cell, each. The progeny cells further divide to produce a few normal cells at a time. The normal cells, in turn, undergo cell division mitotically to produce two daughter cells in each step. The stem/progeny cells have a long life and are primarily responsible for growth, behaving altruistically like queen bees in a honeycomb or in our present algorithm. Normal cells, on the other hand, are like worker bees since they produce only two cells per cell division, unlike stem cells which grow much more rapidly via the progeny cells. In addition, normal cells undergo less cell division cycles

as compared to stem cells due to their relatively shorter lives (normal cells do, indeed, produce daughter cells, unlike worker bees which do not produce offspring). The algorithm has some similarity to the somatic evolution (see Figure 8) in Moran’s process24 in which a daughter cell that is produced replaces another one in the population and so the cell population remains constant. As per the mutator phenotype hypothesis of Loeb et al.,25 normal rates of mutation are not sufficient to cause cancer during the lifetime of human beings since somatic evolution of cells is an extremely slow process (the mutation rate is very low and somatic evolution involves mitotic cell division or replication, and no crossover; see Figure 8). However, cancer can get initiated and grow because of multistage initiation and progression processes. Instabilities like mutation and jumping genes are examples of such multistage initiation. In general, the normal rate of mutation is of the order of 10-7 per gene per cell division.25 Much higher values of the mutation probability are required for solving real-life optimization problems in science and engineering. Because of the similarities between Alt-NSGA-II-aJG and the initiation and growth of cancer, as discussed above, the latter can be studied using, for example, the ZDT4 problem. The similarity between the algorithm and cancer extends still further since the algorithm does not converge for this problem when low (normal) values of the probabilities/rates of the jumping gene operation are used. This is similar to the description of the mutator phenotype hypothesis of Loeb et al. Use of higher values of the probabilities of these operations in the algorithm leads to the introduction of binary

9680


Figure 8. Existing code of Alt-NSGA-II-aJG represents the complete evolution process (shown by the vertical arrow). The cancer problem is represented by the large circle. The present code is modified accordingly.

sequences in the strings that are responsible for initiating cancer. This is similar to the initiation of cancer in the presence of multistage initiation and progression processes when inactivation of the tumor suppressor gene, namely, the adenomatous polyposis coli (APC) gene, takes place. This inhibits normal cell regulation mechanism. The cells formed under these conditions can grow faster than healthy cells and have a higher fitness in evolutionary terms. Alt-NSGA-II-aJG needs to be modified for simulating carcinogenesis. Replication (mitotic cell division of normal cells) replaces crossover for this purpose. A parent string, primarily from the queens (stem cells), and another string, primarily from the (Nqueen + 1) to Np worker bees (normal cells), are just copied/ replicated to give two daughter strings. Here, the actual cell dynamics (e.g., producing a stem cell and a progeny cell from one stem cell, etc.) cannot be mimicked exactly, but the daughters produced can be taken as the cells which would have been produced in actual cell dynamics. The stem cells are preserved through elitism, a requirement of cell dynamics, since their fitness values are superior. The population of cancer cells increases rapidly due to their higher fitness once a cancer cell emerges (loss of apoptosis). However, this is curbed by the limited blood/food supply (angiogenesis). In Alt-NSGA-II-aJG, the use of the concept of crowding distance avoids accumulation of too many close-by Pareto points (fitter cancer cells). The use of crowding distance to spread out the Pareto solutions helps explain metastasis, the spread of cancer cells over the body. Cell kinetics is not applied too well in this algorithm. However, the evolutionary framework of genetic algorithm used in this study appears to be closer to reality, unlike the other algorithms used earlier, e.g., game theory, stochastic processes, etc.,26-29 to describe cancer. In the present simulation, the length of each string in the population is taken as 300 binary numbers. The ZDT4 problem is first solved for 100 generations until (near) convergence. The original Alt-NSGA-II-aJG algorithm (with multiple queens and simple crossovers as used in the one-queen adaptations) incorporating crossover (with Pcross ) 0.9), replication (Prepli ) 0.1), mutation, and aJG is used for this. This is akin to the evolution of the system from single cell animals like protozoa to the level of mammals. In general, the first 100 generations (taken arbitrarily) is equivalent to reaching the point of birth of a mammal offspring (the end of meiosis) from single cell protozoa. Thereafter, the algorithm is run up to 10 000 generations in the complete absence of crossover (with Pcross ) 0, Prepli

Figure 9. Typical string used in Alt-NSGA-II-aJG. The cancer producing sequence is shown by the bold 00001 in the string. Table 2. Different Cases Studied for Cancer Simulationa case 1 modified objective functions, objcancer assumed cancer sequence Pcross (0-100 generations) Pcross (100-10 000 generations) Prepli (0-100 generations) Prepli (100-10 000 generations) Pmute PJG a

case 2

case 3

case 4

0.95objorig 0.95objorig 0.95objorig 1.00objorig 00001 0.9 0.0 0.1 1.0 10-7 10-2

00001 0.9 0.0 0.1 1.0 10-7 10-4

00000 0.9 0.0 0.1 1.0 10-7 10-2

00001 0.9 0.0 0.1 1.0 10-7 10-2

Pmute and PJG are the same throughout.

) 1, and with mutation and aJG). The latter step is akin to following the entire life of the mammal child in which there is only mitosis, mutation, and JG (large circle in Figure 8). Any one string in the initial (first 100 generations) nearconverged Pareto set is selected (see the typical chromosome in Figure 9). The sequence of binary numbers, 00001, is found to occur at the end of this chosen string. This end-sequence of binaries is assumed to represent the inactivated APC gene, and so is the binary sequence responsible for cancer. Strings having this sequence at the end will have a greater chance of survival because of elitism and will get selected as a queen. In fact, any string that has the same sequence as that present in a string on the near-converged Pareto front (not necessarily a cancer sequence) will get selected as a queen and can be considered as a stem cell. In addition, the computed values, objorig, of both the objective functions, f1 and f2, of such cancerous cells are artificially decreased (since the code maximizes f) in the modified (for cancer simulation) Alt-NSGA-II-aJG algorithm so that the values, objcancer, for the cancerous cell are 0.92objorig. This is akin to what is done in the selfish genes algorithm. Thus, such cells have a higher value of the fitness function in the code.


9681

Figure 10. (a) Behavior of a population of cancer cells for the four cases of Table 2. The ordinate is computed by taking the number of strings having 00001 at the end, per 100 strings. (b) Corresponding optimization results for these four cases.

Figure 11. Pareto fronts at the (a) 4 000th (18 circles) and the (b) 10 000 (26 circles) generations (case 1). The straight lines indicate the number of replicates for the Pareto solutions shown by circles for the two cases, respectively.

This further ensures that cancerous cells are selected in the next generation because of elitism and so are susceptible to faster growth. Four cases are studied for carcinogenesis. These are given in Table 2. In cases 1, 2, and 4, the sequence, 00001, at the end of a string, is assumed as the gene responsible for cancer. Case 3, however, works with a different sequence, 00000, at the end of a string as responsible for cancer. In the nearconverged (100th generation) Pareto set, no string has this end-sequence and so the chances of propagation of such strings by elitism are low. The objective functions of strings with the end sequence, 00000, are decreased for this case, too. However, these strings still have to compete with strings having the end sequence, 00001 (which have higher fitness values inherently since this end sequence is present in the converged optimal Pareto set) for the process of elitism. The probabilities, Pjump, are varied in the four cases studied to see their effect. Figure 10a shows that if the genetic sequence (00001 at the end) responsible for cancer is found in the stem cells and if

reasonably large values are taken for Pjump (sufficient amount of instabilities are present; case 1), then cancer progresses. In contrast, if the cancer sequence (00000 at the end; case 3) with reduced objective functions is found in normal cells or the rate of the instabilities is not sufficiently large (case 2), then it eventually dies out over time/generations. In cases 1-3, we have assumed a higher fitness value of the cancerous cells. In case 4, the fitness of cancerous cells is not increased and it is observed that cancer does not propagate under these conditions. It is also observed that the cancer progresses even if PJG (similar to rate of string loss) is ∼10-2. This is in agreement with the experimental findings of Lengauer et al.30,31 It is also observed (results not shown) that the rate of convergence improves as PJG increases above 10-4. Faster convergence for cancer simulation appears to parallel similar observations for the optimization of the test problem, ZDT4, as shown in Figure 10b. This shows the similarity of optimization and the progress of cancer using GA. Figure 11a and b shows the Pareto sets at the 4 000th and the 10 000th generations for case 1 (of the carcinogenesis

9682


problem), respectively. The numbers of replicates (same values of f1 and f2) of each solution are also shown. It is observed that there is some clustering of the cancer cells at the 4 000th generation. This is because as soon as a cancer cell evolves, several copies of it get formed (similar to the loss of apoptosis), both by elitism and replication. As the number of generations increases to 10 000, use of the crowding distance (angiogenesis) spreads out the solutions to 26 points, representing 26 cancer sites. This spreading of the solutions over the entire domain is similar to metastasis. The close resemblance of the algorithm to the process of cancer evolution makes GA an attractive choice for simulating cancer, as compared to other algorithms. 7. Conclusions The present study highlights the use of altruism as present in bee colonies in optimization. The developed algorithm is used to explain the social conflicts present in bee colonies, used for solving optimization problems in science and engineering and also used for describing complex phenomena like carcinogenesis. In bee evolution, it explains the condition in which altruistic behavior is favored over natural selection. The presence of multiple queens in the other species of hymenopterans is also mimicked to improve the two initial one-queen algorithms. The final algorithm, Alt-NSGA-II-aJG, increases the speed of convergence of three benchmark optimization problems studied in this work, as compared to other algorithms based on GA. It also performs better than conventional NSGA-II-aJG for a reallife industrial MOO problem involving PA tubular reactors. The origin and behavior of cancer are explained using evolutionary principles like reproduction, mutation, jumping gene, and altruism.

Figure A-1. Flowchart for Alt-NSGA-II-aJG.

Acknowledgment We are delighted to dedicate this paper to Dr. B. D. Kulkarni. Partial financial support from the Department of Science and Technology, Government of India, New Delhi [through grant SR/S3/CE/46/2005-SERC-Engg, dated November 29, 2005], is gratefully acknowledged. Appendix: Alt-NSGA-II-aJG, the Binary-Coded Altruistic Elitist Nondominated Sorting Genetic Algorithm with the aJG Operation [Note: The following assumes that we are minimizing all of the m objective functions, fi (i ) 1, 2,..., m)] 1. Generate box, P (see Figure A-1), of Np Parent Strings Each string comprises of a set of lchr binary numbers, with lchr/Nd (integral) binary numbers representing each of the Nd decision variables. A string is made by generating several random-numbers, RNs, sequentially, using a random-number code. These strings are given a sequence (position) number, as generated. 2. Nondominated Fronts Classify these strings into fronts based on nondomination,4 as follows: (a) Create a new (empty) box, P′, of size Np (each of the Np locations can accommodate only one string). (b)Transfer (and remove) the ith string from P to P′, starting with i ) 1.

Figure A-2. Schematic of the simple two-mate crossover strategy (occurring in actual genetics) used in the one-queen adaptation.

(c) Compare string i with eVery member, say, j, already present in P′, one at a time. (d) If i dominates4 over j (i.e., i is superior to or better than j in terms of all the objective functions), remove the jth string from P′ and put it back in its original location in P. (e) If i is dominated by j, remove i from P′ and put it back in its original position in P. (f) If i and j are nondominated (i.e., there is at least one objective function of i that is superior to that of j, while the remaining objective functions are inferior), keep both i and j in P′ (in sequence). Test for all j present in P′. (g) Repeat for the next string in P (in sequence, without going back), until all Np have been tested. P′ now contains a subbox (of size e Np) of nondominated strings (a subset of P) referred to as the first front or sub-box. Assign these strings a rank number, Irank, of 1.


9683

(d) Assign large values of Idist to solutions at the boundaries (the convergence characteristics could be influenced by this choice). 4. Generate a Mating Pool (Box P′′) (a) One-Queen Adaptation (a1) One-Good-Queen-Alt-NSGA-II-aJG. Solve the original problem using NSGA-II-aJG until convergence and select one of the solutions of the final Pareto optimal front. Insert this good solution in the initial population as a good queen. For the following generations, follow steps 4b-e with only one queen. (a2) One-Bad-queen-Alt-NSGA-II-aJG. Follow the procedure given in 4b-e with only one queen from the starting generation. (b) Multiple-Queen Adaptation Generate an integral random number (IRN) between 1 and Np to select a string (randomly) from box P′. Generate a random number (RN) in [0, 1] If RN lies in [0, 0.9]: If the string selected (IRN) is from rank 1, accept it (make a copy without deleting it from P′). If the string is not from rank 1, discard it and repeat with a new IRN/RN. If the RN lies in [0.9, 1]: Accept the selected string irrespective of its Irank. The ith string is, thus, selected. Thus, there is a larger than 90% probability that the selected string has a rank of 1. (c) Repeat this procedure to select a second (jth) string. (d) Select the better of these two strings and put it as string no. 1 in box P′′. Do this a total of Nqueen times and mark and store the strings as nos. 2, 3,..., Nqueen, in P′′. At this point, we have Nqueen selected strings numbered 1 to Nqueen in P′′, of which over 90% have Irank ) 1. (e) String numbers (Nqueen + 1) to Np in P′′ are selected using a procedure that is only slightly different. Identify two strings randomly but from the entire Np in P′ and select the better of the two. Only some of these will be from front 1.

5. Crossover (Replication for Cancer Simulation)

Figure A-3. Schematic of the new (a) two-mate and (b) three-mate crossover strategies. The hybrid mate, HM, is first generated. This then undergoes crossover with mate 1, as shown in part a.

(h) Create subsequent fronts in (lower) sub-boxes of P′, using step 2b above on the strings remaining in P. Compare these members only with the members present in the current subbox, and not with those in the earlier (better) sub-boxes. Assign these values of Irank ) 2, 3,... Finally, we have all Np strings in P′ boxed into one or more fronts. 3. Crowding Distance Evaluate the crowding distance, Idist, for the ith string in any front, j, of P′ using the following procedure: (a) Rearrange all strings in front j in ascending order of the values of any one (say, the rth) of the m objective (fitness) functions. This provides a sequence and, thus, defines the nearest neighbors of any string in front j. (b) Find the largest hypercube (rectangle for two fitness functions) enclosing string i that just touches its nearest neighbors in the f-space. (c) Idist ≡ (sum of m independent sides of this hypercube).

(a) Generate an RN in [0, 1] If RN lies between [0, Pqueen () 0.9)]: Generate an IRN between [1, Nqueen] and select a string (mate 1). Similarly, generate two IRNs between [(Nqueen + 1), Np] and select mates 2 and 3, respectively. Do not remove from P′′ (for the one-queen adaptation, Nqueen ) 1). If RN lies between [Pqueen, 1.0]: Generate three IRNs between [1, Np] and select three strings (name these sequentially as mates 1, 2, and 3). (b) Simple Crossover for the One-Queen Adaptation. Using Pcross and an RN in [0, 1], decide if crossover has to be performed. If yes, then perform a two-mate crossover (see step b1 below) between mates 1 and 2. The two daughter strings produced are put in box D. If crossover is not to be performed, select mates 1 and 2 as the daughter strings. Repeat until box D has Np daughter strings. (b1) Simple Two-Mate Crossover. Generate an IRN between [1, lchr] to give the crossover location. Carry out crossovers at location, IRN, between mates 1 and 2 (mate 3 is not used). See Figure A-2. (c) New Crossover for the Multiple-Queen Adaptation. Using Pcross and an RN in [0, 1], decide if crossover has to be performed. If yes, then generate another RN in [0, 1]. If this is between [0, 0.9], then perform a two-mate crossover (see step c1 below) between mates 1 and 2. If the RN is between [0.9, 1], perform a three-mate crossover (see step c2 below) between mate 1 and mates 2 and 3. Two daughter strings are produced in either case.

9684


These are put in box D. If crossover is not to be performed, select mates 1 and 2 as the daughter strings. Repeat until box D has Np daughter strings. (c1) Two-Mate Crossover. Generate an IRN (IRN1) between [1, (lchr/2); lchr is an even number] to give the first crossover location. Similarly, generate another IRN (IRN2) between [(lchr/2) + 1, lchr] to give the second crossover location. Generate an RN in [0, 1] (see Figure A-3a). If RN is in [0, 0.5]: Carry out crossovers at locations IRN1 and IRN2 between mates 1 and 2 (mate 3 is not used). See Figure A-3a(i). If RN is in [0.5, 1.0]: Select crossover sites as IRN1 and IRN2 in mate 1. For mate 2, select the crossover sites as [(lchr/2) + IRN1)] and [IRN2 - (lchr/2)]. Carry out crossover as shown in Figure A-3a(ii). (c2) Three-Mate Crossover. First generate a single hybrid mate (HM) from mates 2 and 3, having its first half as the first half of mate 2 and its remaining half as the second half of mate 3 (see Figure A-3b). Carry out a two-mate crossover as in step c1 above. (d) Replication for the Cancer Problem. Here Pcross ) 0 and Prepli ) 1. Thus, both mates 1 and 2 are copied unchanged as daughters (for the starting 100 generations, a small amount of replication and simple crossover is used with Prepli ) 0.1 and Pcross ) 0.9).

6. Mutation Select a string (sequentially) from D. Then select a binary number (sequentially) in this string. Check if the mutation operation is to be carried out, using Pmute and an RN in [0, 1]. If yes, alter the selected binary number from 0 to 1 or from 1 to 0 (if no, the binary number remains unchanged). 7. Fixed-Length JG Operation (aJG) Select a string (sequentially) from D. Check if the aJG operation is needed, using Pjump and an RN in [0, 1]. If yes (else go to the next string): (a) Generate an IRN between [1, lchr - fb]. This represents the location of the beginning of the aJG. (b)The end of the aJG is at IRN + fb. (c) Replace the set of binaries between these locations by a randomly generated set of new binaries. 8. Elitism Copy all the Np good parents from box P′′ and all the Np daughters (after crossover, mutation and the JG operations) from box D into a new box, PD. Box PD has 2Np locations. (a) Reclassify these 2Np strings into fronts (in box PD′) using only nondomination (as described in step 2). (b) Take the best Np (lowest front numbers) from box PD′ and put into box P′′′. Calculate and use the crowding distance if a full sub-box needs to be broken up. 9. Stopping This completes one generation. Stop if appropriate criteria are met, e.g., the generation number is larger than a userspecified maximum number, maxgen, of generations. Copy P′′′ into the starting box, P. Go to step 2. Nomenclature cin ) concentration of OX in feed used in phthalic anhydride reactor, g OX/(m3 air at NTP) f1, f2 ) objective functions

G ) mass flux of the gas phase, kg/(m2 s) ∆Hi ) enthalpy of the ith reaction, J/mol L ) length of the reactor tube used in phthalic anhydride reactor, m Li ) length of ith catalyst zone used in phthalic anhydride reactor; i ) 1, 2,..., 9, m Lcat ) total length of all catalyst zones in phthalic anhydride reactor, m Ltube ) total length of the reactor tube in phthalic anhydride reactor, m m ˙ ) coolant flow rate used in phthalic anhydride reactor, kg/s Np ) population size Nqueen ) number of queen chromosomes present in the population Pi ) probability of (i ) cross) crossover, (i ) mut) mutation, (i ) JG) jumping gene, (i ) repli) replication ri ) rate of the ith reaction; i ) 1, 2,..., mol/(m3 s) Si ) spacing after the ith zone of catalyst, m T ) bulk gas temperature at any axial location in phthalic anhydride reactor, K TF,in ) feed temperature in phthalic anhydride reactor, K Tc,in ) coolant temperature in phthalic anhydride reactor, K Xi ) mass-based yield of i in phthalic anhydride reactor, kg species i produced per kg of OX consumed (also, variables used in problem 2) Greek Symbols σ2 ) mean sum-of-square errors, with respect to the final converged results

Literature Cited (1) Kirkpatrick, S.; Gelatt, C. D.; Vecchi, M. P. Optimization with Simulated Annealing. Science 1983, 220, 671–680. (2) Dorigo, M. Optimization, Learning and Natural Algorithms. Ph. D. Thesis, Politecnico di Milano, Italy, 1992. (3) Holland, J. H. Adaptation in Natural and Artificial Systems; University of Michigan Press: Ann Arbor, MI, 1975. (4) Deb, K. Multi-objectiVe Optimization Using EVolutionary Algorithms; Wiley: Chichester, UK, 2001. (5) McClintock, B. The Collected Papers of Barbara McClintock; Garland: New York, 1987. (6) Kasat, R. B.; Gupta, S. K. Multiobjective Optimization of Industrial Fluidized-bed Catalytic Cracking Unit (FCCU) using Genetic Algorithm with the Jumping Genes Operator. Comput. Chem. Eng. 2003, 27, 1785– 1800. (7) Coello Coello, C. A.; Lamont, G. B.; Veldhuizen, D. A. V. EVolutionary Algorithms for SolVing Multi-objectiVe Problems, 2nd ed.; Springer: New York, 2007. (8) Pospı´chal, J.; Kvasnicˇka, V. A Study of Altruism by Genetic Algorithm. In AdVances in Soft Computing - Engineering Design and Manufacturing; Springer: London, 1999. (9) Kvasnicˇka, V.; Pospı´chal, J. Simulation of Baldwin Effect and Dawkins Memes by Genetic Algorithm. In AdVances in Soft Computing Engineering Design and Manufacturing; Springer: London, 1999. (10) Glomba, M.; Filak, T.; Kavasnicka, H. Discovering Effective Strategies for the Iterated Prisoner’s Dilemma using Genetic Algorithm. In Proceeding of the 2005 5th International Conference on Intelligent System Design Application (ISDA′05), IEEE Computer Society: Wroclaw, 2005. (11) Bukhari, S.; Adnan, H. A. S. Using Genetic Algorithm to Develop Strategies for the Prisoner’s Dilemma. Asian J. Infor. Technol. 2006, 5, 866–871. (12) Ramteke, M.; Gupta, S. K. Biomimetic Adaptation of the Evolutionary Algorithm, NSGA-II-aJG, using the Biogenetic Law of Embryology for Intelligent Optimization. Ind. Eng. Chem. Res., published online Feb 17, http://dx.doi.org/10.1021/ie801592c. (13) Hamilton, W. D. The Genetical Evolution of Social Behaviour I. J. Theor. Biol. 1964, 7, 1–17. (14) Hamilton, W. D. The Genetical Evolution of Social Behaviour II. J. Theor. Biol. 1964, 7, 17–52. (15) Gadagkar, R. SurViVal Strategies of Animals: Cooperation and Conflicts; Harvard University Press: Cambridge, MA, 1997. (16) Kerr, W. E. Sex Determination in Bees (Apinae and Meliponinae) and its Consequences. Braz. J. Genet. 1997, 20, 601–611.

Ind. Eng. Chem. Res., Vol. 48, No. 21, 2009 (17) Mattila, H. R.; Seeley, T. D. Genetic Diversity in Honey Bee Colonies Enhances Productivity and Fitness. Science 2007, 317, 362–364. (18) Wolpert, D. H.; Macready, W. G. No Free Lunch Theorems for Optimization. IEEE Trans. EVol. Comput. 1997, 1, 67–82. (19) Koppen, M. On the Benchmarking of Multiobjective Optimization Algorithm. In Knowledge-Based Intelligent Information and Engineering Systems; Springer: Berlin, 2003. (20) Bhat, G. R.; Gupta, S. K. MO Optimization of Phthalic Anhydride Industrial Catalytic Reactors using Guided GA with Adapted Jumping Gene Operator. Chem. Eng. Res. Des. 2008, 86, 959–976. (21) Skrzypek, J.; Grzesik, M.; Galantowicz, M.; Solinski, J. Kinetics of Catalytic Air Oxidation of o-Xylene over a Commercial V2O5-TiO2 Catalyst. Chem. Eng. Sci. 1985, 40, 611–620. (22) Wodarz, D.; Komarova, N. L. Computational Biology of Cancer: Lecture Notes and Mathematical Modeling; World Scientific: Singapore, 2005. (23) Komarova, N. Does Cancer Solve an Optimization Problem. Cell Cycle 2004, 3, 840–844. (24) Moran, P. The Statistical Processes of EVolutionary Theory; Clarendon, Oxford, 1962. (25) Loeb, L. A.; Springgate, C. F.; Battula, N. Errors in DNA Replication as a Basis of Malignant Changes. Cancer Res. 1974, 34, 2311– 2321.

9685

(26) Gatenby, R. A.; Gawlinski, E. T. The Glycolytic Phenotype in Carcinogenesis and Tumor Invasion: Insights through Mathematical Models. Cancer Res. 2003, 63, 3847–3854. (27) Gatenby, R. A.; Vincent, T. L. An Evolutionary Model of Carcinogenesis. Cancer Res. 2003, 63, 6212–6220. (28) Gatenby, R. A.; Vincent, T. L. Application of Quantitative Models from Population Biology and Evolutionary Game Theory to Tumor Therapeutic Strategies. Mol. Cancer Ther. 2003, 2, 919–927. (29) Moolgavkar, S. H.; Knudson, A. G., Jr. Mutation and Cancer: A Model for Human Carcinogenesis. J. Natl. Cancer Inst. 1981, 66, 1037– 1052. (30) Lengauer, C.; Kinzler, K. W.; Vogelstein, B. Genetic Instability in Colorectal Cancer. Nature 1997, 386, 623–627. (31) Lengauer, C.; Kinzler, K. W.; Vogelstein, B. Genetic Instability in Human Cancer. Nature 1998, 396, 643–649.

ReceiVed for reView March 24, 2009 ReVised manuscript receiVed May 19, 2009 Accepted May 20, 2009 IE9004817

Biomimicking Altruistic Behavior of Honey Bees in Multi-objective

Recommend Documents