Parallel Multiobjective Evolutionary Algorithms for ... - ACS Publications

Protection Agency, 26 W. Martin Luther King Drive, Cincinnati, Ohio 45268 ... range of the Pareto front than other multiobjective genetic algorithms. ...
0 downloads 0 Views 206KB Size
Download PDF
Ind. Eng. Chem. Res. 2004, 43, 2669-2679

2669

PROCESS DESIGN AND CONTROL Parallel Multiobjective Evolutionary Algorithms for Waste Solvent Recycling Ki-Joo Kim† and Raymond L. Smith* National Risk Management Research Laboratory, Office of Research and Development, U.S. Environmental Protection Agency, 26 W. Martin Luther King Drive, Cincinnati, Ohio 45268

Waste solvents are of great concern to the chemical process industries and to the public, and many technologies have been suggested and implemented in the chemical process industries to reduce waste and associated environmental impacts. As reported in this article, we have developed a novel parallel multiobjective steady-state genetic algorithm (pMSGA) for designing environmentally benign and economically viable processes for waste solvent recycling. This pMSGA can efficiently solve this complex multiobjective design problem and provide accurate and uniform Pareto-optimal (i.e., tradeoff) solutions. In addition, it can approximate a wider range of the Pareto front than other multiobjective genetic algorithms. As a case study, acetic acid recovery from aqueous waste mixtures is investigated under simultaneous maximization of the total profit and minimization of the potential environmental impacts (PEIs). At low acetic acid feed contents (xF ) 0.25), many of the Pareto-optimal solutions are economically infeasible and also provided minimal PEI reduction. However, at medium and high feed contents (xF ) 0.30 and 0.35), the total profit is very large, and the PEI reduction is significant as well. 1. Introduction Solvents are widely used as dissolving, separating, drying, and cooling agents and as reaction media in the chemical processing industries. Solvents are essential in making these chemical processes economically feasible; however, waste solvents from the chemical process industries are a main source of pollution in the environment if they are not properly controlled. In addition, waste solvents reduce economic performance because loss of solvents means loss of valuable materials. Increasing public awareness of environmental, health, and safety issues drives governments and regulatory agencies to tighten environmental regulations and boost the number of new regulatory laws and policies. To cope with these social and political trends, the chemical process industries have applied many innovative ideas and technologies to minimize waste solvent emissions. One of the commonly practiced methods is the recovery and recycling of useful solvents from waste mixtures with or without entrainers.1-4 However, the act of designing a waste solvent recovery process does not guarantee a benefit to the environment. Even if a proposed process is economically favorable, it might generate more environmental impacts than the original wastes, as shown in Figure 1. For example, low recovery rates and high energy usage might increase the total environmental impacts while allowing the process to remain economically feasible. The opposite scenario, negative profit and low environmental impacts, might * To whom correspondence should be addressed. Tel.: 513569-7161. Fax: 513-569-7111. E-mail: [email protected]. † ORISE Post-Doctoral Fellow. Tel.: 513-569-7644. E-mail: [email protected].

Figure 1. Example plot of profit and environmental impacts.

also be possible. Therefore, the design of environmentally benign processes should incorporate both economic and environmental benefits simultaneously, which requires an efficient multiobjective programming technique. Much of the recent research on waste solvent recovery has considered environmental quality as well as economic performance. For example, Kim and Diwekar2 developed a new multiobjective optimization model for simultaneous integration of solvent selection and solvent recovery and applied it to acetic acid recovery using a continuous separation process. The design problem had four objectives with continuous and integer design parameters, and the multiobjective optimization problem was solved using the -constraint method. Charkraborty and Linninger3 developed a two-stage combinatorial flowsheet synthesis method for plantwide waste management and simultaneously minimized total cost

10.1021/ie0343162 CCC: $27.50 © 2004 American Chemical Society Published on Web 04/21/2004

2670 Ind. Eng. Chem. Res., Vol. 43, No. 11, 2004

and pollution index using goal programming. The second stage for designing waste solvent recovery processes was expressed as a NLP (nonlinear programming) problem. Because of the complexity and multiobjective nature of this design problem, we have applied a genetic algorithm (GA),5 a type of evolutionary algorithm (EA), as an efficient and robust design technique. Unlike deterministic optimization techniques that are repeated to find multiple solutions at each iteration, GAs can handle multiple solutions at each generation because the GA method is population-based. This feature leads to the easy implementation of GAs in multiobjective optimization problems. Another feature of the GA approach is that it requires neither convexity nor continuity of the given design problems, which are critical conditions for deterministic optimization techniques (especially for derivative-based methods). Because of this capability, GAs can be used to solve complex nonlinear, nonconvex design problems, such as those commonly encountered in real engineering and science optimization problems. For more comparisons between stochastic techniques (e.g., GAs and simulated annealing) and deterministic techniques, refer to the recent book by Diwekar.6 Even though the GA approach can be applied to complex multiobjective design problems and modern computers are very fast, one of the key drawbacks of using a GA is the computational effort necessary because of its population-based nature. If evaluation of an objective function is very expensive, using a GA might be questionable. One of the remedies for this drawback is to use a parallel multiobjective evolutionary algorithm (pMOEA). In parallel computing, the computational burden is (evenly) distributed among the participating node processors, and thus, the overall computational time can be significantly reduced. In this paper, we report the development of a novel multiobjective programming technique, called the parallel multiobjective steady-state genetic algorithm (pMSGA), for designing environmentally benign and economically viable processes for waste solvent recycling. An industrial case study investigated in this paper is acetic acid recovery from aqueous waste mixtures generated by the pharmaceutical industries. Section 2 provides basic background on GAs, multiobjective optimization problems, and MOEAs and then explains our novel MOEA implementation. This implementation is expanded to the parallel computing domain in section 3. Sections 4 and 5 provide case studies with different waste feed compositions. Section 6 concludes this paper with a discussion of future research topics.

Table 1. Pseudocode of the GA 1 2 3 4 5 6 7 8

Generate population Evaluate population while (* stopping criteria) Do selection/crossover/mutation Evaluate population Update population end while Stop

sents how well an individual can survive under given evolutionary conditions (including stresses). More-fit individuals transmit their genetic properties to subsequent generations through their offspring. To create new child individuals, selection, crossover, and/or mutation operators are generally applied. Selection chooses a parent from two randomly selected individuals on the basis of their fitness values. After two parent individuals have been selected, some of their genes are exchanged with the given crossover probability (Prcrossover). These new individuals can undergo a mutation operation to increase genetic diversity. The mutation probability (Prmutation) is generally very low, as it is in actual genetic systems. Finally, over many generations, a population gradually evolves into a group of best-fit individuals. This general procedure of the GA technique is summarized in Table 1. 2.2. Multiobjective Evolutionary Algorithm. MOEAs have been applied in various fields of science and engineering. In the chemical and environmental fields, MOEAs are used in molecular design,7 chemical process simulation,8 cyclone separator design,9 groundwater pollution remediation,10 water quality control,11 water distribution networks,12 and more. A good review paper on the applications of MOEAs in chemical engineering has been published by Bhaskar et al.13 These applications can be described as multiobjective optimization problems (MOPs), formulated as shown below

[]

f1(x) f (x) min z ) f (x) ) 2 ∈Z l fk(x) s.t.

(1)

x ) (x1, x2, ..., xm) ∈ X

where z is called the objective vector, Z the objective space, x the decision vector, and X the parameter space. A decision vector a ∈ X is said to dominate a decision vector b ∈ X if and only if

∀i ∈ {1, 2, ..., k}: fi(a) e fi(b) and

2. Multiobjective Evolutionary Algorithm 2.1. Genetic Algorithm. The GA approach is a population-based optimization technique in which a group of individuals, called a population, undergoes some evolutionary operations such as selection, crossover, and mutation in a search to find the best-fit individual and eventually the best-fit population. An individual (also called a chromosome or string) is used to represent the decision variables of the given design problem and to evaluate the objective function of the design problem. Then, a fitness value of an individual is calculated in terms of the objective function. A fitness value can be a raw objective function value, or it can be a linearly or nonlinearly transformed value that repre-

∃j ∈ {1, 2, ..., k}: fj(a) < fj(b)

(2)

If there is no vector in X that dominates the decision vector a, then a is said to be nondominated or Paretooptimal, after the French-Italian economist and sociologist Vilfredo Pareto.14,15 The entire set of nondominated decision vectors is called the Pareto-optimal set P*, and the corresponding solutions in the objective space are called the Pareto front, denoted by PF*. Hence, the goal of MOP techniques is to accurately represent the true Pareto front. Many mathematical programming techniques are available for multiobjective optimization problems. However, these techniques have two real drawbacks in

Ind. Eng. Chem. Res., Vol. 43, No. 11, 2004 2671 Table 2. Pseudocode of NSGA-II 1 2 3 4 5 6 7 8 9 10 11 12 13 a

Generate population Evaluate population Do nondominated sorting Do fitness sharing while (* stopping criteria) {Intentional blank}a Do selection/crossover/mutation Evaluate population Do nondominated sorting Do fitness sharing Update population end while Stop

{Intentional blank} will be replaced with our implementation.

implementation. The first is that they tend to generate members of the Pareto-optimal set one at a time. The other is that most of them are very sensitive to the shape (e.g., convexity and continuity) of the Pareto front. These limitations can be easily overcome if a GA is applied. Because all GAs are population-based, they find several members of the Pareto-optimal set in a single iteration. In addition, GAs are less sensitive to the shape of the Pareto front because the GA technique does not rely on derivatives of the MOP. Finding and maintaining multiple nondominated solutions requires that some special features be introduced because GAs are normally designed to converge to the best-fit single individual in a population. First, fitness assignment and selection must handle multiple solutions to search toward the Pareto front. A Paretobased approach, which uses the idea of Pareto dominance or nondominated sorting, is the most popular one for this task.16 As a second feature, the MOEA should maintain a diverse population and achieve a welldistributed nondominated set. Fitness sharing,17 which is the most frequently used technique, aims at maintaining well-distributed solutions by degrading the fitness values of more crowded individuals. Another feature to be considered is elitism, introduced by De Jong,18 in which the MOEA keeps some of the previously known best solutions during evolution. The usefulness of elitism in achieving better convergence was thoroughly tested by Zitzler et al.,19 and most modern MOEAs have an elitist strategy. Many MOEAs have been developed, and they mainly differ in their implementations of the above three features. Only two recently reported MOEAs are briefly described here. Srinivas and Deb20 implemented nondominated sorting and fitness sharing in their MOEA, called the nondominated sorting genetic algorithm (NSGA), and Deb et al.21 improved this algorithm with a new fitness-sharing scheme and elitism, denoted their technique by NSGA-II. Table 2 shows pseudocode for the NSGA-II algorithm, where nondominated sorting and fitness sharing are the two distinctive procedures, as compared to the single-objective GA shown in Table 1. Zitzler and Thiele22 proposed a strength Pareto evolutionary algorithm (SPEA) that externally archives the nondominated solutions and performs clustering to reduce the size of the external storage. In this algorithm, nondominated solutions are defined not in terms of distance, but rather by Pareto dominance. They also published a modified version of SPEA, called SPEA2,23 in which they applied a more accurate method to account for Pareto dominance and a new density function.

Figure 2. Cuboid length.

2.3. Multiobjective Steady-State Genetic Algorithm. Our implementation, multiobjective steady-state GA (MSGA), is based on the nondominated-sorting (to search toward the Pareto front) and fitness-sharing (to maintain well distributed nondominated solutions) approaches used in NSGA-II. Nondominated sorting determines the Pareto front for each individual. The overall nondominated solutions are classified as being in the first Pareto front (PF1), and then the next nondominated solutions (after removal of the first Pareto front) are classified as being in the second Pareto front (PF 2), and so on to the last Pareto front. Note that the nondominated solutions after PF1 are not real nondominated solutions. After many generations, the first Pareto front, PF1, should accurately represent the true Pareto front PF*. Fitness sharing in NSGA-II is based on the cuboid length between the adjacent individuals in the same Pareto front. Figure 2 shows the concept of cuboid length in two-dimensional space, where the fitness value of an individual i is the sum of the cuboid lengths between the neighboring individuals (i - 1) and (i + 1). Mathematically, the fitness value Fi is defined as k

Fi )

∑ p)1

k

di,p )

i+1 | f i-1 ∑ p - f p | p)1

(3)

where k is the number of objectives and f ip is the function value of objective p at individual i. If an individual i is largely separated from individuals (i - 1) and (i + 1), then the individual i has a large fitness value and thus has a high chance of survival during the selection operation. Because the fitness value of individual i is determined by the size of the hypercube (e.g., an area in the twodimensional case, as in Figure 2) created by the individuals (i - 1) and (i + 1), Fi does not consider the actual distance between i and (i - 1) [or (i + 1)]. Thus, this fitness sharing might assign rather high fitness values to very similar or identical individuals that should have very low fitness values, resulting in a problem in keeping well-distributed solutions. Using the “sharing-again” approach in our MSGA can easily eliminate this drawback. This new approach, described in Table 3, can be inserted as line number 6 in the pseudocode of NSGA-II in Table 2 and is designed for removing the hypercubes of the least-fit individuals in the first Pareto front. This module is active when Nfirst is greater than Npareto (also called the Pareto size). If active, the least-fit individual in PF1 is removed from

2672 Ind. Eng. Chem. Res., Vol. 43, No. 11, 2004 Table 3. Pseudocode of Sharing Again in MSGAa while (Nfirst > Npareto) Remove the least-fit individual Reevaluate fitness values Nfirst ) Nfirst - 1 end while

6a 6b 6c 6d 6e

aN 1 first ) number of solutions in the first Pareto front PF . Npareto ) number of solutions in the final PF1.

the population, and the fitness values of the remaining individuals in PF1 are reevaluated. This process is repeated until Nfirst equals Npareto. By sequentially removing the least-fit individuals from the first Pareto front, PF1 (and thus PF*) can have well-distributed nondominated solutions. As the name implies, another important feature of our MSGA is that it is a steady-state MOEA, whereas many of the published MOEAs are generational MOEAs. In generational MOEAs, the entire child population is created from the current parent population. However, in steady-state MOEAs, one additional genetic parameter is used, namely, Nrep, which is the number of individuals that should be replaced during each generation and is generally 10% of the population size (Npop). Only Nrep least-fit individuals in the parent population are replaced with new child individuals. All of the remaining better-fit individuals in the parent population are automatically transmitted to the child population without any genetic operations. This feature can significantly reduce the number of function evaluations, and it acts as an internal elitist archive for nondominated solutions (whereas many other MOEAs have external elitist archives). Thus, the Pareto size in our MSGA becomes Npop - Nrep. To test the sharing-again module in our MSGA, the following simple but famous two-objective test function developed by Schaffer24 is used

Problem P1

{

f (x) ) x2 min 1 f2(x) ) (2 - x)2 s.t.

(4)

-10 e x e 10

where Pareto-optimal solutions lie in x ∈ [0, 2]. As the tested three algorithms (NSGA-II, SPEA, and MSGA) have different elitist strategies, their actual population sizes are different. For fair comparison, Npareto, which

is the number of nondominated solutions in the final PF1, is fixed to 30. Figure 3 shows one of the best results of the three MOEA algorithms. Because all of the solutions are on the given objective functions, they are all exact Pareto-optimal solutions. Therefore, the differences among the algorithms are solution size and uniformity. For a quick comparison, we have defined two simple performance metrics: solution size (S) and uniformity (U). Solution size is defined as the total cuboid length of the sorted individuals and measures the length of the Pareto front. Solution uniformity is the ratio of the standard deviation of the cuboid lengths to the average cuboid length, and thus, lower U values means better uniformity. From 30 independent runs, the average S values of the NSGA-II, SPEA, and MSGA are 8.00, 7.86, and 8.00, respectively, and the average U values are 0.67, 0.27, and 0.18, respectively. These two performance indices indicate that MSGA provides both a slightly longer Pareto front (than SPEA) and better solution uniformity (than NSGA-II and SPEA). As described, NSGA-II and MSGA are based on the cuboid length for fitness sharing. Cuboid length evaluation requires sorting of the individuals in the same Pareto front, and this works very well for two-objective optimization problems. However, if the number of objectives is more than two, this fitness sharing might not properly represent the neighboring individuals, resulting in a poor distribution of Pareto-optimal solutions. This critical drawback of NSGA-II and MSGA is removed by defining a new fitness-sharing function. This new function is based on the Euclidean distances from an individual i in the same Pareto front and is defined by nd

Fi )

∑ R(q) δ(i,q)

(5)

q)1

where Fi is the fitness value of an individual i. nd is the number of Euclidean distances from individual i, and its maximum value is Nfirst - 1. R(q) is a weighted constant, and δ(i,q) is the qth minimum distance from individual i. Instead of using all Nfirst - 1 distances, the maximum value (nd) of q used in this research is 4, and the corresponding R(q) values are as follows: R(1) ) 102, R(2) ) 100, R(3) ) 10-2, and R(4) ) 10-4. This practical setting of R(q) values is very useful in speeding up the selection operation.

Figure 3. One set of Pareto-optimal solutions for problem P1 (Npareto ) 30, number of generations ) 5000, string length ) 15, Prcrossover ) 1, Prmutation ) 1/15).

Ind. Eng. Chem. Res., Vol. 43, No. 11, 2004 2673

To test the solution distribution of the new MSGA, we designed a three-dimensional test function as shown below

{

f1(x) ) x1 min f2(x) ) x2 f3(x) ) x3

Problem P2

(6)

g1(x) ) x12 + x22 + x32 g 1

s.t.

where 0 e x e 1. The difficulty of this problem arises from its combinatorial (because of the multiple objective functions) and nonconvex (because of the constraint) nature. Deb et al.25 proposed a similar test function, but their test function is numeric, not combinatorial, and they do not include a constraint. Figure 4 shows typical results obtained for problem P2. Both SPEA and MSGA perform well in approximating the true Pareto front compared to NSGA-II. NSGA-II fails to evenly distribute Pareto-optimal solutions over the Pareto front because of the many duplicated solutions. Because uniformity U cannot be used here, another performance metric, spacing,26 is used to evaluate solution uniformity. Spacing is defined as

spacing )

x

Npareto

1

Npareto - 1

(d h - di)2 ∑ i)1

(7)

where k

| f ip - f jp|), ∑ p)1

di ) min( j

j ) 1, 2, ..., Npareto, j * i

d h ) average of all di and it measures the distance variance of neighboring individuals. Spacing will be 0 if all of the solutions are equally distributed. The resulting spacing values are 0.503 for NSGA-II, 0.020 for SPEA, and 0.015 for MSGA, which shows a better solution distribution for MSGA. 3. Parallel Multiobjective Genetic Algorithms Our MOEA implementation, MSGA, is expanded to a parallel computing domain to reduce computational time, which is the main drawback of GAs. Fortunately, any GA can be easily parallelized because it is a population-based technique. The commonly used GA parallelization model is the master-servant model in which the master processor assigns unevaluated individuals to the servant processors and obtains the results from the processors. This parallelization model is especially efficient if function evaluations use most of the computer time. We have integrated our MSGA into PGAPack,27 which is a master-servant genetic algorithm package for single-objective optimization problems. The resulting parallelization code, pMSGA, was then performed using a Beowulf cluster with 16 Linux machines (Red Hat Linux 7.1). Each cluster machine has an Athlon XP 2000 CPU and 256 MB of RAM. The implemented library for the message passing interface (MPI) is MPICH (version 1.2.1).28

Figure 4. Pareto-optimal solutions of the problem P2 (Npareto ) 100, number of generations ) 5000, string length ) 45, Prcrossover ) 1, Prmutation ) 1/45).

The cluster performance is measured in terms of speedup and efficiency. Speedup is defined as the time it takes to complete an algorithm with one processor divided by the time it takes to complete the same algorithm with np processors and is used to measure relative program performance as the number of processors increases. In other words, speedup provides an indication of the effective number of processors utilized.

speedup )

runtime for single processor runtime for np processors

(8)

≈ np for an ideal case Efficiency is defined as the speedup with np processors divided by the number of processors np and is a measure

2674 Ind. Eng. Chem. Res., Vol. 43, No. 11, 2004

environmental impacts (PEIs). The MOP can be formulated as follows

Problem P3

min

{

f1(x,y) ) -P f2(x,y) ) PEI

(10)

subject to

g1(x,y) ) xHOAc g 0.99 g2(x,y) ) B g 5 0.5 e x1 ) teq e 2.5 1 e x2 ) t e 20 Figure 5. Speedup and efficiency of the cluster.

5 e x3 ) RB e 40 of effective hardware use. An efficiency of 100% means that all of the processors are being fully used at all times.

efficiency )

speedup np

(9)

A small subset of the case study, described in the next section, was used to evaluate the cluster performance. The testing conditions were as follows: Npop ) 100, Nrep ) 50, and number of iterations ) 50. The crossover and mutation probabilities are commonly used values. Figure 5 shows speedup and efficiency results for this test problem. The square of the correlation coefficient of the speedup is 0.9972, meaning that the cluster system is linearly scalable. However, the efficiency drops from 100 to 85%. This means that only 85% of the nodes are actually involved in calculations, resulting in lower speedup numbers (e.g., a speedup of 13.5 at 16 nodes). The main reason for this low efficiency is the masterservant parallelization model, in which the ratio of the number of unevaluated child individuals to the number of servant processors is important. If the remainder of this ratio is 0, then the efficiency would be 100% (except for the efficiency loss due to communication). Otherwise, there are some idle servant processors at each generation. Even if the ratio is chosen to have zero remainder (i.e., Nrep is a factor of np processors), the GA does not guarantee child individuals that are different from their parents. Especially when the two parent individuals are very similar, the child individuals might be very similar or even identical to the parents. If they are identical, the child individuals inherit a tag that indicates completed evaluation, and some of the servant processors remain idle at the last distribution of individuals. For these two reasons, an efficiency loss might always occur in the master-servant modeling of MOEAs. 4. Case Study: Waste Solvent Recycling In this industrial case study, the aim is to recover acetic acid from aqueous waste mixtures from pharmaceutical manufacturing. In the pharmaceutical industries, batch distillation is a preferred separation technique, so this technique is used in this paper. Our goal is to design a waste solvent recycling process for maximizing total profit (P) and minimizing potential

20 e x4 ) V e 120 7 e y1 ) N e 38 where x is a decision vector for continuous variables and y is a decision vector for discrete variables. More specifically, teq is the equilibration or startup time (h), t is the batch time (h), RB is the reboil ratio, V is the boil-up rate (kmol/h), and N is the number of trays. The purity constraint on acetic acid, xHOAc mole fraction in the product, is equivalent to the glacial grade on the market and is required by the upstream process units. Then, one can estimate how much fresh acetic acid (purchased from an external source) can be saved and construct a profit function based on the amount of recovered acetic acid and the purchase (or sales) cost at this grade. The second constraint, bottom production rate B in kilomoles, is rather relaxed to generate more feasible solutions during the evolutionary processes. This promotes rapid and accurate convergence to the Pareto front as the problem is very complicated and highly constrained. The total profit function is in the simple form

P)

BCHOAc TAC ts + teq + t HRs

(11)

The first term represents recovery profit, and the variables CHOAc and ts are the cost of the product in dollars per mole$/mol and the setup time of 2 h, respectively. The second term represents capital and operating costs as described below

TAC K1V ) HRs

0.5

N 0.8 + K2V 0.65 + K3V HRs

(12)

where TAC is the total annual cost; HRs is the number of hours per year, at 90% capacity; K1V 0.5N 0.8 is the cost of the column; K2V 0.65 is the cost of the reboiler and condenser; and K3V is the cost of utilities. The cost coefficients, based on 316 stainless steel using Guthrie’s cost data,29 are K1 ) 3200, K2 ) 20 000, and K3 ) 200 (see Appendix A for details). As the calculation is based on a small tray number, low vapor boil-up rate, and large cost correction factor, these cost coefficients are very conservative.

Ind. Eng. Chem. Res., Vol. 43, No. 11, 2004 2675 Table 4. Normalized Impact Scores of Acetic Acid and Coal Energy impact category, j human toxicity potential by ingestion (HTPI) human toxicity potential by inhalation or dermal exposure (HTPE) terrestrial toxicity potential (TTP) aquatic toxicity potential (ATP) ozone depletion potential (ODP) global warming potential (GWP) photooxidation chemical potential (POCP) acidification potential (AP)

ψj (kg-1)

ψEj (MJ-1)

0.1065 7.83 × 10-5 0.0117 1.22 × 10-6 0.1065 7.83 × 10-5 0.0107 2.65 × 10-4 2.03 × 10-9 1.93 × 10-4 7.07 × 10-8 5.98 × 10-3

The second objective function (PEI), representing potential environmental impacts30 is given by Figure 6. Numbers of Pareto and feasible solutions for xF ) 0.30.

PEI ) PEIresidue + PEIenergy )

∑i ∑j Miψij + QR∑j ψEj

(13)

design variable

As seen in eq 13 above, PEI has two contributions: PEI from residues after distillation and PEI from energy consumption during distillation. Mi is the amount of component i in the residue, and ψij is a normalized impact score of component i in impact category j. QR is the reboiler heat duty, and ψEj is a normalized impact score of energy in impact category j. It is assumed that the weighting factors between impact categories are set to 1. The normalized impact scores of acetic acid and coal energy are summarized in Table 4.31 If the waste mixture is 100 kmol with 30 mol % of acetic acid, the initial PEI of the waste mixture (no PEIenergy contribution yet) is 423.7 impact scores. Of course, batch distillation decreases PEIresidue by separating highly pure acetic acid from the waste mixture. However, it also generates PEIenergy to operate this waste solvent recovery process. Thus, the total PEI after acetic acid recovery should be lower than the initial value, 423.7, for the process to be claimed as environmentally benign. To model this batch waste solvent recycling process, MultiBatchDS,32 a professional multibatch distillation simulation package, was used. The dynamic rigorous model in this package is based on the assumptions of constant molar liquid hold-up, negligible vapor hold-up, adiabatic operation, theoretical plates, and finite-difference approximations for the enthalpy changes. The complete column dynamics are summarized in Table 1 of the paper published by Kim and Diwekar.33 For vapor-liquid equilibrium calculations, the UNIFAC equation was used. Even though variable-reflux (or reboil ratio) mode and optimal-reflux (or reboil ratio) mode for better performance are already included in the rigorous dynamic model, these operating policies are not used in this case study for simplicity. For example, a linear reboil ratio policy can be expressed as

RB(t) ) c1 + c2t

Table 5. Gene Lengths, Step Sizes, and Lower and Upper Bounds of Design Variables

(14)

where c1 and c2 are the parameters to be optimized. This reboil ratio policy not only increases the number of design parameters but also increases the population size and the number of generations needed because of very wide ranges in the c variables.

gene length step size lower bound upper bound

teq

t

RB

VB

N

10 0.001 96 0.5 2.5

18 0.000 07 1 20

18 0.000 13 5 40

15 0.003 05 20 120

5 1 7 38

5. Results and Discussion Multiobjective optimization problem P3 was solved using the developed pMSGA. The genetic parameters used are as follows: string length ) 66 (binary bits), Npop ) 200, Nrep ) 100, number of iterations ) 1000, Prcrossover ) 1.0, and Prmutation ) 1/66. In general, Nrep is 10% of the population size. However, if the problem is complicated, Nrep should be increased so that more regions in the objective domain can be explored at each generation. In this case study, Nrep of 100 is used to generate more feasible solutions in early evolutionary stages, and the resulting Pareto size (Npareto ) Npop - Nrep) is 100. The crossover probability, Prcrossover, is 1.0 to force all of the selected parent individuals to undergo the crossover operation. The crossover type is uniform, and its probability is set to 0.5. Because the string length is 66, the mutation probability, Prmutation, is set to 1/66, which is a typical value. The gene length, step size, and lower and upper bounds of each design variable are listed in Table 5. 5.1. Case 1: Base Case. The nominal feed content of acetic acid in the waste mixture is 0.30 mole fraction, and the waste mixture has 423.7 PEI if the waste mixture is 100 kmol. A total of 17 processors (1 master and 16 servant nodes) were used in the application of pMSGA to problem P3, and this complex multiobjective design problem was completed within 5.8 h. All of the results reported for the case studies are taken from several independent runs. The numbers of Pareto and feasible solutions are shown in Figure 6. The number of solutions in the first Pareto front reaches Npareto ()100) after 500 generations and remains almost constant. The number of feasible solutions is around 150 after 100 generations. Stringent product purity makes for a high rate (about 25%) of infeasible solutions and no distinctive regions for feasible and infeasible solutions in the objective space. This high rate of infeasible solutions requires many generations to approximate the true Pareto front.

2676 Ind. Eng. Chem. Res., Vol. 43, No. 11, 2004

Figure 7. Pareto-optimal solutions for xF ) 0.30.

Figure 8. PEIresidue, PEIenergy, and bottom production for xF ) 0.30.

Infeasible solutions arise for two main reasons. One is improper design for waste solvent recycling. For example, a high boil-up rate with a long batch time might dry up the reboiler, resulting in an error during the batch distillation simulation. The number of improper designs at the initial population is approximately 10% of the population and rapidly diminishes as evolution proceeds. The other reason for infeasible solutions is the strict purity constraint. Even if the given design conditions simulate the waste solvent recycling process successfully, the product purity might not be satisfied, marking them as infeasible. Figure 7 shows the final Pareto front, where only 50 Pareto-optimal solutions (after the final solutions are shared again) are shown here for graphical tidiness. The Pareto front is well curved, except for a discontinuous region around -47 $/h on the x axis. Because of the tight purity constraint, new child individuals from the feasible parent individuals around this discontinuous region might be infeasible. The total profit ranges from $42 to $83 per hour (on average $500,000 per year), and the PEI abatement is 21.5-24.5% from the original PEI. Figure 8 shows the trends of PEI and bottom production. Surprisingly, the acetic acid recovery rate is less than 50% for xF ) 0.30, and this low recovery rate is partially caused by the difficulty of separating acetic acid from water due to their close boiling points.

Figure 9. Design variables for the Pareto-optimal solutions for xF ) 0.30.

Figure 10. Recovery profit and costs for xF ) 0.30.

It should be noted that the low recovery rate is not directly caused by the low constraint value of g2 in eq 10. A higher acetic acid recovery rate might be possible, but in that case, the recovery profit (i.e., the first term in eq 11) might be worse because of a longer batch time, and the TAC (i.e., the second term in eq 11) might also be worse because of the need for larger equipment and severe operating conditions. Thus, these solutions can be easily dominated by other solutions that have a lower recovery rate. As the proposed Pareto-optimal solutions are economically and environmentally viable, closer investigation of designs and detailed outputs will help engineers and decision makers choose better designs for further analysis. Figure 9 shows the members of the Paretooptimal set, in which t and RB change little while N and V experience significant variations. Generally, N and V vary in inverse proportion to each other to maintain product purity. High profit is obtained at conditions of low tray number and high boil-up rate. Even though the recovery profit is not a maximum at these design conditions, the TAC is lower (but not the minimum) than those of many other members in the Pareto-optimal set (see Figure 10), resulting in a high profit in total. Low profit at high tray number is mainly attributed to the longer batch time and higher column cost. Under these design conditions, the boil-up rate is rather mild (i.e., low heat exchanger cost) compared to the values for many other

Ind. Eng. Chem. Res., Vol. 43, No. 11, 2004 2677

Figure 11. Pareto-optimal solutions for xF ) 0.25.

members in the Pareto-optimal set. As N increases and V decreases, the column and heat exchanger costs cross each other, as shown in Figure 10. 5.2. Case 2: Low Waste Feed Composition. Because of variability in pharmaceutical batch processes, the compositions of the waste can vary significantly. In this situation, the effect of feed composition changes on economic and environmental performance and feasibility is of great concern and interest. In this case study, the acetic acid feed mole fraction is reduced to 0.25, and all other genetic parameters are kept the same. Figure 11 shows the final Pareto front, in which many of the Pareto-optimal solutions have negative profit (i.e., positive values along the x axis), even though the proposed Pareto-optimal solutions are environmentally favorable. In addition, the highest profit is just $11 per hour (less than $100,000 per year), which is mathematically feasible but might be infeasible in practice. The PEI of the original waste mixture is 353.1, and the PEI abatement is just 4.8-8.8%. From these results, we can conclude that, for this low feed concentration, waste solvent recycling might not be a good alternative to the current unrecycled option. Adding an entrainer for easy separation of acetic acid might be another option, but it requires proper selection/design of an entrainer and several different operating modes for batch distillation. This is a more complicated multiobjective optimization problem. 5.3. Case 3: High Waste Feed Composition. In contrast to the previous case study, the feed mole fraction of acetic acid in this case is increased to 0.35 to investigate the effect of feed composition change. Figure 12 shows the final Pareto front, where the profits and PEI abatements are large compared to those for the base-case feed composition. The profit spans $120-160 per hour (on average, $1,100,000 per year), and the PEI abatement reaches 34.4-37.1%. The Pareto front of Figure 12 has two main parts: solutions above and solutions below -138 $/h. The solutions with profits above 138 $/h (to the left of -138 $/h in Figures 12-15) undergo small changes in RB and V but incremental changes in N. Small tray number is the main factor in reducing column cost and obtaining the highest profit, which is shown in Figure 14. The solutions with profits below 138 $/h (to the right of -138 $/h) have qualitatively different designs. In this case, N is fixed at 38, the upper bound of N, but all other decisions are now changed to explore different regions

Figure 12. Pareto-optimal solutions for xF ) 0.35.

Figure 13. Design variables for the Pareto-optimal solutions for xF ) 0.35.

Figure 14. Recovery profit and costs for xF ) 0.35.

of the objective space. Of course, any time an optimal design lies on a constraint (such as N ) 38), easing the constraint can provide additional optimal solutions. In a comparison of various feed compositions, the large profits for xF ) 0.35 are due to the large difference between the recovery profit and the TAC. For xF ) 0.30, the recovery rate is less than 50%, whereas for xF ) 0.35, the recovery rate is slightly greater than 50% (actually 50-55%), resulting in large recovery profit values (see Figures 14 and 15). One would expect that even higher recovery rates would not be economically feasible be-

2678 Ind. Eng. Chem. Res., Vol. 43, No. 11, 2004

Figure 15. PEIresidue, PEIenergy, and bottom production for xF ) 0.35.

cause of the need for larger equipment and severe operating conditions. In addition, PEIenergy increases might overwhelm any gains in PEIresidue, and thus higher recovery rates might not be environmentally beneficial either. 6. Conclusions A novel parallel multiobjective steady-state genetic algorithm (pMSGA) for designing environmentally benign and economically viable systems for recycling waste mixtures has been developed. pMSGA has the following main advantages: (i) it can handle multiple solutions at each iteration because of its populationbased nature, (ii) it requires neither convexity nor continuity of the design problem, and (iii) it can significantly reduce the computational burden because of its parallelization. This new technique was applied to three different scenarios varying in terms of waste compositions and provided Pareto-optimal solutions within several hours of computation time. The number of trays (N) and the boil-up rate (V) are the most significant factors in determining total profit and environmental potential impacts. At low acetic acid feed composition, the margins for total profit and environmental potential impacts are relatively small because of the difficult separation between acetic acid and water. In a worst case, the total profit can be negative because of additional resource expenditures such as time and energy. At medium and high feed compositions, the Pareto-optimal solutions provide very high profits and low potential environmental impacts. Further analysis including controllability, flexibility, and safety might be required to choose the best member in the Pareto-optimal set. In addition, waste solvent recovery design coupled with various wastewater treatment scenarios will provide a more complete context for environmentally conscious designs. Acknowledgment We thank ORISE (Oak Ridge Institute for Science and Education) for their financial support for this project. We also thank Dr. Paul Harten, Mr. Jerry Waterman, and Mr. Greg Tucker for building and maintaining the Beowulf cluster. Notation B ) bottom production (kmol) CHOAc ) sales cost of acetic acid ($/kmol)

di,p ) cuboid distance of individual i at objective p F ) feed (kmol) f pi ) pth objective function value of individual i Fi ) fitness value of individual i HRs ) operation hours per year (h) Mi ) residue of component i (kmol) N ) number of stages Npareto ) pareto size Npop ) population size Nrep ) replacement size np ) number of processors P ) total profit ($/h) P* ) Pareto-optimal set PEI ) potential environmental impact PF* ) Pareto front Pr ) probability QR ) reboiler heat duty (MJ) RB ) reboil ratio t ) batch time (h) TAC ) total annual cost ($) teq ) equilibration time (h) ts ) setup time (h) V ) boil-up rate (kmol/h) xF ) feed mole fraction Greek Letters δ(q) ) qth minimum Euclidean distance ψij ) normalized impact score of component i in impact category j ψEj ) normalized impact score of energy in impact category j

Appendix A. Cost Coefficients The base case for cost estimation is described as follows: V ) 50 kmol/h, N ) 10, column cost ) $100,000 × 2 (instead of 1.8, which is the cost factor between 316 SS and carbon steel), heat exchanger cost ) $200,000 × 2, and utility cost ) $20,000. Then, the Guthrie’s cost parameters are K1 ) 3200, K2 ) 20 000, and K3 ) 200. Because the base condition has a small number of trays, a low boil-up rate, and large cost correction factor, the resulting cost coefficients are very conservative. Literature Cited (1) Lee, Y. G.; Malone, M. F. Ind. Eng. Chem. Res. 2000, 39, 2035-2044. (2) Kim, K.-J.; Diwekar, U. M. Ind. Eng. Chem. Res. 2002, 41, 4479-4488. (3) Chakraborty, A.; Linninger, A. A. Ind. Eng. Chem. Res. 2002, 41, 4591-4604. (4) Kim, K.-J.; Diwekar, U. M.; Tomazi, K. Chem. Eng. Commun., manuscript accepted. (5) Goldberg, D. Genetic Algorithms in Search, Optimization, and Machine Learning; Addison-Wesley: Reading, MA, 1989. (6) Diwekar, U. M. Introduction to Applied Optimization; Kluwer Academic Publishers: Dordrecht, The Netherlands, 2003. (7) Jones, G.; Brown, R. D.; Clark, D. E.; Willett, P.; Glen, R. C. In Proceedings of the Fifth International Conference on Genetic Algorithms; Morgan Kaufmann Publishers: San Mateo, CA, 1993; pp 597-602. (8) Hinchliffe, M.; Willis, M.; Tham, M. In Proceedings of the Third Annual Conference on Genetic Programming; Morgan Kaufmann Publishers: San Mateo, CA, 1993; pp 134-139. (9) Ravi, G.; Gupta, S. K.; Ray, M. B. Ind. Eng. Chem. Res. 2000, 39, 4272-4286. (10) Erickson, M.; Mayer, A.; Horn, J. In First International Conference on Evolutionary Multi-Criterion Optimization; Lecture Notes in Computer Science No. 1993; Zitzler, E., Deb, K., Thiele, L., Coello Coello, C. A., Corne, D., Eds.; Springer-Verlag: New York, 2001; pp 681-695.

Ind. Eng. Chem. Res., Vol. 43, No. 11, 2004 2679 (11) Reed, P. M.; Minsker, B. S.; Goldberg, D. E. J. Hydroinf. 2001, 3 (2), 71-89. (12) Halhal, D.; Walters, G. A.; Quazar, D.; Savic, D. A. J. Water Resources Planning Manage. ASCE 1997, 123 (3), 137-146. (13) Bhaskar, V.; Gupta, S. K.; Ray, A. K. Rev. Chem. Eng. 2000, 16 (1), 1-54. (14) Pareto, V. Cours d’Economie Politique; Lausanne: Rouge, Switzerland, 1896; Vol. I. (15) Pareto, V. Cours d’Economie Politique; Lausanne: Rouge, Switzerland, 1897; Vol. II. (16) Van Veldhuizen, D. A.; Lamont, G. B. Multiobjective Evolutionary Algorithm Research: A History and Analysis; Technical Report 9803; Air Force Institute of Technology: WrightPatterson AFB, OH, 1998. (17) Goldberg, D. E.; Richardson, J. In Genetic Algorithms and Their Applications (ICGA’87); Grefenstette, J. J., Ed.; Lawrence Erlbaum Associates, Publishers: Mahwah, NJ, 1987; pp 41-49. (18) De Jong, K. An Analysis of the Behavior of a Class of Genetic Adaptive Systems. Ph.D. Thesis, University of Michigan, Ann Arbor, MI, 1975. (19) Zitzler, E.; Thiele, L.; Deb, K. Evol. Comput. 2000, 8 (2), 173-195. (20) Srinivas, N.; Deb, K. Evol. Comput. 1994, 2 (3), 221-248. (21) Deb, K.; Agrawal, S.; Pratap, A.; Meyarivan, T. A Fast Elitist Nondominated Sorting Genetic Algorithm for Multi-objective Optimization: NSGA-II; Technical Report KanGAL-2000/01; Indian Institute of Technology; Kanpur, India, 2000. (22) Zitzler, E.; Thiele, L. INFORMS J. Comput. 1997, 9 (3), 231-250. (23) Zitzler, E.; Laumanns, M.; Thiele, L. SPEA2: Improving the Strength Pareto Evolutionary Algorithm; Technical Report 103; Swiss Federal Institute of Technology: Zurich, Switzerland, 2001.

(24) Schaffer, J. D. Some Experiments in Machine Learning Using Vector Evaluated Genetic Algorithms. Ph.D. Thesis, Vanderbilt University, Nashville, TN, 1984. (25) Deb, K.; Thiele, L.; Laumanns, M.; Zitzler., E. Scalable Test Problems for Evolutionary Multi-objective Optimization; Technical Report TIK-112; Swiss Federal Institute of Technology: Zurich, Switzerland, 2001. (26) Schott, J. R. Fault tolerant design using single and multicriteria genetic algorithm optimization. Master’s Thesis, MIT, Cambridge, MA, 1995. (27) Levine, D. Users Guide to the PGAPack Parallel Genetic Algorithm Library; Technical Report ANL-95/18; Argonne National Laboratory: Argonne, IL, 1996. (28) MPICHsA Portable Implementation of MPI; Argonne National Laboratory: Argonne, IL, 2003 (available at http://wwwunix.mcs.anl.gov/mpi/mpich/). (29) Guthrie, K. M. Chem. Eng. 1969, 76, 114. (30) Young, D. M.; Cabezas, H. Comput. Chem. Eng. 1999, 23, 1477-1491. (31) Martin, T.; Young, D. Users Guide for Waste Reduction Algorithm Graphics User Interface. U.S. Environmental Protection Agency: Cincinnati, OH, 2001. (32) Diwekar, U. M. Batch DistillationsSimulation, Optimal Design, and Control; Taylor & Francis: Washington, DC, 1996. (33) Kim, K.-J.; Diwekar, U. M. Rev. Chem. Eng. 2001, 17, 111164.

Received for review December 16, 2003 Revised manuscript received February 20, 2004 Accepted March 3, 2004 IE0343162