A protein inspired RNA genetic algorithm for parameter estimation in hydrocracking of heavy oil

Chemical Engineering Journal 167 (2011) 228–239

Contents lists available at ScienceDirect

Chemical Engineering Journal journal homepage: www.elsevier.com/locate/cej

A protein inspired RNA genetic algorithm for parameter estimation in hydrocracking of heavy oil Kangtai Wang, Ning Wang ∗ National Laboratory of Industrial Control Technology, Institute of Cyber-Systems and Control, Zhejiang University, Hangzhou 310027, China

a r t i c l e

i n f o

Article history: Received 21 August 2010 Received in revised form 26 November 2010 Accepted 13 December 2010 Keywords: Protein inspired RNA genetic algorithm (PIRGA) Biological computing Parameter estimation Hydrocracking of heavy oil Kinetic modeling

a b s t r a c t Hydrocracking is a crucial process in refineries and suitable model is useful to understand and design hydrocracking processes. Simulating the procedure from RNA to protein, a protein inspired RNA genetic algorithm (PIRGA) is proposed to estimate the parameters of hydrocracking of heavy oil. In the PIRGA, each individual is represented by a RNA strand and a new fitness function combining traditional fitness value and individual ranking is employed to maintain population diversity. Furthermore conventional crossover operators are replaced by RNA-recoding operator and protein-folding operators to improve the searching ability. An adaptive mutation probability in the PIRGA makes the algorithm have more chance to jump out of local optima. Numerical experiments on seven benchmark functions indicate that the PIRGA outperforms other genetic algorithms on both convergence speed and accuracy greatly. 10 parameters are obtained by the PIRGA and the kinetic model for hydrocracking of heavy oil is established. Experimental results reveal that the predictive values are in good agreement with the experimental data with relative error less than 5%. The effectiveness and the robustness of the model are also validated by experiments. © 2010 Elsevier B.V. All rights reserved.

1. Introduction The conflict of growing demands of middle distillates with overabundance of heavy crude oils makes people focus on hydrocracking technique and hydrocracking is increasingly becoming a crucial secondary petroleum refinery processes to treat heavy oil, such as vacuum resid. The main objective of this process is to convert heavy molecules like vacuum gas oil (VGO) into lighter and more valuable fractions such as naphtha, gasoline, and hydrogen gas. Different approaches have been utilized for kinetic modeling of hydrocracking, varying from the most common and used lumping technique to more complex models based on continuous mixture or single events. A comprehensive review of hydrocracking modeling can be found in [1]. Accurate analytical or numerical modeling of these new upgrading processes is essential, in order to correctly interpret experiment measurements and to lead to a better understanding and design of industrial-scale processes. However, many undesirable and complicated reactions in hydrocracking process bring an arduous task to make a reasonable tradeoff between accuracy and complexity of modeling. Usually, people consider these reactions dominant and significant in all reaction and neglect tiny ones. After getting

∗ Corresponding author. E-mail address: [email protected] (N. Wang). 1385-8947/$ – see front matter © 2010 Elsevier B.V. All rights reserved. doi:10.1016/j.cej.2010.12.036

mechanism model of hydrocracking processes, determination of model parameters is also difficult by conventional deterministic optimization methods [2]. Most deterministic optimization methods encounter some deficiencies such as sensitive to initial value and/or requiring of differentiable information of optimized problems. Genetic algorithms (GAs) first developed by Holland are stochastic search technique based on the mechanism of natural selection and survival of the fittest [3,4]. Those individuals with higher fitness value are assigned higher survival probability. Genetic algorithms do not tackle problems directly but in a chromosome space, i.e., candidate solutions of an optimization problem are converted to chromosome firstly. Due to no requirement of prior information about search space and owning excellent global search ability, GAs have been applied widely to address complicated real-world problems with non-differentiable, non-convex and non-linear [5–10]. Despite successful applications of GAs to optimization problems, traditional genetic algorithms suffer from slow convergence speed and poor local searching ability. Furthermore, simple operators of GAs limit their searching efficiency in searching space and redundant search dominates the whole procedure. In order to overcome these drawbacks of GAs, some improved genetic algorithms simulating biological behavior emerged. Inspired by DNA molecular structure and operators, Tao et al. proposed a RNA genetic algorithm, which adopts quaternary encoding and three RNA molecular operators to improve the searching capability of GAs [11]. Wang et al. presented a novel RNA genetic algorithm to

K. Wang, N. Wang / Chemical Engineering Journal 167 (2011) 228–239

Nomenclature N population size M individuals number in mating pool l length of RNA individuals Fi original fitness value of individual i new fitness value of individual i Fi Ri ranking number of individual i Nbest population size of neutral group number of left-shift steps n1 n2 number of right-shift steps PRNA-recoding probability of RNA-recoding operator Pmutual-folding probability of protein mutual-folding operator Pself-folding probability of protein self-folding operator mutation operator probability Pm Gmax number of maximal generations D dimensions of benchmark functions g generation number VGO vacuum gas oil reaction rate of product composition i ri ωi mass fraction of product composition i E relative error T0 starting time of collecting data stopping time collecting data Tt yi (nt) experimental data of composition i at time of nt yˆ i (nt) predictive value of composition i at time of nt

make a reasonable tradeoff between searching accuracy and computational efforts by introducing a new RNA molecular operator, stem-loop operator, and two different fitness values. Computational results showed that this RNA genetic algorithm solved three parameter estimation problems of dynamic systems successfully [12]. Chen et al. developed another DNA genetic algorithm to address parameter estimation problem in hydrogenation reaction well [13]. Tao et al. and Chen et al. suggested two different hybrid genetic algorithms incorporating DAN double helix genetic algorithm and sequential quadratic programming (SQP) for gasoline blending scheduling problems [14,15]. Simulating RNA molecular operators and procedure from DNA to protein in biological cell, a protein inspired RNA genetic algorithm (PIRGA) is proposed in this work. In this algorithm, we first encode each individual with a strand of nucleotide bases, RNA strand. Then RNA-recoding operator and protein-folding operators are designed to replace conventional crossover operators to improve the performances of GAs. Apart from encoding procedure, RNA strands are first translated into amino acids ones, protein strands, according to triplet codons in decoding procedure. Numerical solutions on seven benchmark functions, varying from two-dimensional to ten-dimensional, show the superiority of the PIRGA in contrast to other genetic algorithms. The parameters estimation problem in hydrocracking kinetic model is also well addressed by the proposed genetic algorithm. 2. Protein inspired RNA genetic algorithm (PIRGA) 2.1. Representing and decoding DNA is the major genetic material for life and contains plentiful genetic information. DNA molecular owns a double helix structure and during DNA replication and transcription, the double helix must be separated transiently and reversibly. Because the separated antisense strand, RNA individual, contains almost all useful information of DNA and is simpler on structure than DNA molecular, RNA is also regarded as main genetic material of living cells. RNA

229

Table 1 Relationship between triplet codons and 64-ary integers. First nucleotide

U U U U C C C C A A A A G G G G

Second nucleotide

Third nucleotide

U

C

A

G

0 1 2 3 16 17 18 19 32 33 34 35 48 49 50 51

4 5 6 7 20 21 22 23 36 37 38 39 52 53 54 55

8 9 10 11 24 25 26 27 40 41 42 43 56 57 58 59

12 13 14 15 28 29 30 31 44 45 46 47 60 61 62 63

U C A G U C A G U C A G U C A G

contains 4 kinds of nucleotides: Adenine (A), Uracil (U), Guanine (G) and Cytosine (C). In the PIRGA, all individuals are represented by



l

RNA strands and the encoding space is E = A, U, G, C , where l is the length of RNA strands (All strands are the same length in the PIRGA). In genetic code, a three-letter codes, triplet codon, decides a amino acid, i.e., three nucleotides in RNA strand decide a amino acid by the translating operator. On the decoding procedure of the PIRGA, RNA strands are first translated into amino acid strands, protein ones. In biological field, 20 kinds of common amino acids are recognized by codons, and some different triplet codons decide the same amino acid. Though the degeneracy of triplet encoding in cell plays an important role in reducing replication error and abnormal mutation, in GAs we hope individuals and candidate solutions are one–one relationship which can reduce computational costs in searching procedure. Hence on decoding procedure, RNA strands are first translated into 64-ary. The relationship between triplet codons and 64-ary integers is shown in Table 1. After translating operator individuals are represented by amino acid ones, i.e., integer sequences between [0,63] which are easy to convert to real-values in domain of an optimization problem. This indirect decoding method first translates RNA individuals into amino acid sequences and then these 64-ary amino acid strands are converted to real-values. The length of RNA individuals is decided by variable number of problem and precision. For example, there are n variables and each variable is represented by m 64-ary integer. The length of RNA individuals is l = 3 × n × m. The decoding procedure is shown in Fig. 1. Encoding procedure is a mapping from solution space to RNA encoding space, i.e., RNA strands are adopted to represent individuals in the PIRGA. Subsequently, these RNA individuals change their forms in the encoding space by some RNA and protein molecular operators. When the algorithm terminated, Optimization results are obtained by decoding procedure, which is an inverse mapping from RNA space to problem space. In other words, the encoding and the decoding procedures complete a mutual mapping between solution space and RNA operator space.

Fig. 1. Decoding procedure.

230

K. Wang, N. Wang / Chemical Engineering Journal 167 (2011) 228–239

Fig. 2. RNA-recoding operator.

2.2. Selection strategy Selection operator decides which individuals are kept to next generation and how many individuals are replicated to offspring. Generally, those individuals with higher fitness value have a greater chance to survival or reproduce. Many types of selection methods were put forward, such as proportionate selection, Boltzmann selection, rank selection, tournament selection etc. The proportionate selection and tournament selection are the most popular and useful selection methods. Individual numbers reproduced by the proportionate selection method are proportional to their fitness values. When the fitness values of several individuals are larger than the others greatly, those individuals will domain all offspring individuals, which lead premature convergence easily. Tournament method chooses several individuals each time and the best one is selected as one of offspring individuals. Therefore the tournament method only relates to ranking of fitness values in current population and cannot reflect discrepancy of fitness values. We hope to maintain more useful individuals to next generation and it is necessary to design a suitable selection strategy. Combining the advantages of above two selection methods, a new fitness function is constructed in the PIRGA. The details of this fitness function and the selection method are described as follows. To fitness value Fi and ranking number Ri , we define the new fitness value Fi as: Fi i=M F i=1 i

Fi = ˛1 

Ri i=M R i=1 i

+ ˛2 

(1)

where M is the number of candidate individuals, ˛1 and ˛2 are weights of the original fitness value and the ranking number respectively. This fitness value considers the original fitness values and the ranking number at the same time and the new fitness value does not change individual original ranking. After calculating the novel fitness values of all candidate individuals, tournament selection method is adopted according to the new fitness values. In order to simplify parameter setting, we set ˛1 = ˛2 = 1 in this work. In the PIRGA, all individuals in mating pool are first divided into two classes according to their original fitness values. The best Nbest individuals in terms of the original fitness values are neutral group and the left ones are bad group. Since neutral individuals have higher fitness values than the others, keeping those individuals to

next generation is hoped to generate better solutions. On the contrary, introducing those individuals in bad group to next generation cannot improve final solutions generally. However, maintaining bad individuals can increase population diversity, which is also an important aspect to avoid premature convergence. In order to address the confliction of maintaining more neutral individuals and increasing population diversity, we change ranking numbers of some individuals when they are selected. If one reproduced individual i belongs to neutral group, the ranking number of that is unchanged. However, if one reproduced individual i belongs to bad group, we hope to give more chances to other individuals for increasing the diversity of population. So if one bad group individual is reproduced, we decrease ranking number of this individual, Rinew = Ri − 1, and the new fitness value of this individual is also reduction. By this selection strategy, we make a tradeoff between keeping more neutral group individuals and maintaining population diversity. 2.3. Crossover operators Crossover operator in GAs simulates a process of sexual reproduction that offspring individuals inherit parts of parents respectively and crossover operator is the most important operator in GAs. However, conventional crossover operators, such as singlepoint, multi-points and arithmetic crossover operators, are all simple abstraction on biological concepts, which weakens searching performance of GAs. In order to decrease useless searching, exploit more useful information of current individuals and make algorithm have a better ability to explore searching space, we design some novel molecular operators to replace these conventional crossover operators. Inspired by procedure translating from RNA strand to protein, RNA-recoding and protein-folding operators are adopted in the proposed PIRGA. The details of these operators are depicted as follows. 2.3.1. RNA-recoding operator In biological field, before a RNA strand translated into a protein one, the RNA strand often change its form for correcting errors of reproduction, which is called RNA-recoding [16]. Simulating this biological phenomenon, we design the RNA-recoding operator as shown in Fig. 2.

K. Wang, N. Wang / Chemical Engineering Journal 167 (2011) 228–239

1

2

3

4

5

6

7

8

9

10

11

12

231

13

14

15

Original strand 1

U C A A GC U C U A CGGGC A U GU A GC A AGA C A C C U U U A UGGA U C C A GC A Original strand 2

A GCGA C C GU GC U A A GC U U GGA C A C U GUGC C A UGU C U GGGC A U C GC

U C A AGC U C U A CGGGC A A A GC U U C U G G A GC GA C CGU

U C A AGC U C U A CGA A GA GGC C U U C U G G A GC GA C CGU

A C G U U U A U GGA U C C A C GGGC C U A A C AG C U A U C GC A U GU A GC G U A C C GA C A C U GUG

A C G U U U A U GU GGC C A C A A U C C U A A C AG C U A U C GC A U GU A GC G U A C C GA C A C U GUG

Final strand 1

Final strand 2

Fig. 3. Protein mutual-folding operator.

Current RNA individual splits into two subsequences at a random location. Left subsequence makes left-shift n1 steps and right one makes right-shift n2 steps. Subsequently, corresponding nucleotides of another random RNA sequence are inserted to the blank between two subsequences. In order to keep the same length, superfluous nucleotides of current RNA are cut and discarded in RNA-recoding operator. 2.3.2. Protein-folding operators After RNA strands translated into protein ones, some other operators such as modification operator and protein-folding operator are needed to form living proteins in biological field [16]. In the PIRGA, we design protein-folding operators to improve searching ability. Protein-folding operators can be classified into two groups, the mutual-folding operator and the self-folding operator. The mutual-folding operator takes place between two different individuals. However the self-folding operator occurs in single one. Mutual-folding operator: two amino acid strands, protein ones, fold each other and some amino acids of corresponding location connect and influence mutually. That is to say that some amino acids belonging to different individuals communicate their information by this mutual-folding operator. To two nucleotides influenced, if they are the same, each nucleotide is replaced by another different randomly. Otherwise, we exchange the two nucleotides. The detail of mutual-folding operator is shown in Fig. 3.

From Fig. 3, we can see that the fifth and the thirteenth amino acids are influenced by another strand and we take the thirteenth amino acid for example. Original nucleotides of the thirteenth amino acid strands are GAU and GGG respectively. Because the two nucleotides of the first location are the same, G, we use two different nucleotides generated randomly, U and A to replace them. To the second and the third nucleotides, we exchange them directly because they are different. Self-folding operator: a protein strand folds and some nucleotides influence each other. To two nucleotides influenced, if they are the same, each nucleotide is replaced by another different randomly. Otherwise, we exchange the two nucleotides. The detail of the selffolding operator is shown in Fig. 4. Because stochastic searching algorithm cannot insure the best individual of next generation is better than current one, elitism is applied to avoid degradation in the PIRGA, i.e., the best individual of current generation is reproduced to next generation directly at the end of selection operator. On the procedure of crossover operator, new individuals generated are added into mating pool directly and do not replace father individuals in mating pool. Hence, after crossover operator, the number of individuals in mating pool increases. Because RNA-recoding operator plays an important role in crossover procedure, the probability of RNA-recoding operator is relatively larger than that of protein-folding operators. Through many experiments, we draw the conclusion that setting the probability of RNA-recoding operator as [0.6,0.9] is suitable in the PIRGA.

232

K. Wang, N. Wang / Chemical Engineering Journal 167 (2011) 228–239

U C A AGC U C U A C GGGC A U GU A GC A A GA C A C C U U U A UGGA U C C A GC A

Folding

C A U GU A G G C G A G A C A GA C C A C U U U A U GGA U C C A GC A U C A A GC U C U

C A U GU A G G C G A A GC A GA U G C A U U U A U GGA U C C A GC A U C A A GC U C U Fig. 4. Protein self-folding operator.

To decrease computational efforts, protein-folding operators are only taken place among those father individuals that RNA-recoding operator is not executed. For example, if we set the probability of RNA-recoding operator as 0.8 (PRNA-recoding = 0.8), the probability of protein-folding operator is 0.2 (Pmutual-folding + Pself-coding = 0.2). As discussed above, Pmutual-folding and Pmutual-folding are relatively small and we set Pmutual-folding = Pmutual-folding in the PIRGA to decrease the difficulty of parameter setting. Therefore, the probabilities of mutual-folding operator and self-folding operator are equal to 0.1 when the probability of the RNA-folding operator is 0.8. To each individual in mating pool, generates a random number between 0 and 1. If the random number is smaller than PRNA-recoding , RNA-recoding operator is executed in this individual. Otherwise, regenerates a random number. If the number regenerated is smaller than 0.5 (because the probability of mutual-folding operator is equal to that of self-folding operator), mutual-folding operator is carried out. If that is greater than or equal to 0.5, self-folding operator is executed. In a word, every individual in mating pool must be executed one of the three crossover operators, RNA-recoding operator, mutualfolding operator and self-folding operator, one time. Therefore, if population size is N, there are 2N individuals in mating pool after crossover operators.

current optima. The adaptive mutation probability in the PIRGA is shown in Eq. (2) Pm = a1 +

a2 1 + e(a3 −g)/a4

(2)

where a1 , a2 , a3 , a4 are constants, g is the number of current generation. Constant a1 decides a minimum of mutation probability at the beginning stage, a2 denotes the variation range of mutation probability, a3 reveals the place where most change of mutation takes place, a4 denotes the slope degree of the location where most change mutation probability occurs. In this work, these parameters are set as: a1 = 0.1, a2 = 0.1, a3 = 0.5 × Gmax , a4 = 0.1 × Gmax . Gmax is the maximal iteration number. The curve of the adaptive mutation probability is shown in Fig. 5 when Gmax is equal to 1000. Mutation operator takes place among all individuals in mating pool with adaptive mutation probability and new individuals generated by mutation operator replace old ones. Therefore, mutation operator does not change individual number in mating pool.

0.2 0.19

2.4. Mutation operator Mutation operator plays an important role when algorithms trap into local optima. The mutation operator can guide algorithms to jump out of local optimum when most individuals in population become similar. The mutation operator in the PIRGA is the same as binary GA, i.e., using a random different code replaces original one. Mutation probability has a significant influence on quality of final results and setting a suitable value of mutation probability is also an arduous task. In the PIRGA, an adaptive mutation probability is adopted to improve performance of global exploration. At the beginning stage of evolutions, a small value of mutation probability is enough because of a good distribution or diversity of population. However, a comparatively large one is necessary at the end stage of running in order to get more useful information and jump out of

Mutation probability

0.18 0.17 0.16 0.15 0.14 0.13 0.12 0.11 0.1

0

100

200

300

400

500

600

700

Generation Fig. 5. Mutation probability curve.

800

900

1000

K. Wang, N. Wang / Chemical Engineering Journal 167 (2011) 228–239

233

Fig. 6. Procedure of the PIRGA.

Table 2 Formula of benchmark functions. Test functions f1 =

2

2 ) + (1 − x1 ) − x 2

100(x12

x2 +x2

sin

f2 = 0.5 + f3 = (x12 +

−0.5

2

2

1.0+0.001(x2 +x2 ) 1 2 0.25 2 x22 ) [sin (50(x12

D 

f4 = −

1

2



xi sin(

+

0.1 x22 ) )

+ 0.1]

|xi |)

Bounds

Optimal value

−2.048 ≤ xi ≤ 2.048

0

−100 ≤ xi ≤ 100

0

−100 ≤ xi ≤ 100

0

−500 ≤ xi ≤ 500

−418.9829 × D

−5.0 ≤ xi ≤ 5.0

0

−100 ≤ xi ≤ 100

0

−50 ≤ xi ≤ 50

0

i

f5 =

D 

[xi2 − 10 cos (2xi ) + 10]

 ⎞2  D  D   1 1

⎝ ⎠ −0.2 xi − exp cos(2xi ) + 20 + e f6 = −20 exp N N ⎛

i

i=1

f7 =

1 4000

D 

i=1

xi2 −

 

D 

cos i=1

i=1

x

√i

i

+1

234

K. Wang, N. Wang / Chemical Engineering Journal 167 (2011) 228–239

2.5. Procedure of the PIRGA To solve complex optimization problems with the PIRGA, we first initialize a population of RNA individuals. Then, crossover and mutation operators are employed to improve their fitness values. After each generation, some candidate individuals are reproduced to combine next population. Based on the above encoding method, selection strategy and genetic operators, the procedure of the PIRGA can be summarized as follows. Seen in Fig. 6.

Step1: Initialize a population with N individuals. Step2: Calculate fitness value of each individual and divide individuals into two groups. Step3: Select N individuals to combine mating pool with the selection strategy. Step4: Judge whether or not meets the condition of the RNArecoding operator. If affirmative, execute the RNA-recoding operator and jump to Step 6. Otherwise, jump to Step 5. Step5: Judge whether or not meets the condition of the protein mutual-folding operator. If affirmative, execute the protein mutual-folding operator and jump to Step 6. Otherwise, carry out protein self-folding operator and jump to Step 6. Step6: execute mutation operator with adaptive mutation probability. Step7: Repeat Steps 2–6 until termination conditions are met and the final solution is found.

3. Computational experiments To investigate the performance of the protein inspired RNA genetic algorithm, a test environment is constructed in terms of seven typical benchmark functions. Many different kinds of optimization problems are covered by these benchmark functions. The details of the benchmark functions are shown in Table 2. Functions f1 − f3 are 2-dimentional and f4 − f7 are 10 dimensional. In order to compare results with those in literatures published, we do 30 independent runs, the same times as literatures [17,18], to each function. In the experiments, the parameters of the PIRGA are set as: population size N = 80, the probability of RNA-recoding operator PRNA-recoding = 0.8, the number of best group Nbest = 10, the number of maximal generation Gmax = 5000. The seven benchmark functions own some typical characteristics for real-world optimization problems. f1 is the well-known Rosenbrock function with saddle shape. This function has a single optimum at the point x* = (1, 1)T and finding the optimum is not very difficult even using a simple local search method. Therefore, the function f1 is often used to test local searching ability. The global optimum of f2 locates at the center of search space and many local optima are around the global optimum symmetrically. That is to say, algorithms are likely to fall into a local optimum in all directions. Hence, obtaining the global optimum is a tremendous task for this function. Function f3 owns many local optima around global optimum and locating the global optimum with sufficient precision is not so easy. The Schwefel function f4 is a multimodal problem, and its global minimum is near the bound of domain and geometrically far from the second minimum points. Therefore, algorithms are potentially prone to converge to in a wrong direction. Function f5 is a famous test function, Rastrigin function, for validating anti-deceptive searching ability of algorithms. Inseparable and multimodal characteristics of function f6 make it more challenging to find the global optimum. Griewank function f7 is also a widely used multimodal problem with thousands of regularly distributed minima.

Table 3 Experimental results of f1 − f3 . Method

Function

Best

Average

SD

SSGA [17]

f1 f2 f3

1.57592e−10 6.94578e−12 1.86689e−3

2.5781e−3 1.7474e−3 1.5436e−2

7.6389e−3 3.6466e−3 9.8162e−3

PfGA [17]

f1 f2 f3

3.77082e−11 2.77778e−13 8.18957e−4

7.6235e−4 6.5457e−3 3.6606e−2

1.7286e−3 5.2392e−3 4.0253e−2

IGA [17]

f1 f2 f3

6.30557e−12 2.77556e−13 7.30713e−4

1.2780e−4 9.4619e−3 7.3071e−4

2.2475e−4 1.8021e−2 0

PIRGA

f1 f2 f3

4.0107e−21 0 1.8461e−8

3.4376e−8 0 1.1e−3

1.7935e−7 0 1.2e−3

To a random searching algorithm, average result of many times is required to validate the robustness. Therefore, all experiments on benchmark functions are done for 30 runs to get the best solution, average solution and standard deviation. The detailed results of 2-dimensional functions f1 − f3 are listed in Table 3 and those of 10-dimensional functions f4 − f7 are listed in Table 4. Some results reported in literatures [17–20] are also listed in Tables 3 and 4 for comparing. From Table 3, the best solution of 30 runs on f1 − f3 is better than SSGA [17], PfGA [17] and IGA [17] greatly. Especially to function f2 , the proposed algorithm finds the true global optimum every time, which indicates the perfect capabilities of the PIRGA on global exploration and local exploitation. The average solution and the standard deviation are crucial indexes to verify robustness of stochastic algorithms. From Table 3, it can be concluded that the performances of the PIRGA on searching accuracy and robustness about the three 2-dimensional functions are better than the other genetic algorithms. To high-dimensional function f4 , the PIRGA gets the optimum of −4189.7 and average optimum of −4185.8, which are also better than that of the other algorithms. The standard deviation of the PIRGA on f4 is bad than that of flc-aGA [19]. However, the superiorities on ‘best solution’ and ‘average solution’ show that even the worst solution of the PIRGA is likely better than the best one of the other algorithms. This conclusion can also be applied to functions f5 − f7 . Final population distribution is relative to the performance of searching algorithm because bad distribution of final population leads a premature convergence easily. Therefore, some experiments about the distributions of final population are done. The details of initial and final populations of the PIRGA on f1 − f3 are shown in Fig. 7. Table 4 Experimental results of f4 − f7 . Method

Function

Best

Average

flc-aGA [19]

f4 f5 f6 f7

−4166.058 3.181 11.88 0.041

−4151.103 5.366 13.539 0.090

SD 0.0045 0.0537 0.1353 0.0009

stGA [20]

f4 f5 f6 f7

– – – –

−4178.400 3.020 12.970 0.087

– 0.0302 0.1297 0.0009

atGA [18]

f4 f5 f6 f7

−4181.433 1.520 7.266 0.038

−4157.679 2.855 9.492 0.076

– 0.0285 0.0949 0.0008

PIRGA

f4 f5 f6 f7

−4189.7 5.93668e−7 6.34952e−4 6.28619e−8

−4185.8 0.967286 0.091424 0.069948

13.6558 1.186335 0.21590 0.038790

K. Wang, N. Wang / Chemical Engineering Journal 167 (2011) 228–239

1.5

1

1

0.5

0.5

x2

1.5

x2

0

0

-0.5

-0.5

-1

-1

-1.5

-1.5 -2

-2 -2

-1.5

-1

-0.5

0

0.5

1

-2

1.5

-1.5

-1

-0.5

x1

1.5

80

60

60

40

40

20

20

x2

x2

1

100

80

0

0

-20

-20

-40

-40

-60

-60

-80

-80

-80

-60

-40

-20

0

20

40

60

80

100

-100 -100

-80

-60

-40

-20

0

20

40

60

x1

x1

c) Initial population of f 2

d) Final population of f 2

100

80

100

100

80

80

60

60 40

20

20

x2

40

0

0

-20

-20

-40

-40

-60

-60

-80 -100 -100

0.5

b)Final population of f1

100

-100 -100

0

x1

a)Initial population of f1

x2

235

-80

-80

-60

-40

-20

0

20

40

60

80

100

-100 -100 -80

-60

-40

-20

0

20

40

60

80

100

x1

x1

e) Initial population of f 3

f) Final population of f3

Fig. 7. Population distribution of the PIRGA on f1 − f3 .

Because initial individuals are generated randomly, the individuals are uniformly distributed over the whole searching space. With the running of the algorithm, more and more individuals locate to a narrow space near the global optimum. As shown in (b), (d) and (f) of Fig. 7, some individuals find the optimum and most of the others distribute around the optimum. At the same time a few individuals locate far from global optimum, which is very useful to jump out of local optima and avoid premature convergence.

Convergence speed is also an important index of optimization algorithms. The convergence curves of the proposed algorithm on function f4 and f5 are shown in Fig. 8. From the figure, we can learn that the proposed algorithm can converge quickly to near optimum less than 1000 generations. Especially to function f4 , the global optimum is obtained less than 500 generations, which is perfect to a 10 dimensional multimodal function.

236

K. Wang, N. Wang / Chemical Engineering Journal 167 (2011) 228–239 120

-1000 Best Average

-2000 -2500 -3000 -3500

80 60 40 20

-4000 -4500 0

Best Average

100

Function value

Function value

-1500

1000

2000

3000

4000

5000

0

0

1000

Generation

2000

3000

4000

5000

Generation

b) Convergence curves of f 5

a) Convergence curves of f 4

Fig. 8. Convergence curves of f4 and f5 .

4. Parameters estimation problem of hydrocracking of heavy oil Reactions of hydrocracking of heavy oil are cumbersome and getting a suitable model is a huge task for engineers and researchers. Many models of hydrocracking processes are established using various approaches [1,21]. After analyzing detailed mechanisms of different reactions, experimental data based parameter estimation becomes a crucial step for a suitable model of hydrocracking of heavy oil. Moreover complex compositions and reactions, as well as numerous parameters, increase the difficulty of modeling and parameter estimation of hydrocracking. The kinetic model used in this work is proposed by Sanchez et al. [21] and researched by Hassanzadeh and Abedi [22]. The temperature of the reactor is maintained at the desired level by using a three-zone electric furnace. A Ni/Mo commercial catalyst was employed for all experiments. A total of 100 ml of catalyst was loaded into the reactor and in situ activated by sulfiding with a hydrodesulfurized naphtha containing 0.8 wt% CS2 at the following operating conditions: pressure of 54 kg/cm2 , H2 /oil reaction of 2000 ft3 /bbl, and LHSV of 3.2/h. The model includes five lumped pseudo-components: unconverted residue, vacuum gas oil (VGO), distillates, naphtha, and gases. For a very high excess of hydrogen, all reactions are considered as the pseudo first order. The relationship of five pseudo-components is shown in Fig. 9. The reaction rates are given by Eqs. (3)–(7) [21]. rR = −(k1 + k2 + k3 + k4 )ωR

(3)

rVGO = k1 ωR − (k5 + k6 + k7 )ωVGO

(4)

Residue : VGO :

Distillates : Naphtha : Gases :

rD = k2 ωR + k5 ωVGO − (k8 + k9 )ωD

(5)

rN = k3 ωR + k6 ωVGO + k8 ωD − k10 ωN

(6)

rG = k4 ωR + k7 ωVGO + k9 ωD + k10 ωN

where ri = dωi /dt is reaction rate of each fraction. All the experimental data are generated according to the parameters of the literature [22] under the conditions of 3.45 Mpa and 320, 350, 380 ◦ C. The feed composition and kinetic parameters else, such as pre-exponential factor and activation energy, can also be found in the literature [22]. The feed API gravity for the three experiments are 10.996, 10.875 and 11.311 respectively. All experiments of parameter estimation are performed with residence times of 3–72 h and all data are collected every 5 h. This model contains 10 unknown parameters and five first order differential equations. Because reaction rate is low enough, defining the searching scope of each parameter as [0,0.1] is wide enough. The objective function is to minimize the sum of square errors between the experimental data and calculated ones of the model estimated. To simulate real situation, all experiment data are added normal distribution noise with mean of zero and variance of 0.01. All experiments run 30 times in MATLAB with a double-core computer, 1.8 Hz main frequency. Table 5 summarizes all obtained parameter values of the five reactions under the conditions of 320, 350, 380 ◦ C. Then we calculate the product composition by solving Eqs. (3)–(7) according to these parameter values given in Table 5. The calculated product compositions are compared with the experimental data, and the results are shown in Fig. 10. From Fig. 10, we can see clearly that almost all experimental data distribute near the prediction curves. Average optima of objective function at the three different temperatures are 7.505 × 10−3 , 8.514 × 10−3 and 7.304 × 10−3 respectively. The average optima of objective function indicate that the errors between the experimental data and the predictive data are very small and the predictive data fit experimental data well. Standard deviations of the optima at the three temperatures are 9.369644 × 10−4 at 320 ◦ C, 5.169989 × 10−4 at 350 ◦ C and 2.398 × 10−3 at 380 ◦ C, which demonstrates an excellent robustness of the proposed genetic algo-

(7) Table 5 Kinetic parameters estimated by the PIRGA. Kinetic parameters

Fig. 9. Reaction kinetic model for the upgrading process of bitumen [21].

k1 k2 k3 k4 k5 k6 k7 k8 k9 k10

Temperature 320 ◦ C

350 ◦ C

380 ◦ C

0.001515 0.000268 0.000303 0.000200 0.000275 0.000027 0.000298 0.001953 0.000002 0.001131

0.005070 0.006437 0.000289 0.000476 0.000005 0.001673 0.000048 0.008160 0.000263 0.005005

0.027763 0.034379 0.025036 0.005320 0.007916 0.004741 0.001956 0.004589 0.001449 0.001533

K. Wang, N. Wang / Chemical Engineering Journal 167 (2011) 228–239 0.8

0.8 Experimental data Residue VGO Distillates Naphtha Gases

0.7 0.6

0.4 0.3

0.6 0.5 0.4 0.3

0.2

0.2

0.1

0.1

0

0 0

10

20

30

40

50

60

-0.1

70

Time

75

80

85

90

95

100

Time

a)Predicted values comparing with experimental data at 320ºC

a)Predicted values comparing with experimental data at 320ºC

0.8

0.8

Experimental data Residue VGO Distillates Naphtha Gases

0.7 0.6

0.4 0.3

Experimental data Residue VGO Distillates Naphtha Gases

0.7 0.6 0.5

Mass frction

0.5

Mass frction

Experimental data Residue VGO Distillates Naphtha Gases

0.7

Mass frction

Mass frction

0.5

-0.1

237

0.2

0.4 0.3 0.2

0.1

0.1

0

0 -0.1

0

10

20

30

40

50

60

70

-0.1

Time

75

80

85

90

95

100

Time

b)Predicted values comparing with experimental data at 350 ºC

b)Predicted values comparing with experimental data at 350ºC 0.8 Experimental data Residue VGO Distillates Naphtha Gases

0.8 0.7 0.6

Mass frction

0.5 0.4 0.3

0.7 0.6 0.5

Mass frction

Experimental data Residue VGO Distillates Naphtha Gases

0.4 0.3 0.2

0.2

0.1 0.1

0 0

-0.1 -0.1

0

10

20

30

40

50

60

70

Time

c)Predicted values comparing with experimental data at 380ºC Fig. 10. Predicted product composition comparing with experimental data at the times of 3–72 h. (a) Predicted values comparing with experimental data at 320 ◦ C. (b) Predicted values comparing with experimental data at 350 ◦ C. (c) Predicted values comparing with experimental data at 380 ◦ C.

75

80

85

90

95

100

Time

c)Predicted values comparing with experimental data at 380ºC Fig. 11. Predictive values comparing with experimental data at the times of 72–102 h. (a) Predicted values comparing with experimental data at 320 ◦ C. (b) Predicted values comparing with experimental data at 350 ◦ C. (c) Predicted values comparing with experimental data at 350 ◦ C.

238

K. Wang, N. Wang / Chemical Engineering Journal 167 (2011) 228–239

0.8 Experimental data Residue VGO Distillates Naphtha Gases

0.7 0.6

Mass frction

0.5 0.4 0.3 0.2 0.1 0 -0.1

0

10

20

30

40

50

60

70

Time

a)Predicted values comparing with experimental data at 320ºC 0.8 Experimental data Residue VGO Distillates Naphtha Gases

0.7 0.6

Mass frction

0.5 0.4

  y (nt) − yˆ i (nt)    yi (nt)

5 Tt

E=

i=1

i nt=T0 5 Tt nt=T0 i=1



(8)

where yi (nt) is the experimental data of composition i at the time of nt, yˆ i (nt) is the predictive value of composition i at the time of nt. T0 and Tt are starting time and stopping time of experiments. E reflects a mismatch degree of predictive model and real model. If E is large, we confirm that at least some predictive values are far away experimental data and the model we established is not fit for the real one. In experiments, we calculate the E at three different temperatures and the results are 3.8%, 4.0% and 3.3% respectively. Considering the experimental data are added normal distribution noise with variance of 0.01, the precision of the estimated model is acceptable. In order to verify the effectiveness of the model obtained, some experimental data with times of 71–102 h are gathered to compare with predictive data of the model and the results are shown in Fig. 11. Fig. 11 reveals that though those experimental data are not used to estimate the parameter of the model, they make a reasonable match to predictive values. To a practical model, robustness is common important to tiny disturbance. Hence, we add some disturbance, normal distribution noise with mean of zero and variance of 0.01, to feedstock to investigate the robustness of the model. Fig. 12 is the result that the predictive data of the model with noise feedstock compare with the experimental data at the times of 3–72 h. The curves reveal a suitable match when feed compositions change.

0.3

5. Conclusions

0.2 0.1 0 -0.1

0

10

20

30

40

50

60

70

Time

b)Predicted values comparing with experimental data at 350ºC 0.8 Experimental data Residue VGO Distillates Naphtha Gases

0.7 0.6 0.5

Mass frction

rithm. In order to validate the model more clearly, we define a relative error E as follows.

0.4

Kinetic modeling is an important issue for optimization of chemical processes and complex reactions make the parameter estimation problem in modeling become an arduous task. In this work, a protein inspired genetic algorithm (PIRGA) is proposed for global optimization problems as well as parameter estimation of hydrocracking of heavy oil. The algorithm takes RNA based encoding together with some simulating molecular operators for improving the searching performance. A new fitness function combining original fitness value and individual ranking number is established for making a tradeoff between increasing population diversity and maintaining more useful individuals. Experimental results on test functions demonstrate the power of the proposed algorithm in terms of optimization quality, convergence speed and the distribution of final population. The model we obtained by PIRGA is fit for experimental data well with relative error less than 5% and the effectiveness and the robustness of the model are also validated by some experiments.

0.3

Acknowledgements

0.2

This paper has been supported by the National Natural Science Foundations of China under grants Nos. 60874072 and 60721062 and the National High Technology Research & Development Program of China under grants No. 2008AA04Z209.

0.1 0 -0.1

0

10

20

30

40

50

60

70

Time

c)Predicted values comparing with experimental data at 350ºC Fig. 12. Predictive valus comparing with experimental data with feedstock distribution. (a) Predicted values comparing with experimental data at 320 ◦ C. (b) Predicted values comparing with experimental data at 350 ◦ C. (c) Predicted values comparing with experimental data at 380 ◦ C.

References [1] J. Ancheyta, S. Sanchez, M.A. Rodriguez, Kinetic modeling of hydrocracking of heavy oil fractions: a review, Catal. Today 109 (2005) 76–92. [2] B. Mansoornejad, N. Mostoufi, F. Jalali-Farahani, A hybrid GA–SQP optimization technique for determination of kinetic parameters of hydrogenation reactions, Comput. Chem. Eng. 32 (2008) 1447–1455. [3] J.H. Holland, Adaption in Natural and Artificial Systems, University of Michigan Press, Ann Arbor, 1975. [4] D.E. Goldberg, Genetic Algorithms in Search, Optimization and Machine Learning, Addison-Wesley, New York, 1989.

K. Wang, N. Wang / Chemical Engineering Journal 167 (2011) 228–239 [5] E. Rezaei, S. Shafiei, Heat exchanger networks retrofit by coupling genetic algorithm with NLP and ILP methods, Comput. Chem. Eng. 33 (2009) 1451–1459. [6] H. Kordabadi, A. Jahanmiri, Optimization of methanol synthesis reactor using genetic algorithms, Chem. Eng. J. 108 (2005) 249–255. [7] T. Dias Jr., L.F. Milanez, Optimal location of heat sources on a vertical wall with natural convection through genetic algorithms, Heat Mass Transfer 49 (2006) 2090–2096. [8] I.R.S. Victorino, J.P. Maia, E.R. Morais, M.R. Wolf Maciel, R. Maciel Filho, Optimization for large scale process based on evolutionary algorithms: genetic algorithms, Chem. Eng. J. 132 (2007) 1–8. [9] A.C. Papes Filho, R. Maciel Filho, Hybrid training approach for artificial neural networks using genetic algorithms for rate of reaction estimation: application to industrial methanol oxidation to formaldehyde on silver catalyst, Chem. Eng. J. 157 (2010) 501–508. [10] E.M. Faghihi, A.H. Shamekhi, Development of a neural network model for selective catalytic reduction (SCR) catalytic converter and ammonia dosing optimization using multi objective genetic algorithm, Chem. Eng. J. 165 (2010) 508–516. [11] J.L. Tao, N. Wang, DNA computing based RNA genetic algorithm with application in parameter estimation of chemical engineering processes, Comput. Chem. Eng. 31 (2007) 1602–1618. [12] K.T. Wang, N. Wang, A novel RNA genetic algorithm for parameter estimation of dynamic systems, Chem. Eng. Res. Des. 88 (2010) 1485–1493. [13] X. Chen, N. Wang, A DNA based genetic algorithm for parameter estimation in the hydrogenation reaction, Chem. Eng. J. 150 (2009) 527–535.

239

[14] J.L. Tao, N. Wang, DNA double helix based hybrid genetic algorithm for the gasoline blending recipe optimization problem, Chem. Eng. Technol. 31 (2008) 440–451. [15] X. Chen, N. Wang, Optimization of short-time gasoline blending scheduling problem with a DNA based hybrid genetic algorithm, Chem. Eng. Process. 49 (2010) 1076–1083. [16] D.P. Clark, Molecular biology: understanding the genetic revolution, Academic Press, New York, 2005. [17] R.L. Wang, K. Okazaki, An improved genetic algorithm with conditional genetic operators and its application to set-covering problem, Soft Comput. 11 (2007) 687–694. [18] L. Lin, M. Gen, Auto-tuning strategy for evolutionary algorithms: balance between exploration and exploitation, Soft Comput. 13 (2009) 157– 168. [19] Y.S. Yun, M. Gen, Performance analysis of adaptive genetic algorithm with fuzzy logic and heuristics, Fuzzy Optim. Decis. Mak. 2 (2003) 161–175. [20] V.K. Koumousis, C.P. Katsaras, A saw-tooth genetic algorithm combining the effects of variable population size and reinitialization to enhance performance, IEEE Trans. Evol. Comput. 10 (2006) 19–28. [21] S. Sanchez, M.A. Rodriguez, J. Ancheyta, Kinetic model for moderate hydrocracking of heavy oil, Ind. Eng. Chem. Res. 44 (2005) 9409–9413. [22] H. Hassanzadeh, J. Abedi, Modelling Parameter estimation of ultra-dispersed in situ catalytic upgrading experiments in a batch reactor, Fuel 89 (2010) 2822–2828.