Multiobjective Operation Optimization of Naphtha Pyrolysis Process

Subscriber access provided by University | of Minnesota Libraries

Article

Multi-objective Operation Optimization of Naphtha Pyrolysis Process using a Parallel Differential Evolution Xianpeng Wang, and Lixin Tang Ind. Eng. Chem. Res., Just Accepted Manuscript • DOI: 10.1021/ie401954d • Publication Date (Web): 11 Sep 2013 Downloaded from http://pubs.acs.org on September 12, 2013

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Industrial & Engineering Chemistry Research is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 30

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

Multi-objective Operation Optimization of Naphtha Pyrolysis Process Using a Parallel Differential Evolution

Xianpeng Wang, Lixin Tang *

Liaoning Key Laboratory of Manufacturing System and Logistics, The Logistics Institute, Northeastern University, Shenyang, 110819, China

*

Corresponding author. E-mail address: [email protected] (L. Tang); Tel/Fax: 86-24-83680169 1

ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Abstract A steam cracking furnace plays an important role in the naphtha pyrolysis process that produces ethylene and propylene. These two products are the key monomers in petrochemical industry. This paper investigates the multi-objective operation optimization of the naphtha pyrolysis process in order to maximize the yields of ethylene and propylene. To solve this problem efficiently, a multi-objective parallel differential evolution with competitive evolution strategies (MOPDE-CES) is presented. In the MOPDE-CES, the population is divided into four parts, each part having its own evolution strategy. During evolution, the four parts of the population compete with each other to survive and evolve. At the same time, they share search information with each other by communicating with an external archive which stores the obtained non-dominated solutions. The MOPDE-CES is compared to other powerful multi-objective algorithms, and computational results show that the MOPDE-CES is superior based on benchmarking and practical problems. In addition, the computational results of MOPDE-CES indicate that the operation of an ethylene plant can be improved by increasing the yields of ethylene and propylene. Keywords: Operation optimization; naphtha pyrolysis; multi-objective parallel differential evolution.

2


Page 2 of 30

Page 3 of 30

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


1. Introduction Ethylene is one of the most important monomers because it is the base for organic material production in the petrochemical industry. Although light hydrocarbons have been the major raw material for ethylene production in America and Middle East, in Asia the naphtha pyrolysis mechanism is the most important process in ethylene production.1 Another important product in the naphtha pyrolysis process is propylene, however, its yield is much smaller than the ethylene’s yield. Due to the dramatically increasing demand for propylene’s derivative products, how to increase its yield is now attracting more attentions. The key equipment of the naphtha pyrolysis mechanism in an ethylene plant is the steam cracking furnace. The pyrolysis process in the cracking furnace is shown in Figure 1. The liquid naphtha is first fed into the convection section to be preheated, and then a hot dilution steam is added into the convection section and mixed with the heated naphtha feed. In the convection section, the mixture of naphtha feed and dilution steam is completely vaporized. The vaporized mixture subsequently enters the cracking coil which is laid out in a U-type manner in the radiant heating section. The pyrolysis reaction occurs immediately along the cracking coil. The products of the reaction flow out to downstream equipment in which each product constituent is separated.

Figure 1. Illustration of the naphtha pyrolysis process in the steam cracking furnace

In the ethylene plant, the steam cracking furnace has the largest production capacity and consequently has the highest energy consumption. The operation level of the furnace can significantly impact the production cost of the whole ethylene plant and thus even a little 3



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

improvement in its operation can increase the economic benefits.2 In the practical operation of the steam cracking furnace, there are four key control variables: coil outlet temperature (COT), steam-to-naphtha ratio, naphtha feed flow (FLOW), and coil outlet pressure (COP). These control variables together determine the yield rates of ethylene and propylene. The objectives of the operation optimization of naphtha pyrolysis mechanism is to maximize the yield rates of ethylene and propylene. However, the two optimization objectives cannot be achieved simultaneously because the yields of ethylene and propylene are in conflict with each other. Increasing ethylene’s yield will reduce propylene’s yield. So this problem is a multi-objective optimization problem, in which there exists strong nonlinearity because the reaction network is very complex. In light of these observations, an algorithm should be able to solve this problem quickly and at the same time provide a good spread of non-dominated solutions for the decision maker to choose. Due to its great importance in reducing the production cost, the operation optimization of naphtha pyrolysis process has drawn continuing attentions from researchers. However, in literature there are relatively few research papers that deal with this problem. Most research efforts are devoted to analyzing the pyrolysis reaction mechanism. Once activated, thousands of reactions occur concurrently in the cracking coil in a very short time, consequently there is no reaction kinetics model that can precisely describe the pyrolysis process. In literature, many kinds of reaction kinetics models are proposed for the naphtha pyrolysis mechanism, and they can be classified into three categories according to the review of Gao et al.3: free-radical reaction kinetics models, molecular reaction kinetics models, and experiential models. Among the kinetics models, the free-radical reaction kinetics model4,5 has the most precise prediction of product yields; however, solving this model requires a great deal of computational resources, which in turn makes its application very limited. On the other hand, experiential models6 require less computational resources, but they are based on experience and can only be applied to specific naphtha feeds and furnaces. Therefore, the molecular kinetics models are most frequently used due to their simplicity, less demand for computational resources, and high preciseness of prediction. For more details on molecular kinetics models, the reader can refer to Hirato and Yoshioka7, Van Damme et al.8, and Kumar and Kunzru9. 4


Page 4 of 30

Page 5 of 30

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


As mentioned above, most research efforts are centered on the mechanism modeling of the naphtha pyrolysis process, but few papers focus on its operation optimization. Tarafder et al.2 investigated the operation optimization of an ethane cracking furnace and proposed a multi-objective optimization model for this problem. The non-dominated sorting genetic algorithm (NSGA-II) proposed by Deb et al.10 was used to solve this problem. Li et al.11 dealt with operation optimization of the naphtha pyrolysis process and adopted a multi-objective particle swarm optimization (MOPSO) to achieve the trade-off between the yields of ethylene and propylene. Gao et al.3 further considered the coke process and proposed a periodic operation optimization model for the naphtha pyrolysis process. The successive quadratic programming (SQP) was incorporated into the NSGA-II to construct a parallel algorithm to efficiently solve this multi-objective optimization problem. Nabavi et al.12 proposed an improved version of NSGA-II, named NSGA-II-aJG, in which a jumping gene operator was incorporated for the multi-objective operation optimization of an industrial LPG thermal cracker. Later Nabavi et al.13 further tackled the multi-objective design optimization of the LPG thermal cracker. In recent years, the research on multi-objective evolutionary algorithms (MOEAs) has made considerable progress. Many new MOEAs have been proposed and show improved performances in comparison to previous canonical MOEAs; such as NSGA-II and MOPSO14. The applications of MOEAs in many industries are very successful and have brought great benefits for many enterprises15-20. For the multi-objective operation optimization problem (MOOP) in chemical process, the MOEAs such as NSGA-II and MOPSO are often used by many researchers3,11,21,22. Differential Evolution (DE) is a new member of the swarm-based evolutionary algorithms, but it has become one of the most powerful stochastic search algorithms.23 Due to its simplicity and ease of implementation, DE has performed well in many single-objective optimization problems including the optimization problems in chemical processes24,25. Although few research efforts are devoted to designing the multi-objective DE (MODE), computational results reported in published papers show that MODE is very promising and that its performance is superior to NSGA-II and MOPSO26-28. The main reason for this superiority is that MODE has a better convergence speed than NSGA-II and a better diversity maintenance than MOPSO. The application of MODE in the 5



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

operation optimization in the electrical industry and the chemical industry also reveals that this algorithm is very powerful and effective29-31. Therefore, in this paper we adopt the MODE to solve the operation optimization of the naphtha pyrolysis process. The proposed MODE is different from the ones in literature because we adopt different kinds of evolutionary strategies to improve the algorithm’s robustness. This method has a parallel structure because the whole swarm is divided into several parts and each part has its own evolution strategy. During evolution, these strategies compete with each other to improve their evolutionary opportunities. The allocation of evolving opportunities is based on the suitability and contribution of each strategy and is implemented by the dynamical adjustment of each part’s size. The competition makes the algorithm adaptively select the most appropriate strategy, which in turn ensure the search robustness. In addition, these strategies also cooperate with each other by sharing their search information extracted from the external archive. This guarantees that the algorithm can converge quickly to the Pareto optimal solutions. We refer to the proposed parallel MODE in this paper as the parallel MODE with competitive evolution strategies (MOPDE-CES). One major concern in the design of the MOPDE-CES is the calculation of the objective values, i.e., the yields of ethylene and propylene. Since the pyrolysis process is complicated, one simulation of this reaction network will requires a significant computational effort, for example, solving ordinary differential equations based on the molecular kinetics model of Kumar and Kunzru9 takes about 6 seconds.32 In addition, the molecular kinetics models are constructed based on the assumption that the reaction is in a steady state. However, in reality production is dynamic because the source feed often changes and the operating conditions are constantly changing, so the adopted algorithm should provide the recommended operational settings as quickly as possible. Therefore, the steady molecular kinetics models have some limitations especially when applied to dynamic operating conditions. In a practical operation process, a large amount of production data is stored in the database, this provides a good selection on which to build data-based models that will predict the yields of ethylene and propylene for a given set of operation variables. The search results of Bhutani et al.18 have shown that the adoption of data-based prediction models in the operation optimization of chemical processes has many advantages over the first principle models. So in this paper, we 6


Page 6 of 30

Page 7 of 30

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


adopt the data-based model to predict the objective values (i.e., yields of ethylene and propylene) for a given set of decision variables. We chose the least square support machine (LSSVM) hybrid with a particle swarm optimization in our method. Our choice is different from the artificial neural network (ANN) that was adopted by Li et al.11 and Bhutani et al.18. The reason for choosing the LSSVM is that it has a better generalization ability and needs much less sample data than the ANN. Since the performance of the LSSVM is dependent on parameters of the adopted kernel functions, we take the parameter determination as an optimization problem and proposed a particle swarm optimization method to solve it. The structure of this paper is organized as follows. The procedure of the MOPDE-CES is described in Section 2. In Section 3, the computational experiments are carried out to evaluate the performance of the proposed MOPDE-CES. It is compared to the other powerful MOEAs using benchmark multi-objective optimization problems. In addition, a case study is provided to illustrate its efficiency in the application of the operation optimization of naphtha pyrolysis process, in which the LSSVM model is introduced and tested. This paper is finally concluded in Section 4. 2. Multi-objective parallel differential evolution with competitive evolution strategies Instead of finding a single optimal solution in the single objective optimization problem, the task of multi-objective optimization problem (MOP) is to obtain a set of Pareto optimal solutions with a good spread. For two solutions x and y in a MOP to minimize objectives, x is said to dominate y if and only if each objective of x is smaller or equal to that of y and at the same time there exists at least one objective of x that is smaller than that of y. If they cannot dominate each other, then they are called non-dominated solutions. For a MOP, if there exists no other solution x' that can dominate x, then x is called a Pareto optimal solution. All the Pareto optimal solutions are referred to as the Pareto set and their corresponding objective vectors in the objective space are called the Pareto front. When designing multi-objective algorithms to solve MOPs, four important tasks should be considered: (1) the Pareto front obtained by the algorithm should be close enough to the true Pareto front, (2) the Pareto front obtained by the algorithm should converge as quickly as possible to the true Pareto front, (3) the distribution of the non-dominated solutions obtained

7



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

by the algorithm should be as uniform as possible, and (4) the Pareto front obtained by the algorithm should be able to cover as many regions as possible of the true Pareto front, especially its extreme points. Begin: Initialization: 1. Set the iteration number i =0 and the external archive A to be empty. 2. Randomly generate the swarm with N solutions according to the diversity method of Nebro et al. (2008), divide it into k (number of evolution strategies) parts with equal size, and then assign an evolution strategy to each part Pj, j=1, …, k. 3. Assign each part with a given evolution strategy to a processor. 4. Evaluate the solutions in each part Pj, j=1, …, k. 5. Add the nondominated solutions in the whole swarm into A. while (the termination criterion is not reached) do 1. Set i = i +1. 2. Update each part Pj, j=1, …, k, with the given differential evolution strategy. 3. Evaluate the solutions in each part Pj, j=1, …, k. (Note that the idle processors will wait until all parts finish the evolution.) 4. After the evaluation of all parts, update A with the new nondominated solutions in the whole swarm. 5. Perform local search to A 6. If i is a multiple of ncompete (frequency of competition) 6.1 Calculate the contribution of each evolution strategy, and then determine the new size of each part Pj, j=1, …, k. 6.2 Adjust each part according to its new size. 7. Share the search information between different parts. End while Report the obtained non-dominated solutions in A. End Figure 2. The overall framework of the proposed MOPDE-CES

8


Page 8 of 30

Page 9 of 30

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Start

Set the generation counter i=0 and the external archive A to be empty

Generate the initial population

Divide the population into k parts and assign an evolution strategy to each part

Evaluate the solutions in each part and then update the external archive A

Set i=i+1

Allocate the evolution of each part to a different processor

Update each part with the given differential evolution strategy

Idle processors wait

N

All parts have been updated Y

Evaluate the solutions in each part and then update the external archive A

Apply local search on A

N

i is a multiple of ncompet e Y

Calculate the new size of each part and then adjust each part based on its new size

Share the search information between parts

Termination criterion is reached

N

Y Output the obtained Pareto front

Stop

Figure 3. The flowchart of the proposed MOPDE-CES

9



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 30

With these tasks in mind, we developed the MOPDE-CES algorithm. With respect to the first two tasks relating to the exploitation ability, the MOPDE-CES adopts multiple evolution strategies that compete and cooperate with each other during evolution. Their competition ensures that most evolution opportunities are adaptively allocated to only evolutionary strategies that are most suitable for the current MOP. On the other hand, their cooperation through sharing search information can ensure that the search direction is guided to promising regions. The competition and cooperation together guarantee a very close distance and quick convergence speed of the obtained Pareto front to the true Pareto front. For the last two tasks concerning the exploration (diversity) ability, the MOPDE-CES divides the whole swarm into several parts and assigns an evolution strategy to each part. Since each evolution strategy has its own search behavior, the diversity of the obtained Pareto front can be ensured. In addition, the parallel computation is adopted by allocating each evolution strategy to a different processor since the CPU generally has multiple processors in practical industry. The overall framework and the flowchart of the MOPDE-CES are given in Figure 2 and Figure 3, respectively. In the following, each component of the MOPDE-CES are described in details. 2.1 Update of each part with different evolution strategies For a target solution Xi with dimension n, the generation of a new solution consists of three consecutive steps: mutation, crossover, and selection. There are several mutation strategies in literature and each of them has its own distinct suitability for different MOPs. A so called SaDE algorithm was developed by Qin et al.33 in which different kinds of mutation strategies are combined together and SaDE adaptively selects the most appropriate strategies to generate new solutions. Motivated by the promising results obtained by SaDE, we extended this idea to a multi-objective scenario and made some modifications such as the adoption of multiple swarms and the competition and cooperation of different strategies to achieve further improvements. 2.1.1 Mutation In our MOPDE-CES, the swarm is divided into four parts (i.e., k=4) and each part is assigned a mutation strategy from a set of four strategies. For simplicity, these strategies, namely DE/rand/1, DE/best/1, DE/rand-to-best/1, and DE/best/2, are assigned to part 1, part 2, part 3, and part 4, respectively. The definitions of these strategies are as follows. 10


Page 11 of 30

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


DE/rand/1: Vi = Xr1 + F  (Xr2 – Xr3);

(1)

DE/best/1: Vi = Xbest + F  (Xr1 – Xr2);

(2)

DE/ rand-to-best /1: Vi = Xi + F  (Xbest – Xi) + F  (Xr1 – Xr2);

(3)

DE/best/2: Vi = Xbest + F  (Xr1 – Xr2) + F  (Xr3 – Xr4).

(4)

In the above equations, Vi is called the perturbed or mutated vector, F is the control parameter, and Xr1, Xr2, Xr3, Xr4 are mutually different solutions selected from each part (note that they are also different from Xi). 2.1.2 Crossover After the perturbed vector Vi is generated in the mutation step, the crossover operator is performed on Vi = (v1,i, v2,i, …, vn,i) and the target solution Xi = (x1,i, x2,i, …, xn,i) to generate a v j ,i , if rand j  Cr or j  jrand trial vector Ui = (u1,i, u2,i, …, un,i) according to u j ,i   for  x j ,i , otherwise

j=1, …, n, where randj is a random number in [0, 1], jrand is a random integer number in [1, n] to ensure that the trial solution Ui is different from Xi for at least one dimension, and Cr is the crossover probability. 2.1.3 Selection The selection step is used to guarantee the elitism of evolution process, that is, if the trial vector Ui dominates Xi or they are non-dominated, then the target vector Xi is replaced by Ui; otherwise, Xi will survive to the next generation. 2.2 Competition of different evolution strategies The competition of different evolution strategies is illustrated by the size adjustment of parts they are assigned to. As shown in Figure 2, this procedure is performed every ncompete generations instead of every generation because a sufficient number of evolutions should be carried out by the strategies before their performances are evaluated. Once called by the algorithm, this procedure first determines the new size of each part based on its contribution to the obtained external archive A. The contribution of a part Pj (or the evolution strategy assigned to it) is defined as the percentage of non-dominated solutions obtained by this part in the external archive A (denoted as pj). Then the new size of this part is calculated by  N  p j  . If the total sum of new sizes is less than N, then the smallest size is 11



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 30

increased to ensure that there is a total of N solutions in the whole swarm. Given the new size of each part, we then begin to adjust the solutions in each part. For a part Pj, if its size decreases, we first sort the solutions in this part using the fast non-dominated sorting method of NSGA-II and then iteratively removing the worst solution until the new size is reached. On the contrary, if its new size increases, we first randomly select a solution from this part and another one from A, and then generate a new solution using the SBX operator34 based on them. This generation of new solutions will be iteratively performed until the new size is reached. 2.3 Cooperation of different evolution strategies The cooperation of different evolution strategies is implemented in two ways. The first one lies in the mutation procedure. For example, in the last three strategies, Xbest is not the best solution in each part, but a non-dominated solution that is randomly selected from A. Since Xbest may be obtained by the other strategies, the search information can be shared by different strategies. The second way of sharing search information shown in Figure 2 is performed as follows. At each generation, we randomly select two parts, then randomly select a good solution in each part, and finally exchange them. This exchange is carried out twice for each generation. Through this exchange, the good search information can be transferred between different parts, which in turn can help to improve the convergence speed of our algorithm and the quality of obtained solutions. 2.4 Update of the external archive As did in many MOEAs, we use the method of NSGA-II to update the external archive. That is, if a new solution X is dominated by a certain solution in the external archive, then solution X is discarded. On the contrary, if solution X is not dominated by any solution in the external archive, then solution X will be inserted into the external archive. After the insertion, all the solutions that are dominated by solution X will be removed from the external archive. In addition, if no solution is removed and the size of the external archive exceeds its limit, we will remove the most crowded solution from the external archive. 2.5 Local search to improve the external archive To further improve the convergence speed and at the same time maintain diversity, a 12


Page 13 of 30

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


local search is performed after the external archive A is updated by the new swarm. The local search is based on the search mechanism of genetic algorithm, which can inherit the good characteristics of parents and at the same time maintain good diversity. In the local search, a total of |A| new solutions are generated, and the generation of each new solution is as follows. First, we randomly select two non-dominated solutions from A. Next, we generate a new solution by performing the SBX operator on them. Last, the external archive is updated with these |A| new solutions. Since the adjustment of each part in the swarm is based on its contribution to A, we should give a mark to each new solution so as to calculate the contribution of each part. In our algorithm, this mark on a new solution is randomly selected to help improve the search diversity of the whole swarm. 2.6 Handling constraints In most practical production processes, the operation optimization problem generally contains some constraints, e.g., constraints on control variables, constraints on production process, and so on. To handle these constraints, we adopt the method proposed by Deb et al.10. For two solutions X1 and X2, X1 is said to constrained-dominate X2 if any of the following three conditions is met: (1) both solutions are feasible and X1 dominates X2; (2) solution X1 is feasible but solution X2 is infeasible; and (3) both solutions are infeasible and the constraint violation of X1 is smaller than that of X2. 3. Case studies To test the performance of the proposed MOPDE-CES, two case studies were carried out. The first one is based on the benchmark MOPs selected from the literature, and the second one is the operation optimization of naphtha pyrolysis process. The setting of parameters is as follows: the size of the whole swarm is 100, the control parameter F=0.5, the crossover probability Cr=0.1, and the frequency of competition ncompete = 10. To avoid the scenario that the size of one part takes an overwhelming majority, each part has a minimum size of 10. Since in practical operation it is commonly required that the optimization algorithm should converge as quick as possible, the stopping criterion is set to a relatively small number of function evaluations, i.e., a maximum of 15000 function evaluations. Please note that this stopping criterion includes all the relative function 13



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 30

evaluations in the algorithm, i.e., swarm update and local search. 3.1 Performance metrics The performances of the tested algorithms are evaluated by two metrics that are often adopted in the literature. The first one is the general distance (GD) metric which is used to evaluate how close the obtained Pareto front is to the true Pareto front. The GD is defined as GD 



A i 1

di2 A , where |A| is the size of the external archive and di is the Euclidean

distance in the objective space between every solution in the obtained Pareto front and its nearest solution in the true Pareto front. The second one is the hyper-volume (HV) metric that can evaluate not only the distance between the obtained Pareto front and the true Pareto front but also the distribution of the obtained Pareto front. The value of HV is calculated as the cumulative area covered by a set of non-dominated solutions in comparison to a given reference point R. An illustration of the HV for a bi-objective minimization optimization problem is shown in Figure 4. In the calculation of HV, the jth objective value of the Pareto front obtained by each algorithm is normalized by f j '  ( f j  f jmin ) ( f jmax  f jmin ) where f jmax and f jmin are respectively the maximum and minimum values of the jth objective of the true Pareto front. The reference point is then selected as (1, 1) for the bi-objective MOPs and (1, 1, 1) for the three-objective MOPs. Please note that a large value of HV results in a better performance. f2

Reference point

nondominated solution f1 Figure 4. Illustration of the calculation of Hypervolume metric

3.2 Case study I: benchmark multi-objective optimization problems As shown in Table 1, eight difficult benchmark multi-objective optimization problems 14


Page 15 of 30

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


are selected from the literature to compare the MOPDE-CES with three powerful MOEAs such as MOEA/D proposed by Zhang and Li35, NSGA-II and SMPSO proposed by Nebro et al.36. Please note that SMPSO is shown to be superior to previous MOPSO proposed by Coello Coello et al.14. The source codes of the three rival algorithms are implemented in Java using jMetal, a framework that can be downloaded from http://jmetal.sourceforge.net/. Table 1. Definitions of test problems used in the experiment Problem

Definition

ZDT3 (Zitzler et al.37)

Kita (Kita et al.38)

Constraints

f1 ( x)  x1 , f 2 ( x)  g ( x)[1  x1 g ( x)  x1 sin(10 x1 ) g ( x)], g ( x)  1  9 i  2 xi (n  1) n

f1 ( x, y)  x  y, f 2 ( x, y)  0.5x  y  1

x/6 + y – 6.5 ≤ 0, x/2 + y – 7.5 ≤ 0, 5x + y – 30 ≤ 0, 0 ≤ x, y ≤ 7

f1 ( x)  x1 , f 2 ( x)  (1  x2 ) x1

9x1+ x2 ≥ 6, 9x1– x2 ≥ 1, 0.1 ≤ x1≤ 1.0, 0 ≤ x2≤ 5

2

Constr (Deb et al.10)

n = 30, 0 ≤ xi ≤ 1, i = 1,…, n

f1 ( x)  (1  g ( x)) cos(0.51 ) cos(0.5 2 ), f 2 ( x)  (1  g ( x)) cos(0.51 )sin(0.5 2 ),

DTLZ6 (Deb et al.39)

n = 12, 0 ≤ xi ≤ 1, i = 1,…,n

f3 ( x)  (1  g ( x))sin(0.51 ), g ( x)   i 3 xi0.1 , n

i   (1  2 xi g ( x)) (4(1  g ( x)))

Viennet4 (Viennet et al.40)

f1(x, y) = (x – 2)2/2 + (y +1)2/13 + 3, f2(x, y) = (x + y – 3)2/175 + (2y – x)2/17 – 13, f3(x, y) = (3x –2y + 4)2/8 + (x – y +1)2/27 + 15

y < –4x + 4, x > –1, y > x – 2, –4 ≤ x, y ≤ 4

Binh4 (Binh Ulrich41)

f1(x, y) = 1.5 – x(1– y), f2(x, y) = 2.25 – x(1– y2), f3(x, y) = 2.625 – x(1– y3)

–x2–(y–0.5)2+9≤0, (x–1)2+(y–0.5)2 –6.25≤0, –10 ≤ x, y ≤ 10

f1(x, y, z) = x, f2(x, y, z) = y, f3(x, y, z) = z

x2 + y2 + z2 ≤ 1, 0 ≤ x, y, z ≤ 1

and

Tamaki (Tamaki et al.42)

LZ09_F6 (Li and Zhang43)

f1 ( x)  cos(0.5 x1 ) cos(0.5 x2 ) 

2 2 x j  2 x2 sin(2 x1  j n)  ,   jJ1 | J1 |

f 2 ( x)  cos(0.5 x1 )sin(0.5 x2 ) 

2 2   x j  2 x2 sin(2 x1  j n)  , | J 2 | jJ 2

f3 ( x)  sin(0.5 x1 ) 

2 2   x j  2 x2 sin(2 x1  j n)  , | J 3 | j J 3

J1   j | 3  j  n, and j  1 is a multiplication of 3 ,

n = 10, 0 ≤ x1 , x2 ≤ 1, –2 ≤ xi ≤ 2, i = 3,…, n

J 2   j | 3  j  n, and j  2 is a multiplication of 3 , J 3   j | 3  j  n, and j is a multiplication of 3

In the experiments, we run each algorithm for 100 duplications of each test problem and collect the average result and the interquartile range (IRQ) as measures of location and 15



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 30

statistical dispersion. In the comparison, a better value of the average result is preferred. 3.2.1 Efficiency of adopting the competitive evolutionary strategies In this subsection, we test the efficiency of the proposed competitive evolutionary strategies. We modified the MOPDE-CES by allowing only a single evolutionary strategy and by removing the dynamic control of each part’s size, and subsequently we obtained four variations, namely MOPDE-S1 (with only DE/rand/1), MOPDE-S2 (with only DE/best/1), MOPDE-S3 (with only DE/ rand-to-best /1), and MOPDE-S4 (with only DE/best/2). The comparison results of GD and HV metrics are shown in the following Table 2 and Table 3, respectively. In the tables, the better results are shown in bold style, and the last column is the significance result to show that whether the performance difference between our MOPDE-CES and the best one among the others is significant (denoted as “+”) or not (denoted as “–”) with a confidence interval of 95%. Table 2. Comparison results of GD metric between MOPDE-CES and other variations Problems

MOPDE-CES

MOPDE-S1

MOPDE-S2

MOPDE-S3

MOPDE-S4

sig

average

IQR

average

IQR

average

IQR

average

IQR

average

IQR

ZDT3

3.358e-05

3.341e-06

9.495e-05

6.882e-05

3.477e-05

4.043e-06

1.007e-04

8.524e-05

3.562e-05

4.288e-06

Kita

3.807e-04

1.344e-03

3.710e-03

3.910e-03

1.936e-03

1.617e-03

3.966e-04

1.866e-03

2.138e-03

1.927e-03

Constr

2.761e-04

6.678e-05

6.795e-04

2.277e-04

6.272e-04

1.959e-04

6.272e-04

1.321e-04

6.535e-04

2.167e-04

DTLZ6

5.718e-04

3.676e-05

1.254e-03

3.341e-05

6.368e-04

3.030e-05

2.630e-02

3.303e-02

7.925e-04

2.886e-05

Viennet4

3.563e-04

1.292e-04

4.715e-04

1.605e-04

4.444e-04

1.559e-04

4.469e-04

1.368e-04

4.477e-04

1.368e-04

+ + + +

Binh4

1.168e-03

2.629e-04

2.322e-03

3.168e-04

1.180e-03

2.590e-04

1.213e-03

2.553e-04

1.225e-03

3.107e-04

–

Tamaki

2.083e-03

3.840e-04

4.169e-02

3.880e-02

2.201e-03

3.524e-04

1.887e-03

2.508e-04

2.314e-03

3.919e-04

+

LZ09_F6 3.581e-02 3.364e-02 4.424e-02 2.969e-02 4.548e-02 3.437e-02 3.925e-02 2.908e-02 3.630e-02 2.725e-02

–

–

Table 3. Comparison results of HV metric between MOPDE-CES and other variations Problems

MOPDE-CES

MOPDE-S1

MOPDE-S2

MOPDE-S3

MOPDE-S4

sig

average

IQR

average

IQR

average

IQR

average

IQR

average

IQR

ZDT3

0.5876

2.181e-05

0.5847

2.066e-03

0.5872

6.795e-05

0.5789

5.154e-03

0.5875

7.333e-05

Kita

0.6525

4.178e-04

0.6465

2.678e-03

0.6506

8.001e-04

0.6508

6.797e-04

0.6498

1.002e-03

Constr

0.7749

1.583e-03

0.7679

4.471e-03

0.7689

3.453e-03

0.7722

2.929e-03

0.7698

4.278e-03

DTLZ6

0.0949

4.069e-05

0.0949

6.959e-05

0.0949

4.101e-05

0.0657

2.593e-02

0.0949

3.792e-05

Viennet4

0.8575

1.779e-03

0.8568

1.883e-03

0.8575

1.618e-03

0.8575

2.409e-03

0.8574

1.859e-03

Binh4

0.3565

8.304e-04

0.3541

1.878e-03

0.3557

9.374e-04

0.3559

8.183e-04

0.3556

1.154e-03

– – – – – –

Tamaki

0.4525

3.519e-03

0.2507

3.524e-02

0.4509

3.414e-03

0.4508

4.458e-03

0.4499

3.567e-03

+

LZ09_F6

0.2368

3.751e-02

0.2332

3.158e-02

0.2369

3.290e-02

0.2310

3.414e-02

0.2332

3.348e-02

–

Based on the results, it appears that the adoption of competitive evolution strategies can combine the advantages of different strategies and thus the MOPDE-CES has a much better 16


Page 17 of 30

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


performance and robustness than other MOPDEs that only adopts a single evolution strategy. 3.2.2 Comparison with other powerful MOEAs The comparison results are given in Table 4 and Table 5 for the GD metric and the HV metric. In addition, to give a visualization view of these results, we show the graphic results for the problem Kita in Figure 5, the problem DTLZ6 in Figure 6, the problem Binh4 in Figure 7, and the problem Tamaki in Figure 8, respectively. In each figure, the result of each algorithm corresponds to the one with the best GD metric (the results of MOEA/D are not included in these figures because the results have shown that the MOPDE-CES is superior to MOEA/D). Table 4. Comparison results of GD metric between MOPDE-CES and other MOEAs Problems

MOPDE-CES

MOEA/D

SMPSO

NSGA-II

average

IQR

average

IQR

average

IQR

average

IQR

ZDT3

3.358e-05

3.341e-06

2.761e-03

6.339e-04

2.447e-04

8.397e-04

2.484e-04

4.272e-05

Kita

3.807e-04

1.344e-03

1.484e-01

3.603e-05

2.613e-01

4.175e-03

1.247e-03

3.159e-03

Constr

2.761e-04

6.678e-05

3.131e-04

1.310e-07

3.346e-02

3.804e-04

3.033e-04

5.040e-05

DTLZ6

5.718e-04

3.676e-05

2.359e-03

4.282e-02

5.800e-04

5.033e-05

3.534e-01

2.718e-02

Viennet4

3.563e-04

1.292e-04

1.605e-02

2.183e-04

2.798e-02

3.290e-03

4.574e-04

1.837e-04

Binh4

1.168e-03

2.629e-04

9.976e+01

0.000e+00

1.150e+00

5.107e+00

8.170e-03

1.667e-04

Tamaki

2.083e-03

3.840e-04

7.482e-02

0.000e+00

7.482e-01

0.000e+00

3.598e-03

5.434e-04

LZ09_F6

3.581e-02

3.364e-02

7.271e-02

7.182e-03

3.860e+00

2.068e+00

1.951 e+00 1.752 e+00

sig

+ + + + + + + +

Table 5. Comparison results of HV metric between MOPDE-CES and other MOEAs Problems

MOPDE-CES

MOEA/D

SMPSO

NSGA-II

sig

average

IQR

average

IQR

average

IQR

average

IQR

ZDT3

0.5876

2.181e-05

0.2731

4.247e-02

0.5806

3.100e-02

0.5821

1.043e-03

Kita

0.6525

4.178e-04

0.0062

0.000e+00

0.5277

3.961e-01

0.6482

6.496e-04

+ +

Constr

0.7749

1.583e-03

0.5583

1.420e-05

0.7490

8.079e-03

0.7744

3.645e-04

–

DTLZ6

0.0949

4.069e-05

0.0225

9.806e-03

0.0949

1.421e-04

0.7743

0.000e+00

Viennet4

0.8575

1.779e-03

0.3734

1.463e-02

0.7062

1.766e-01

0.8552

2.457e-03

Binh4

0.3565

8.304e-04

0.3787

0.000e+00

0.6151

0.000e+00

0.4425

2.067e-03

Tamaki

0.4525

3.519e-03

0.3787

0.000e+00

0.0000

0.000e+00

0.4339

5.139e-03

LZ09_F6

0.2368

3.751e-02

0.0000

0.000e+00

0.0000

0.000e+00

0.0000

0.000e+00

+ + + + +

Based on these results, it appears that our MOPDE-CES obtains the best GD results for all the eight test MOPs and that its performance is significantly better than MOEA/D, NSGA-II and SMPSO. For the HV metric, our MOPDE-CES obtains the best results for six out of the eight MOPs. Although the NSGA-II obtains the best result for problem DTLZ6, from Figure 6 we can see that the reason for this is that the NSGA-II algorithm does not 17



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 30

converge to the true Pareto front which is a curve. Therefore, our MOPDE-CES is still superior to NSGA-II for DTLZ6. From the graphical results, it can be seen that SMPSO cannot reach the true Pareto front for Kita, and it is actually trapped in local optimum for problems Binh4 and Tamaki. The reason is that PSO has a very quick convergence speed but a very weak ability to get out of local optimum. Although NSGA-II has a better performance than SMPSO, it fails to reach the true Pareto front for DTLZ6, and the distribution of the Pareto fronts obtained by it, is inferior to that obtained by our MOPDE-CES. Therefore, it can be concluded that the proposed MOPDE-CES outperforms the MOEA/D, the NSGA-II and the SMPSO for almost all the test MOPs.

Figure 5. Comparison of Pareto fronts between different algorithms for Kita

Figure 6. Comparison of Pareto fronts between different algorithms for DTLZ6 18


Page 19 of 30

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 7. Comparison of Pareto fronts between different algorithms for Binh4

Figure 8. Comparison of Pareto fronts between different algorithms for Tamaki

19



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 30

3.3 Case study II: multi-objective operation optimization of naphtha pyrolysis process In this Subsection, the operation optimization of naphtha pyrolysis process is investigated. In the ethylene plant considered in this paper, there are three main control variables, namely the COT, the COP, and the naphtha feed flow (FLOW), which have important impact on the production process control. In the practical production of the ethylene plant considered in this paper, the steam-to-naphtha ratio is fixed. It should be noted that, although the coking phenomena is relevant to the yield rates of ethylene and propylene, the coking process is not considered in this paper due to the following two reasons: (1) The aim of this paper is to investigate the multi-objective operation optimization and construct a LSSVM model based on practical data to predict the yields of ethylene and propylene, so the input control variables of the model should be measurable. However, the coke thickness cannot be obtained from the database of the practical production because the coke thickness cannot be measured during production process, and consequently it is impossible to take the coke thickness as one of the input variables in the prediction model. (2) For a given naphtha whose physical property is fixed, the coke thickness is mainly determined by the depth of cracking, the cracking temperature, the dilution steam, and the retention time.44 These factors can be seen as a function of the above three control variables. Therefore, as did in Li et al.11 that adopted the artificial neural network as the prediction model, the coke thickness is not defined as a control variable in this paper. 3.3.1 LSSVM for the prediction of yields of ethylene and propylene As mentioned in the introduction section, we prefer to use a LSSVM to construct the mathematical relationship between the input control variables (COT, COP, and FLOW) and the two output objectives (the yields of ethylene and propylene). Since there are two objectives, two LSSVM models are built for them, respectively. In the LSSVM model, the regularization parameter γ and the kernel width parameter σ2 are the key parameters.45 Since the two key parameters have significant impact on the performance of the LSSVM model, we take the parameter determination as an optimization problem that is defined as follows:

20


Page 21 of 30

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Minimize

1 n ( yˆ i  yi ) 2  n i 1

(5)

s.t. yˆi  f ( ,  2 , COT , COP, FLOW )

(6)

γmin≤ γ ≤ γmax

(7)

σ2min≤ σ2 ≤ σ2max

(8)

where yˆ i is the prediction value obtained by the LSSVM model that is determined by the parameters γ and σ2, yi is the ith sample data, n is the sum of sample data, and f(·) denotes the LSSVM model. In the experiment, we set γmin=0, γmax=1000, σ2min=0, and σ2max=100. To solve the above optimization problem, a scatter particle swarm optimization (PSO) algorithm proposed by Yin et al.46 is adopted. The data used in the construction of the LSSVM model are collected from the practical production process. The preprocessing of these data is carried out as follows: first, the k-mean clustering algorithm is used to check and remove outlier samples; second, the z-score normalization method is adopted to normalize the left samples; at last, we select 80 sample data to train the LSSVM and another 20 sample data to test the constructed LSSVM model. The training and test results based on the PSO algorithm are given in Figure 9, from which it can be found that the maximum prediction errors for both the training and test sets are less than 1.1%. Therefore, it is reasonable to substitute the constructed LSSVM models for the molecular kinetics model so as to obtain a quick response to the dynamic changes in practical production.

Figure 9. Training and test results based on the LSSVM model

21



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 30

3.3.2 Multi-objective operation optimization The task of the operation optimization of naphtha pyrolysis process is to determine the optimal values of control variables so as to maximize two conflicting objectives: the yield of ethylene and the yield of propylene. The mathematical model of this problem can be given as follows. min -YC2 H4

(9)

min -YC3H6

(10)

s.t. 9000 kg/h ≤ FLOW≤ 10666 kg/h

(11)

1093 K ≤ COT ≤ 1123 K

(12)

0.152 MPa ≤ COP ≤ 0.165 MPa

(13)

where YC2 H 4 and YC3 H6 are the yield of ethylene and the yield of propylene, respectively. 3.3.2.1 Sensitivity analysis of decision variables In this section, the sensitivity analysis of each decision variable is carried out by keeping the other decision variables constant. The effect of the naphtha feed flow on the yield of ethylene and the yield of propylene is shown in Table 6 and Figure 10, in which the outlet pressure and the outlet temperature are kept constant (COP = 0.155 MPa and COT = 1108 K). From the results, it can be seen that with the increase of the naphtha feed flow from 9000 kg/h to 10666 kg/h, the yield of ethylene decreases while the yield of propylene increases accordingly. Table 6. Naphtha feed flow to the yield of ethylene and the yield of propylene Naphtha feed flow (kg/h)

Yield of ethylene (%)

Yield of propylene (%)

9000.00

33.2405

16.1663

9416.50

33.1737

16.1676

9833.00

33.1075

16.1919

10249.50

32.9961

16.2163

10666.00

32.8469

16.2652

22


Page 23 of 30

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 10. Naphtha feed flow relative to the yield of ethylene and the yield of propylene

With the setting of FLOW = 9833 kg/h and COT = 1108 K, the results in Table 7 and Figure 11 show that, as the outlet pressure increases from 0.152 MPa to 0.165 MPa, the yield of ethylene first increases quickly but then decreases slowly, while the yield of propylene first decreases quickly but then increases slowly. The sum of the yields of ethylene and propylene also first increases and reaches its maximum value of about 49.41% when COP takes the value of about 0.1585 MPa, and then decreases slowly. Table 7. Outlet pressure to the yield of ethylene and the yield of propylene Outlet pressure (MPa)



0.1520

32.9657

16.2238

0.1553

33.1190

16.1869

0.1585

33.2530

16.1605

0.1617

33.2453

16.1629

0.1650

33.2276

16.1739

Figure 11. Outlet pressure relative to the yield of ethylene and the yield of propylene 23



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 24 of 30

Based on Table 8 and Figure 12, in which the naphtha feed flow and the outlet pressure are kept constant (i.e., FLOW = 9833 kg/h and COP= 0.16 MPa), it can be found that the yield of ethylene continuously increases from 33.0875% to 33.4345%, while the yield of propylene continuously decreases from 16.2911% to 16.1819%, as the outlet temperature increases from 1093 K to 1123 K. In addition, it appears that the increase of the yield of ethylene and the decrease of the yield of propylene are much faster when COT changes from 1093 K to 1108 K. Table 8. Outlet temperature to the yield of ethylene and the yield of propylene Outlet temperature (K)



1093.00

33.0875

16.2911

1100.50

33.2702

16.2347

1108.00

33.3875

16.1902

1115.50

33.4135

16.1837

1123.00

33.4345

16.1819

Figure 12. Outlet temperature relative to the yield of ethylene and the yield of propylene

3.3.2.2 Optimization results The computational results of the proposed MOPDE-CES and the NSGA-II, which has a wide application in chemical process optimization2,3,12,15,17,18, are given in Table 9 and Table 10. Since the true Pareto front of this problem is not known, we run both the MOPDE-CES and NSGA-II for a maximum of 50,000 function evaluations and then take the non-dominated solutions from the union of their obtained results as the reference Pareto front. Although a similar problem was handled by Li et al.11, we do not compare our method with the multi-objective PSO proposed in Li et al.11. The main reasons are as follows: first, the prediction model of Li et al.11 is the ANN model and the authors did not provide detailed 24


Page 25 of 30

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


parameter setting, which makes it very difficult to precisely re-implement this method, second, the multi-objective PSO used in Li et al.11 is very simple in comparison to the MOPSO proposed by Coello Coello et al.14 and the SMPSO proposed by Nebro et al.36. In Table 9, another metric is adopted to show the number of non-dominated solutions obtained by each algorithm (please note that in the experiment both algorithms take a maximum size of 100 for the external archive A). It is clear that a larger value of |A| is preferred because more candidate decisions can be provided to the decision-makers. Table 9. Comparison results for the operation optimization of naphtha pyrolysis MOPDE-CES

NSGA-II

Metrics

significance average

IQR

average

IQR

GD

2.657e-04

1.616e-05

3.039e-04

2.388e-05

+

HV

0.8268

1.521e-04

0.8249

3.327e-04

+

|A|

100.00

0.00

93.59

3.00

+

From Table 9, it is obvious that the proposed MOPDE-CES outperforms NSGA-II for all the performance metrics and that their performance differences are significant. Furthermore, the detailed results obtained by the two algorithms are presented in Table 10. Based on these results, it can be found that the proposed MOPDE-CES obtains a better distribution of both the yield rates and the decision variables than those obtained by the NSGA-II. With the help of our MOPDE-CES, the decision makers can have a wider range of candidate settings for the control variables. Table 10. Optimization results of NSGA-II and MOPDE-CES Metric

MOPDE-CES

NSGA-II

Min

Max

Min

Max

Ethylene

32.6199

33.5072

32.6284

33.4752

Propyle

16.1663

16.7740

16.1681

16.7719

FLOW

9000.00

10394.09

9833.00

10395.49

COP

1.5196

1.6500

1.5205

1.6500

COT

820.00

849.91

827.65

845.11

Yield rate

ne Decision variable

In the ethylene plant, the current yields of ethylene and propylene are 32.727% and 15.332%. With the MOPDE-CES, the yields of ethylene and propylene can be simultaneously

25



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 26 of 30

increased to about 33.43% and 16.19%, respectively. These improvements will be very important and beneficial to the ethylene plant. 4. Conclusions This paper studies the operation optimization of naphtha pyrolysis process to maximize the yields of ethylene and propylene, which usually conflict with each other in practical production. A multi-objective differential evolution algorithm is adopted to solve this problem because the differential evolution has the ability to balance convergence speed and search diversity. To make the proposed MODE more robust, multiple evolution strategies are adopted and they compete and cooperate with each other during the evolution strategies. A local search is also incorporated to further improve its performance. The computational results on both the benchmark problems and the practical operation optimization of naphtha pyrolysis show that the proposed MOPDE-CES can produce acceptable Pareto fronts with good distribution for all testing problems, and that its performance is significantly better than those of the MOEA/D, NSGA-II and SMPSO. In addition, since the proposed MOPDE-CES is a general algorithm for the multi-objective optimization, it may also be applied to other complicated multi-objective optimization problems in the petrochemical industry.

Acknowledgements This research is partly supported by State Key Program of National Natural Science Foundation of China (Grant No. 71032004), the Foundation for Innovative Research Groups of the National Natural Science Foundation of China (Grant No. 71321001), and the National Natural Science Foundation of China (Grant No. 70902065, 61374203).

Supporting Information Available This information is available free of charge via the Internet at http://pubs.acs.org/.

Reference (1) Nakamura, D. N. Global ethylene capacity increases slightly in 2006. Oil. Gas. J. 2007, 105, 46–60. (2) Tarafder, A.; Bennett, L.; Ray, A. K.; etc. Multiobjective optimization of an industrial 26


Page 27 of 30

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


ethylene reactor using a nondominated sorting genetic algorithm. Ind. Eng. Chem. Res. 2005, 44, 124−141. (3) Gao, X. D.; Chen, B. Z.; He, X. R.; etc. Multi-objective optimization for the periodic operation of the naphtha pyrolysis process using a new parallel hybrid algorithm combining NSGA-II with SQP. Comput. Chem. Eng. 2008, 32, 2801−2811. (4) Sundaram, K. M.; Froment, G. F. Modeling of the thermal cracking kinetics. 3. Radical mechanisms for the pyrolysis of simple paraffins, olefins, and their mixtures. Ind. Eng. Chem. Fund. 1978, 17, 174–182. (5) Joo, E.; Park, S. Pyrolysis reaction mechanism for industrial naphtha cracking furnaces. Ind. Eng. Chem. Res. 2001, 40, 2409−2415. (6) Hao, H.; Xiong, G. H.; Zhang, S. F.; etc. The calculation of propylene yield by the thermal cracking of hydrocarbon. Chem. Eng. 1999, 27, 54–56. (7) Hirato, M.; Yoshioka, S. Simulation of the pyrolysis of naphtha, kerosene and gasoil with a tubular reactor. Int. Chem. Eng. 1973, 13(2), 347−354. (8) Van-Damme, P. S.; Froment, G. F.; Balthasar, W. B. Scaling up of naphtha cracking coils. Ind. Eng. Chem. Process. Des. Dev. 1981, 20, 366–376. (9) Kumar, P.; Kunzru, D. Modeling of naphtha pyrolysis. Ind. Eng. Chem. Process. Des. Dev. 1985, 24, 774−782. (10) Deb, K.; Agrawal, S.; Pratap, A.; etc. A fast and elitist multi-objective genetic algorithm: NSGA-II. IEEE Trans. Evolut. Comput. 2002, 6, 182−197. (11) Li, C. F.; Zhu, Q. X.; Geng, Z. Q. Multi-objective particle swarm optimization hybrid algorithm: an application on industrial cracking furnace. Ind. Eng. Chem. Res. 2007, 46, 3602–3609. (12) Nabavi, R.; Rangaiah, G. P.; Niaei, A.; Salari, D. Multi-objective optimization of an industrial LPG thermal cracker using a first principles model. Ind. Eng. Chem. Res. 2009, 48, 9523−9533. (13) Nabavi, R.; Rangaiah, G. P.; Niaei, A.; Salari, D. Design optimization of an LPG thermal cracker for multiple objectives. Int. J. Chem. React. Eng. 2011, 9, Article A80. (14) Coello Coello, A. C.; Pulido, G. T.; Lechuga, M. S. Handling multiple objectives with particle swarm optimization. IEEE Trans. Evolut. Comput. 2004, 8, 256–279. (15) Kurup, A. S.; Hidajat, K.; Ray, A. K. Optimal operation of an industrial-scale Parex process for the recovery of p-xylene from a mixture of C-8 aromatics. Ind. Eng. Chem. Res. 2005, 44(15), 5703-5714. 27



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 28 of 30

(16) Guria, C.; Verma, M.; Mehrotra, S. P.; etc. Multi-objective optimal synthesis and design of froth flotation circuits for mineral processing, using the jumping gene adaptation of genetic algorithm. Ind. Eng. Chem. Res. 2005, 44(8), 2621-2633. (17) Agrawal, N.; Rangaiah, G. P.; Ray, A. K.; etc. Multi-objective optimization of the operation of an industrial low-density polyethylene tubular reactor using genetic algorithm and its jumping gene adaptations. Ind. Eng. Chem. Res. 2006, 45(9), 3182-3199. (18) Bhutani, N.; Ray, A. K.; Rangaiah, G. P. Modeling, simulation, and multi-objective optimization of an industrial hydrocracking unit. Ind. Eng. Chem. Res. 2006, 45(4), 1354 -1372. (19) Sankararao, B.; Gupta, S. K. Multi-objective optimization of pressure swing adsorbers for air separation. Ind. Eng. Chem. Res. 2007, 46(11), 3751-3765. (20) Amorim, P.; Antunes, C. H.; Almada-Lobo, B. Multi-objective lot-sizing and scheduling dealing with perishability issues. Ind. Eng. Chem. Res. 2011, 50(6), 3371-3381. (21) Shastri Y.; Diwekar U. Stochastic Modeling for Uncertainty Analysis and Multiobjective Optimization of IGCC System with Single-Stage Coal Gasification. Ind. Eng. Chem. Res. 2011, 50, 4879-4892. (22) Mandli, A. R.; Modak, J. M. Evolutionary algorithm for the determination of optimal mode of bioreactor operation. Ind. Eng. Chem. Res. 2012, 51, 1796-1808. (23) Das, S.; Suganthan, P. N. Differential Evolution: A Survey of the State-of-the-Art. IEEE Trans. Evolut. Comput. 2011, 15, 4–31. (24) Babu, B. V.; Angira, R. Modified differential evolution (MDE) for optimization of non-linear chemical processes. Comput. Che. Eng. 2006, 30, 989–1002. (25) Zhang, H. B.; Rangaiah, G. P. An efficient constraint handling method with integrated differential evolution for numerical and engineering optimization. Comput. Chem. Eng. 2012, 37, 74-88. (26) Soliman, O.; Bui, L. T.; Abbass H. A memetic coevolutionary multi-objective differential evolution algorithm. In Multi-Objective Memetic Algorithm; Goh, C. K., Ong, Y. S.; Tan, K. C., Eds.; Springer: Berlin, 2009; pp 369-388. (27) Santana-Quintero, L. V.; Hernández-Díaz, A. G.; Molina, J.; Coello, C. A. C.; Caballero, R. DEMORS: A hybrid multi-objective optimization algorithm using differential evolution and rough set theory for constrained problems. Comput. Oper. Res. 2010, 37, 470–480. 28


Page 29 of 30

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(28) Ali, M.; Siarry, P.; Pant, M. An efficient Differential Evolution based algorithm for solving multi-objective optimization problems. Eur. J. Oper. Res. 2012, 217, 404–416. (29) Zhang, H.; Zhou, J.; Fang, N.; etc. Daily hydrothermal scheduling with economic emission using simulated annealing technique based multi-objective cultural differential evolution approach. Energy. 2013, 50, 24-37. (30) Xu, B.; Qi, R.; Zhong, W.; etc. Optimization of p-xylene oxidation reaction process based on self-adaptive multi-objective differential evolution. Chemometr. Intell. Lab. 2013, 127, 55-62. (31) Sharma, S.; Rangaiah, G. P. An improved multi-objective differential evolution with a termination criterion for optimizing chemical processes. Comput. Chem. Eng. 2013, 56, 155-173. (32) Gao, X. D. Study on the simulation and optimization method for ethylene cracking furnaces. Ph.D. Thesis, The Tsinghua University, June 2008. (33) Qin, A. K.; Huang, V. L.; Suganthan, P. N. Differential evolution algorithm with strategy adaptation for global numerical optimization. IEEE Trans. Evolut. Comput. 2009, 13, 398-417. (34) Deb, K.; Agrawal, R. B. Simulated binary crossover for continuous search space. Complex. Syst. 1995, 9, 115–148. (35) Zhang, Q.; Li, H. MOEA/D: A multiobjective evolutionary algorithm based on decomposition. IEEE Trans. Evolut. Comput. 2007, 11(6), 712–731. (36) Nebro, A. J.; Durillo, J. J.; García-Nieto, J.; Coello, C. A. C.; Luna, F.; Alba, E. SMPSO: A new PSO-based metaheuristic for multi-objective optimization. Proceeding of IEEE Symposium on Computational Intelligence in Multicriteria Decision-Making, Nashville, USA, Mar 30-Apr 2, 2009; IEEE Press: Piscataway, 2009. (37) Zitzler, E.; Deb, K.; Thiele, L. Comparison of multiobjective evolutionary algorithms: empirical results. Evolut. Comput. 2000, 8, 173–195. (38) Kita, H.; Yabunoto, Y.; Mori, N.; etc. Multi-objective optimization by means of the thermodynamical genetic algorithm. Proceeding of 4th Workshop on Parallel Problem Solving from Nature – PPSN IV, Berlin, Germany. Sep 22-26, 1996; Werner, E., Hans-Michael, V., Eds.; Springer: Berlin, 1996. (39) Deb, K.; Thiele, L.; Laumanns, M. etc. Scalable test problems for evolutionary multi-objective optimization. In Evolutionary Multiobjective Optimization: Theoretical Advances and Applications; Ajith, A., Lakhmi, J, Rober, G, Eds.; Springer: Berlin, 2005; 29



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 30 of 30

pp 105–145. (40) Viennet, R.; Fontiex, C.; Marc, I. Multicriteria optimization using a genetic algorithm for determining a Pareto set. Int. J. Syst. Sci. 1996, 27, 255–260. (41) Binh, T. T.; Ulrich, K. Multiobjective evolution strategy with linear and nonlinear constraints. Proceeding of 15th IMACS World Congress on Scientific Computation, Modeling and Applied Mathematics, Berlin, Germany, Aug 24-29, 1997; Sydow, A., Eds.; Wissenschaft & Technik Verlag: Berlin, 1997. (42) Tamaki, H.; Kita, H.; Kobayashi, S. Multi-objective optimization by genetic algorithms: a review. Proceeding of IEEE Conference on Evolutionary Computation, Nagoya, Japan, May 20-22, 1996; IEEE Press: Piscataway, 1996. (43) Li, H.; Zhang, Q. Multiobjective optimization problems with complicated Pareto set, MOEA/D and NSGA-II. IEEE Trans. Evolut. Comput. 2009, 13(2), 284–302. (44) Kumar, P.; Kunzru, D. Kinetics of coke deposition in naphtha. Can. J. Chem. Eng. 1985, 63(4), 598-604. (45) Gijsberts, A. Metta, G.; Rothkrantz, L. Evolutionary optimization of least-squares support vector machines. Ann. Inform. Syst. 2010, 8, 277-297. (46) Yin, P. Y.; Glover, F.; Laguna, M.; Zhu, J. X. Scatter PSO – a more effective form of particle swarm optimization. Proceedings of IEEE Conference on Evolutionary Computation, Singapore, Sep 25-28, 2007; IEEE Press: Piscataway, 2007.

30


Multiobjective Operation Optimization of Naphtha Pyrolysis Process

Recommend Documents