Ind. Eng. Chem. Res. 2003, 42, 4723-4731
4723
Kinetic Parameter Estimation in Hydrocracking Using a Combination of Genetic Algorithm and Sequential Quadratic Programming P. Balasubramanian,† S. J. Bettina,† S. Pushpavanam,*,† and K. S. Balaraman‡ Department of Chemical Engineering, Indian Institute of Technology, Madras, Chennai 600 036, India, and Research and Development Centre, Chennai Petroleum Corporation Limited, Manali, Chennai 600 068, India
Hydrocracking is carried out in a multiphase trickle-bed reactor. The reactions that occur in this reactor are the cracking of high molecular weight hydrocarbons into lower molecular weight hydrocarbons followed by their hydrogenation. In a typical hydrocracker, many reactions occur simultaneously. It would be cumbersome to keep track of all of the individual reactions (literally thousands) at the molecular level. This problem is avoided by lumping the reactants together on the basis of a property like the carbon number, etc. The different lumps in the hydrocarbon mixture are characterized by the carbon number (number of carbon atoms in a molecule) or molecular weight of the components or their boiling point. A discrete lumping approach is used to model the evolution of the weight fraction distribution with time using one of these variables as the characteristic parameter. A complete description of the behavior of the lumps would necessitate consideration of all possible reactions the lumps undergo. This can be a large number. However, in a realistic situation only some of these reactions will be dominant or significant. The rate constants of these fast reactions are significantly higher than those of the slow reactions, which can be approximated to be zero. We discuss a technique based on genetic algorithm and sequential quadratic programming to determine the significant reactions and their corresponding rate constants. We consider two choices of characteristic variables (i) carbon number (or molecular weight) and (ii) true boiling point. It is shown that for the latter choice multiple solutions exist for the choice of the important reactions and their rate constants. These solutions identically satisfy the evolution of the weight fraction distribution equation. When considering the carbon number as the characteristic variable, the important reactions and their rate constants are identified uniquely. The method proposed is verified using simulated as well as experimental data. Introduction Hydrocracking is an important refinery operation aimed at the conversion of a complex heavy feedstock like vacuum gas oil (VGO) into lighter fractions such as naphtha, gasoline, and gases in the presence of hydrogen gas. The reactions that occur in hydrocracking are (i) cracking of the C-C bond of high molecular weight compounds giving rise to unsaturated low molecular weight compounds and (ii) hydrogenation of olefins, aromatic rings, sulfur compounds, nitrogen compounds, and oxygen compounds. The catalyst used in the hydrocracking process is dual-functional. A typical catalyst consists of silica-alumina with basemetal components such as Ni, W, and Mo. Here the silica-alumina portion encourages cracking activity and the noble metals such as nickel, tungsten, etc., encourage hydrogenation. The cracking reaction provides the olefins for hydrogenation, while the hydrogenation liberates the heat required for cracking. The cracking reactions are slightly endothermic, while the hydrogenation reactions are highly exothermic. The overall hydrocracking process is highly exothermic in nature. * To whom correspondence should be addressed. Tel.: +91(44)22578218. Fax: +91(44)22570509. E-mail: spush@ che.iitm.ac.in. † Indian Institute of Technology. ‡ Chennai Petroleum Corp. Ltd.
A major part of engineering involved in the design of a commercial hydrocracker is aimed at controlling the adiabatic temperature rise. This temperature rise can be controlled by the use of cold recycle gas as a quench injected at various points between a series of catalyst beds.1 Recently, the demand for middle distillates and light distillates has increased. To meet this product demand, the hydrocracking process is used to upgrade the heavier petroleum fractions. Hydrocracking yields superior quality as well as improved quantity of products as compared to other cracking processes such as fluidized catalytic cracking and thermal cracking. The formation of coke is a major drawback in catalytic cracking and thermal cracking. In the hydrocracking process, coke formation is suppressed by the hydrogenation of the unsaturated cracked products. Estimating the kinetics of the reactions in a hydrocracker is important to predict the performance of the process, to optimize it, and to control it. The kinetic model should ideally take into account all of the reactions that the components in the feedstock undergo. However, in actual practice, it is difficult to consider all individual components and the reactions involved. So, the kinetic models proposed in the past have been based on discrete lumping or continuous lumping approaches. The lumping of individual components can be based on the true boiling point (TBP) curve or on the
10.1021/ie021057s CCC: $25.00 © 2003 American Chemical Society Published on Web 08/28/2003
4724 Ind. Eng. Chem. Res., Vol. 42, No. 20, 2003
carbon number or molecular weight. Here, rather than keep track of individual molecules, the molecules with a similar carbon number, molecular weight, or boiling point are treated as reacting with a particular rate constant. Sullivan et al.2 have analyzed the chemistry of hydrocracking by investigating mixtures of several pure compounds. Discrete lumped models have been proposed to study the hydrocracking process.3-7 In these models, the various individual reactions that molecules with similar properties undergo have been lumped into an overall hydrocracking reaction. The success of this approach depends on the choice of the number of lumps. A large number of lumps results in more accurate predictions. This, however, leads to an increase in the number of kinetic parameters that describe the system. Stangeland8 suggested a simple approach for modeling of hydrocracking kinetics. He used the analogy of cracking with comminution of particles and chose the TBP of the hydrocarbon mixture as his characteristic variable. Krishna and Saxena9 modeled the hydrocracker using an axial dispersion model with two parameters. In their detailed lumped kinetics, the number of kinetic parameters involved is high. An alternative approach in the hydrocracker modeling uses a continuous lumping of hydrocarbon mixtures. Browarik and Kehlen10 considered fragmentation kinetics to describe paraffin hydrocracking. Laxminarasimhan et al.11 used the continuum theory of lumping to describe hydrocracking based on TBP. They considered a skewed Gaussian distribution function as the yield distribution to predict the evolution of the weight fraction distribution. Recently, Martens and Marin12 proposed a kinetic model for hydrocracking of VGO based on carbonium ion chemistry. A one-dimensional heterogeneous model was used to simulate the multiphase adiabatic industrial hydrocracker. They considered the effect of gas-liquid mass transfer in their reactor simulation. Alhumaizi et al.13 formulated two kinds of catalytic reaction schemes for hydrocracking of n-heptane at low temperatures. They considered the reaction mechanism based on the Langmuir adsorption isotherm. They used a one-dimensional pseudo-homogeneous model to estimate the kinetic parameters. In the present work, a discrete lumped model has been formulated for hydrocracking based on two characteristic parameters: (i) carbon number or molecular weight and (ii) the TBP. Compounds in a particular lump can crack and yield product molecules in the lower lumps. The number of reactions depends on the number of lumps chosen. In the earlier studies, researchers have assumed all of the possible reactions to be important and estimated the corresponding rate constants. In this approach the number of rate constants that have to be estimated is very large. This would also require a large amount of experimental data. Besides, the kinetic parameters that are estimated from the experimental data may not be valid outside the range of operating conditions of the experiments. It is hence desirable to have a model where the number of parameters that have to be estimated is minimal. This would increase the reliability of the model. In a typical hydrocracker, only some of the reactions undergone by the lumps may be important. Keeping this in mind, we ask the following question for a given number of significant reactions: can a set of reactions be found that best models the experimental data. For these reactions, estimate the
Table 1. Classification of Typical Lumps Based on the Carbon Number and Boiling Point Ranges lump
name of lump
carbon number range
TBP range (°C)
1 2 3 4
gases light distillate middle distillate residue
C1-C4 C5-C9 C10-C18 C18+
0-30 30-150 150-370 370+
kinetic rate constants that represent the data. The main benefit of this approach is that we increase the reliability on the model because the number of parameters to be estimated is small. In a hydrocracker, it is likely that only some of the reactions are fast and significant. For all practical purposes, the rate constants of the other slow reactions can be approximated as zero. In this work, we identify the significant reactions and their rate constants from a pool of all possible reactions that can best model a system. The reduction in the number of rate constants makes the model tractable and reliable. The important representative lumped reactions are identified using genetic algorithms (GA), and the corresponding rate constants are estimated using sequential quadratic programming (SQP). We have generated data from simulations using a known set of reactions and rate constants. We validate the method proposed on this data set. The approach is applied to experimental data of other researchers. The limitations of the proposed method are discussed. The primary objective of this paper is to demonstrate the applicability and validity of the approach proposed to determine the subset of significant reactions and their rate constants. Theoretical Models We have formulated a discrete lumped model based on carbon number as well as TBP as the characteristic parameters. In the model using carbon number as a basis, we assume that a high molecular weight hydrocarbon with a large number of carbon atoms cracks into two low molecular weight products, each with a lower number of carbon atoms. Here a Cr molecule with r carbon atoms cracks into two molecules, Cm and Cr-m, with m and r - m carbon atoms, respectively. This holds because in a reaction the total number of C atoms is conserved. Here identification of the first product lump automatically identifies the second product lump. In the second lumping approach, the TBP of pseudocomponents is used as the characterization parameter. A typical TBP curve of the hydrocarbon mixture represents the weight fraction distribution with respect to its boiling point range. Here, in general, molecules from a high TBP lump crack into product molecules of lower TBP lumps. The main difference between this and the carbon number basis is that here the product lumps are determined independently. The rth lump can give rise to products in the mth and nth lumps, where m + n * r. Typical lumps considered in the discrete lumped kinetic model based on these parameters are shown in Table 1 .4,12 The values here are only indicative of typical ranges of carbon number and TBP of various lumps. The ranges used in the calculations performed in this paper are different, and this is explained in the relevant section. Two Bases for Hydrocracking. The general stoichiometry of lumped reactions when the carbon number14,15 is the basis for the classification can be represented as
Ind. Eng. Chem. Res., Vol. 42, No. 20, 2003 4725 km,r
Cr 98 Cm + Cr-m
(1)
where r varies from 2 to n and m varies from 1 to r - 1. Here n represents the maximum number of carbon atoms present in the mixture. Here km,r is the rate constant for cracking of species with r carbon atoms into two products with m and r m carbon atoms, respectively. It signifies the rate of consumption of reactants in lump r giving products in lump m. In the cracking event using the carbon number as the basis, two products are formed, and fixing one product uniquely determines the other product. This feature also holds when the molecular weight is chosen as the characteristic property for determining the lumps. When using the TBP as the basis for the lumps, this restriction is relaxed. Here the lumps in which the products are formed during cracking do not depend on each other. The general reaction scheme when TBP is the basis can be written as ki,j,r
Cr 98 Ci + Cj
dt dcr dt
kj,ncn ∑ j)1
n
)2
(7)
r-1
kj,rcr ∑ kr,jcj - ∑ j)r+1 j)1
(8)
where r varies from 2 to n - 1
dc1
n
)2
k1,jcj ∑ j)2
(9)
Here kj,n represents the rate constant of the reaction where molecules with n carbon atoms gives rise to products with j and n - j carbon atoms. Here cj represents the molar concentration of molecules present in the lump j. Here the sum of the molar concentration n ci * constant. However, increases with time; i.e., Σi)1 n Σi)1ici ) constant, implying conservation of the carbon number. TBP Approach. The mass balance equation that describes the time evolution of the molar concentration distribution is now given by
dcn
n-1 n-1
)-
dt dcr dt
n
)2
∑ ∑ ki,j,ncn i)1 j)1
j-1
∑ ∑
(10)
r-1 r-1
kr,i,jcj -
j)r+1 i)1
∑ ∑ki,j,rcr i)1 j)1
(11)
where r varies from 2 to n - 1.
and for n odd (n ) 2M + 1)
(4)
Here the lump n is such that it has molecules containing exactly n C atoms. Similarly, when using TBP as the basis, the number of reactions involved for even n is
2M(2M + 1)(2M + 2) 6
(5)
(2M + 1)(2M + 2)(2M + 3) 6
(6)
and for odd n is
Nk )
n-1
)-
(2)
Nk ) M + 2(M - 1) + 2(M - 2) + ... + 2.1 ) M2 (3)
Nk )
dcn
dt
In general, i and j can go from 1 to r; i.e., molecules in the reactant lump r can give rise to molecules in the product lump r. However, we will assume that i and j can go from 1 to r - 1. This implies that molecules in reactant lump r can give rise to product molecules in lump r - 1 or lower. This reduces the total number of reactions that the lumps can undergo, thereby making the computations very efficient. Hence, throughout this work, r varies from 2 to n and i and j vary from 1 to r - 1. Here ki,j,r is a rate constant that determines cracking of cut r to give two products in cuts i and j, respectively. It represents the rate of consumption of reactants of lump r. When the carbon number is used as the basis and the number of lumps n considered is even (n ) 2M), then the total number of reactions (Nk) is
Nk ) M2 + M
variable in the batch reactor is the time of the reaction, while in the plug-flow reactor, it is the space time. We assume all reactions undergone by the lumps to be isothermal, first-order, and irreversible. Carbon Number. The dependency of the molar concentration distribution of the different lumps is governed by the ordinary differential equations, and it can be represented as
Model Equations The unsteady-state mass balance in a batch reactor is identical to the steady-state mass balance in an ideal plug-flow reactor if volume changes due to reaction are neglected. This is a valid assumption because the reaction occurs in the liquid phase. The independent
dc1 dt
n j-1
)2
∑ ∑k1,i,jcj j)2 i)1
(12)
Here ki,r,j represents the rate constant where lump j gives rise to product molecules in lumps i and r, respectively. Here the sum of the molar concentration n ici * constant. These increases with time and Σi)1 equations arise from a direct extension of ideas from comminution or polymer reaction engineering.15,16 Usually in hydrocracking the experimental data are available in the form of weight fraction distributions. The mass balance equations (7)-(12) are in terms of molar concentration distribution. These are rewritten in terms of wn, the weight fraction distribution, using the fact that wnF ) Mncn, where F is the mass density of the mixture and Mn is the molecular weight of lump n. We assume F is a constant because the reaction occurs in the liquid phase and volume changes are neglected. Carbon Number Basis. The time evolution of the weight fraction distribution (wn) governed by the ordinary differential equations
4726 Ind. Eng. Chem. Res., Vol. 42, No. 20, 2003
dwn
n-1
)-
dt dwr dt
kj,nwn ∑ j)1
n
)2
(13)
r-1
kj,rwr ∑ δr,jkr,jwj - ∑ j)r+1 j)1
(14)
where r varies from 2 to n - 1
dw1 dt
n
δ1,jk1,jwj ∑ j)2
)2
(15)
Here δr,j represents the fraction of weight of lump j added to lump r by virtue of reaction. This is given by Mr/Mj. In this paper this is approximated as r/j. It can be readily shown from eqs 13-15 that the sum of the n weight fractions is equal to unity, i.e., Σi)1 wi ) 1, at all time instances. TBP Basis. The evolution of the weight fractions when TBP is used as the basis is given by
dwn
dt
n
)2
)-
∑ ∑ ki,j,nwn i)1 j)1
j-1
∑ ∑
(16)
r-1 r-1
δr,i,jkr,i,jwj -
∑ ∑ki,j,rwr i)1 j)1
(17)
where r varies as 2, 3, ..., n - 1.
dt
[wpre(i,j) - wexp(i,j)]2 ∑ ∑ i)1 j)1
subject to all k > 0 (19)
j)r+1 i)1
dw1
ntime n
error )
n-1 n-1
dt dwr
This represents a population of three strings. Each string has 10 elements signifying the 10 reactions. Because four of them are 1, each string represents four significant or fast reactions. These strings represent the initial guess of the algorithm. The actual algorithm uses a population of eight strings. 2. The objective function used in the optimization algorithm was to minimize the error between experimental weight fractions and the model prediction results, i.e., minimize
n j-1
∑ ∑δ1,i,jk1,i,jwj j)2 i)1
)2
(18)
Here δr,i,j is the fraction of weight of lump j that goes to a product in lump r. In this work we assume δr,i,j ) r/i + r. The above system of equations also ensures that the sum of the weight fractions is equal to unity at all instances of time. Kinetic Parameter Estimation. The primary aim of this paper is to demonstrate how kinetic parameters of the subset of significant reactions in hydrocracking can be estimated. We have used an optimization algorithm consisting of a combination of GA and SQP. GA identifies the subset of significant reactions, and SQP estimates the rate constants. We now talk briefly about GA and explain the principles on which it is based. Genetic algorithm (GA). GA17-20 is a random search and optimization technique. In this work it is used to determine the significant reactions in the hydrocracking process. The steps involved in the GA are as follows: 1. An initial population, consisting of points in the search space, is randomly created. This is usually a binary representation. The set of reactions for the system being investigated is represented as an m × n dimensional matrix. Here m denotes the population size and n denotes the total number of possible reactions that have to be considered. The elements of each row vector can be either 0 or 1. The former signifies that the reaction is not important while the latter signifies that the reaction is important. For example, when using the TBP as the basis and a population size of three, the 3 × 10 matrix represents three different possible reaction schemes.
Here wpre is the model prediction weight fraction data and wexp is the experimental data. The objective function for each string was evaluated using SQP. Thus, SQP determines the rate constants for a choice of significant reactions. The fitness of each string in the population is evaluated through the objective function. A high value is assigned to a string, which performs better. For problems that minimize the error, a good candidate for a fitness function is
fitness function (F) ) 1/(1 + objective function) (20) 3. A reproduction (selection operator) on the population was performed using the roulette wheel selection scheme. This roulette wheel selection scheme can be simulated using the fitness value Fi of all strings. The probability of selecting string Pi can be calculated as18
Pi )
Fi n
(21)
Fj ∑ j)1 The cumulative probability is calculated, and this determines an interval between 0 and 1. We generate a random number. The location of the random number in the cumulative probability range identifies the string, which is copied to the mating pool. A string with high fitness value is selected more frequently and is copied into the mating pool more often. In this step, no new strings are produced. 4. A search is performed by creating a new string from the selected string in the mating pool through GA operators such as crossover and mutation. In the crossover operation, two consecutive strings are picked up from the mating pool and some portion of the strings are exchanged between the strings. A single-point crossover operator18 is performed by randomly choosing a crossing site along the string and by exchanging all bits on the right side of the crossing site with each other as
Ind. Eng. Chem. Res., Vol. 42, No. 20, 2003 4727 Table 2. Parameters Used in GA Implementation parameter
carbon number
TBP
total number of reactions population size crossover probability mutation probability maximum number of generations
9 8 0.8 0.1 50
10 8 0.8 0.1 50
The mutation operator creates a point in the neighborhood of the current point, thereby achieving a local search around the current solution. The mutation operator changes 1 to 0 and vice versa with small mutation probability18
5. The number of 1’s can be more than or less than the specified number of the significant fast reactions in each string. This is corrected randomly to get the required number of significant reactions. The new strings calculated form the new mating pool. 6. Steps 2-5 are repeated until the population converges, i.e., all strings are identical or the specified number of generations (iteration) is reached. The parameters used in the GA are depicted in Table 2. The kinetic parameters are estimated using the “constr” function in the optimization toolbox of Matlab. In the proposed algorithm, GA and SQP are interconnected. Here first GA is initiated using a population of random strings. GA searches for the important reactions. SQP estimates the required rate constants in the search space. This is used to determine the fitness of each string in the GA. This is then used to generate new strings. So, reaction identification and parameter estimation occur simultaneously in the proposed optimization algorithm. Results and Discussion We now discuss the applicability of the proposed algorithm to different cases. We first discuss kinetic parameter estimation using simulated data. The objective here is to validate the proposed method under ideal conditions. Carbon Number Basis. The behavior of the reactor was simulated assuming that the maximum number of carbon atoms in the compounds making the feed mixture is 6. Each lump consists of molecules with a fixed number of carbon atoms. Thus, lump j consists of molecules with j carbon atoms. The total number of possible cracking reactions is 9. We assume only four reactions to be fast, i.e., significant. The rate constants of the other slow reactions are set to zero. The rate constant vector (k12, k13, k14, k24, k15, k25, k16, k26, k36) is assumed to be (0.1, 0, 0.3, 0, 0, 0.4, 0, 0.5, 0). This set of parameters was used to simulate the behavior of the reactor. Equations 13-15 governing weight fraction distribution from an initial weight distribution of (0, 0, 0, 0.1, 0.4, 0.5) were integrated, and the corresponding weight fraction distribution at different time
Table 3. Simulated Data Obtained Using the Carbon Number as a Basis k12 0.1
Assumed Rate Constants (h-1) k14 k24 k15 k25 k16 0.3 0 0 0.4 0
k13 0
k26 0.5
k36 0
simulated weight fraction
time (h)
lump 1
lump 2
lump 3
lump 4
lump 5
lump 6
0 0.5 1.0 1.5 2.0
0 0.014 0.037 0.062 0.087
0 0.115 0.183 0.220 0.239
0 0.113 0.211 0.292 0.357
0.100 0.186 0.206 0.194 0.168
0.400 0.268 0.180 0.121 0.080
0.500 0.303 0.184 0.112 0.068
Initial Guess Used in the GA reaction scheme
string no.
k12
k13
k14
k24
k15
k25
k16
k26
k36
1 2 3 4 5 6 7 8
0 1 0 1 0 1 1 0
1 0 0 0 1 0 0 1
0 1 1 0 1 0 1 0
0 0 1 0 0 1 0 1
1 0 0 1 0 0 1 0
0 1 0 0 1 0 0 1
0 0 1 0 1 0 1 0
1 1 0 1 0 1 0 0
1 0 1 1 0 1 0 1
Table 4. Estimated Rate Constants Assuming Different Numbers of Significant Reactions Using the Carbon Number as a Basisa no. of significant reactions
k12
k13
k14
k24
k15
k25
k16
k26
k36
2 3 4 5 6 7
0 0 0.10 0.10 0.10 0.10
0 0 0 0 0 0
0 0.32 0.30 0.30 0.30 0.30
0 0 0 0 0 0
0 0 0 0 0 0
0.57 0.38 0.40 0.40 0.40 0.40
0 0 0 0 0 0
0.30 0.50 0.50 0.50 0.50 0.50
0 0 0 0 0 0
estimated rate constant value (h-1)
a The ideal case is where the rate constant of the slow reactions is zero.
instances 0.5, 1, 1.5, and 2 h was obtained. This weight fraction distribution was assumed to be experimental data and is shown in Table 3. The initial guess of the GA strings is also shown here. The optimization algorithm was used to determine the important reactions and the rate constants. The results are shown in Table 4. Here we show how the algorithm performs when the number of important reactions is varied. When the number of reactions is less than 4, the algorithm identifies reactions that are part of the significant set as being important. For n ) 4, the algorithm identifies the important reactions correctly, along with the rate constants. As we increase the number of important reactions to more than 4, the algorithm identifies the same set of four important reactions. Thus, the algorithm is able to accurately determine the set of important reactions and the rate constants for this test case. In a realistic situation, the reaction rates of the slow reactions will be nonzero but small. We want to determine the performance of the algorithm when the reaction rates of the slow reactions are nonzero and small. All slow reactions are assigned a uniform low value. For k ) 0.01, the slow reactions are about 1/30 of the fast reactions. For this the algorithm detects the four important reactions accurately. It also identifies the rate constants of these reactions reasonably accurately. When we evaluate the rate constant for five and six
4728 Ind. Eng. Chem. Res., Vol. 42, No. 20, 2003 Table 5. Influence of Slow Reactions in the Optimization Algorithm Using the Carbon Number as a Basisa no. of significant reactions
a
estimated rate constant value (h-1) k12
k13
k14
k24
k15
k25
k16
k26
k36
4 5 6
Slow Reaction Rate Constant Value ) 0.01 0.10 0 0.30 0 0 0.40 0 0.51 0 0.10 0 0.29 0 0 0.41 0.02 0.50 0 0.10 0.01 0.31 0 0 0.40 0 0.51 0.01
4 5 6
Slow Reaction Rate Constant Value ) 0.05 0.10 0 0.30 0 0 0.40 0 0.56 0 0.19 0 0.28 0.15 0 0.41 0 0.57 0 0.12 0.06 0.32 0 0 0.41 0 0.56 0.04
4 5
Slow Reaction Rate Constant Value ) 0.06 0 0 0 0.22 0.23 0.67 0.58 0 0 0 0.07 0.44 0 0.12 0.29 0 0.58 0
The rate constant of the slow reaction is nonzero.
important reactions, the program identifies the four significant reactions accurately. These results are depicted in Table 5. For a rate constant of 0.05 for the slow reactions, the performance of the algorithm is shown in Table 5. Here again the four fast reactions are identified accurately when we seek more than four rate constants. The algorithm hence is robust and works well when the slow reaction rates are around 15% of the fast rates. When the rate constant for slow reactions is 0.06, the algorithm cannot identify the four important reactions accurately. TBP Approach. The behavior of the reactor using the TBP as the characteristic parameter was simulated assuming that the number of lumps is 4. The total number of possible cracking reactions is 10. Here we preclude the possibility that lump j on cracking can give product molecules in lump j. We assumed only four reactions to be fast and significant. The rate constants of the other slow reactions are set to zero. The initial weight fraction distribution considered is (0, 0.1, 0.4, 0.5). The assumed rate constant vector is (0, 0, 0.4, 0, 0.3, 0, 0.5, 0, 0, 0.6) corresponding to (k112, k213, k223, k113, k314, k324, k334, k214, k224, k114). This rate constant vector was used to simulate the performance of the reactor. The equations governing the weight fraction distribution are
dw4/dt ) -1.7w4
(22)
dw3/dt ) 0.95w4 - 0.4w3
(23)
dw2/dt ) 0.4w3
(24)
dw1/dt ) 0.75w4
(25)
These were integrated, and the corresponding weight fraction distribution at the time instances 0.5, 1, 1.5, and 2 h was obtained. This weight fraction distribution was assumed to be experimental data and is depicted in Table 6. The optimization algorithm discussed above was used to determine the important reactions and their rate constants. Here the algorithm converged to a rate constant set that was different from that assumed. This indicates that the system has multiple solutions. The assumed rate constant set for generating the experimental data and a solution set are depicted in Table 6. It can be easily verified that the estimated rate con-
Table 6. Simulated Data and Estimated Rate Constants Assuming Four Reactions To Be Significant Using the TBP Modela k112 0
k213 0
Assumed Rate Constants (h-1) k223 k113 k314 k324 k334 k214 0.40 0 0.30 0 0.50 0
k224 0
k114 0.60
simulated weight fraction
time (h)
lump 1
lump 2
lump 3
lump 4
0 0.5 1.0 1.5 2.0
0 0.126 0.180 0.203 0.213
0.100 0.190 0.282 0.366 0.438
0.400 0.471 0.446 0.392 0.332
0.500 0.214 0.091 0.039 0.017
k112 0
k213 0
Estimated Rate Constants (h-1) k223 k113 k314 k324 k334 k214 0.4 0 0 0 0.95 0
k224 0
k114 0.75
a The ideal case is where slow reactions have rate constants of zero.
Table 7. Influence of Slow Reactions in the Optimization Algorithm Using the TBP Approacha k112
Different Possible Reaction Rate Constants k213 k223 k113 k314 k324 k334 k214 k224
k114
Slow Reaction Rate Constant ) 0.01 0.01 0.03 0
0 0.01
Assumed Rate Constants (h-1) 0.40 0.01 0.30 0.01 0.50 0.01
0.01
0.60
0 0.03
Estimated Rate Constants (h-1) 0.43 0 0.65 0 0 0 0.38 0 0 0 0.98 0
0 0
0.45 0.76
Slow Reaction Rate Constant ) 0.04 0.04 0 0 a
0.04
Assumed Rate Constants (h-1) 0.40 0.04 0.30 0.04 0.54 0.04
0.04
0.60
0.27 0.27
Estimated Rate Constants (h-1) 0 0 0.69 0 0 0.22 0 0 0.40 0.37 0 0
0 0
0.09 0.38
Slow reactions have a nonzero rate constant.
stants satisfy eqs 22-25 identically. Hence, this method yields multiple solutions when we use TBP as a basis. Here the subset of significant reactions cannot be uniquely identified. As mentioned earlier, in a realistic situation the reaction rates of the slow reaction will be nonzero and small. We would like to see if the algorithm can uniquely identify the fast reactions in the presence of nonzero but small values of the rate constants for the slow reactions. Essentially this would imply that the nonzero values of the slow rate constants can break the multiple solutions. All slow reactions are assumed to have a uniform low value of k. For k ) 0.01 and 0.04, the results are shown in Table 7. Here again the system possesses multiple solutions. We repeated this analysis by using simulated data that consider different sets of reactions to be significant. In all cases we obtained multiple solutions for the set of significant reactions. Hence, when we use TBP as a basis, we cannot discriminate and identify which of the multiple solutions for the significant reactions is actually responsible for generating the experimental data. The multiplicity of the solutions of optimization can be understood by looking at the system of linear equations that have to be solved to obtain k. Here we assume that the rate constants of the slow reactions are zero. For the set of parameters considered, the nonzero rate constants must satisfy
Ind. Eng. Chem. Res., Vol. 42, No. 20, 2003 4729 Table 8. Estimated Rate Constants for Different Catalyst Experimental Data of Ali et al.1 Using the Carbon Number Approach reaction temperature (°C) β + ASA/NiMo catalyst
β + ASA/NiW catalyst
β/NiW catalyst
rate constant (h-1)
380
395
300
410
380
395
400
410
380
395
400
410
k12 k13 k14 k24 k15 k25 k16 k26 k36
0.27 0 0 2.01 0.22 0.12 0 0 0.69
0 0 0.31 2.21 0.36 0.16 0 0.72 0
0.11 0.83 0 0 0.002 0.73 0 0 1.94
0.12 0 0 3.94 0.67 0.33 2.76 0 0
0 0.50 0.09 0 0.05 0.29 0.31 0 0
0 0.80 0.04 0 0.002 0.64 0 0 1.38
0 0 0.27 2.99 0.57 0.25 0 0.88 0
0.12 0 0 5.41 1.02 0.41 0 2.92 0
0 0.56 0.12 0 0.04 0.25 0 0.33 0
0 0 0 3.41 0.43 0.31 0.90 0 0
0 0.78 0 0.35 0 0.77 0.72 0.20 0
0 0 0.04 4.27 0.81 0.41 0 2.22 0
2k134 + k334 + k114 ) 1.7
(26)
1.5k134 + k334 ) 0.95
(27)
0.5k134 + k114 ) 0.75
(28)
This system of equations can be written as
Ak ) b
(29a)
Here
( ) ( ) ( )
k134 2 1 1 1.7 A ) 1.5 1 0 , k ) k334 , b ) 0.95 k114 0.5 0 1 0.75
(29b)
The rank of A is 2, which is less than 3, the order of the matrix. Consequently, the system has an infinite number of solutions. This explains the cause for the multiple solutions. For the case where the carbon number is the basis, we can show that the rank of A is 4, which is equal to its order, and so the system has a unique solution. Application on Experimental Data. In the earlier section, we have verified and validated the applicability of the proposed procedure on simulated data. We now test the applicability of the parameter estimation algorithm on the experimental data of Ali et al.1 Ali et al. conducted experiments in a high-pressure fixed-bed reactor at four different temperatures (380, 395, 400, and 410 °C) with five different catalysts. The operating conditions used are a pressure of 100 kg/cm2, a liquid hourly space velocity of 1 h-1, and a hydrogen-tohydrocarbon ratio of 5900 scfb. Desulfurized VGO was utilized as the feed for the evaluation of the catalysts. We now discuss how the kinetic parameters were estimated using the proposed optimization algorithm on their experimental data. Carbon Number Basis. We consider the parameter estimation problem using the carbon number as the basis. The number of lumps considered is 6 (equal to that considered by Ali et al.1). The initial weight fraction distribution determined experimentally was (0, 0, 0, 0.089, 0.847, 0.064), and its corresponding lumps are gases, naphtha, kerosene, gas oil, VGO, and vacuum residue (VR), respectively. The total number of possible reactions is, hence, 9. The kinetic constant vector is assumed to be of the form (k12, k13, k14, k24, k15, k25, k16, k26, k36). We determined the kinetic constants using the experimental data of catalyst β + ASA/NiMo at a reactor temperature of
380 °C. The parameter estimation algorithm efficiently converges with low error (1.8 × 10-7) when five reactions are assumed to be significant. The estimated rate constants values are 0.27, 0, 0, 2.01, 0.22, 0.12, 0, 0, and 0.69. Here the kinetic constants k13, k14, k16, and k26 are insignificant (i.e., slow reaction rate constants with zero values). We repeated the parameter estimation algorithm for other reaction temperatures of 395, 400, and 410 °C. Here also the algorithm converges with low error (1.7 × 10-7) when five reactions are considered to be significant. The estimated rate constant values are depicted in Table 8. We see that the algorithm identifies different reactions as being significant at different temperatures. Besides, the kinetic constants are not monotonically increasing functions of temperature. These apparently contradictory results arise because each lump measured experimentally contains molecules such that the number of C atoms lie in a range. However, in the formulation of the optimization problem, we assume that the number of C atoms is fixed and is unique for a lump. Consequently, in the experimental data molecules from lump 4 may give molecules in lump 3 or 2. The algorithm has neglected this. It is also likely that at different temperatures the mechanisms of the reactions are different. This can give rise to certain rate constants becoming zero or attaining high values as the temperature changes. A temperaturedependent activation energy can also give rise to this kind of behavior. It must be borne in mind that the temperature ranges over which the experiments were carried out are very narrow. Consequently, variations in the data reported within the standard deviation can also give rise to this anomalous behavior. We applied the proposed kinetic parameter estimation algorithm on the experimental data obtained from other catalysts (β + ASA/NiW and β/NiW). The estimated significant rate constants are depicted in Table 8. Here again the rate constants do not monotonically increase with temperature. However, for each data set, the weight fraction predicted by the model agrees very well with experimental data. TBP Basis. Ali et al.1 have depicted the weight fraction distribution for 6 lumps in their work. When 6 lumps are considered, the total number of possible reactions using the TBP as the basis is 56. This increases the computational complexity. To reduce the computational intensity, we have regrouped the lumps and considered only 4 lumps. The lumps are classified based on the boiling point range of the hydrocarbon mixture, i.e., (i) gases (0-30 °C), (ii) light distillate (30150 °C), (iii) middle distillate (150-343 °C), and (iv)
4730 Ind. Eng. Chem. Res., Vol. 42, No. 20, 2003 Table 9. Estimated Rate Constants for Experimental Data on β + ASA/NiMo and β + ASA/NiW Catalysts of Ali et al.1 Using the TBP as a Basisa reaction temperature (°C) β + ASA/NiMo catalyst 380
395
β + ASA/NiW catalyst
400
410
380
395
400
410
rate constant (h-1)
1
2
1
2
1
2
1
2
1
2
1
2
1
2
1
2
k112 k213 k223 k113 k314 k324 k334 k214 k224 k114
0.59 0 0 0 0.002 0.16 0 0 0.36 0
0 0 0.127 0.002 0 0 0.54 0 0 0.13
0 0 0.25 0 0.25 0 0 0 0.52 0.04
0 0.42 0 0.02 0 0.50 0 0 0 0.06
0 0 0 0.32 0 0.42 0 0.17 0.48 0
0.24 0.32 0 0 0 0.56 0 0 0.55 0
0 0 1.89 0 0.76 0.09 0.29 0 0 0
0 0 1.61 0 0.76 0 0.16 0 0.31 0
0 0 0.18 0.27 0 0.32 0 0 0 0
0 0 0.98 0.15 0.06 0 0.51 0 0 0
0.10 0 0 0 0.22 0 0 0 0.83 0
0.17 0 1.71 0 0.19 0 0.90 0 0 0
0 0.46 0 0 0 0 0.84 0 0.73 0.09
0 0 0 0.07 0.26 0 0 0 1.04 0.09
0 0.53 0.94 0 0.53 0.94 0 0 0 0
0 0 1.37 0 0.39 0.54 0 0.54 0 0
a
Two values are shown for each data set. Solutions are indicated by 1 and 2.
residue (>343 °C). Here kerosene (150-250 °C) and gas oil (250-343 °C) are regrouped as middle distillate and VGO (343-520 °C) and VR (>520 °C) are regrouped as residue. This regrouping reduces the number of lumps and hence the total number of reactions to be considered. The total number of reactions involved when we considered 4 lumps (n ) 4) is 20. To reduce this even further, we neglected reactions in which reactant molecules from a lump give rise to product molecules in the same lump. With this assumption the maximum number of possible reactions is reduced to 10. We have estimated the kinetic constants for the experimental data of catalyst β + ASA/NiMo using the proposed algorithm. When TBP is the basis, the algorithm converges with multiple kinetic constants for the different reaction temperatures. These estimated rate constants when used in the mass balance equations result in the prediction of weight fraction data that match the experimental values very closely. Here also the rate constants are not monotonically increasing functions of temperature. The estimated rate constants obtained using the experimental data of this catalyst for different reaction temperatures are depicted in Table 9. We have also applied the parameter estimation algorithm on the experimental data obtained using the β + ASA/NiW catalyst. Here also we obtain multiple solutions for the set of significant reactions and their rate constants. The estimated rate constants demonstrating the multiplicity are depicted in Table 9. The algorithm proposed yields multiple solutions for the significant reactions and their rate constants when TBP is used as a basis. We now discuss a method to discriminate between these multiple solutions. We illustrate this by applying the algorithm on experimental data reported by Bennett and Bourne.21 They have reported results for four residence times (0.383, 0.952, 1.724, and 2.5 h). The feed had compositions of 0, 0, 0.41, and 0.59 corresponding to the TBP lumps of gases (0-30 °C), gasoline (30-150 °C), middle distillate (150-370 °C), and residue (370 °C+). Because the number of lumps is 4, the total number of reactions is 10. For a residence time of 2.5 h, the weight fraction data obtained were (0.3, 0.26, 0.39, 0.05). The algorithm converges with an error of 1.6 × 10-11 when four reactions were considered significant. The estimated rate constants are (0, 0.02, 0, 0.25, 0, 0, 0.56, 0, 0.43, 0) for (k112, k213, k223, k113, k314, k324, k334, k214, k224, k114).
Figure 1. A comparison of model predictions with experimental data using two identified sets of kinetic constants for various residence times: (a) from eqs 30-33; (b) from eqs 34-37. Symbols: experimental data of Bennett and Bourne.21 Solid line: model predictions.
The generating equations for the weight fractions using these rate constants are
dw4/dt ) -0.99w4
(30)
dw3/dt ) -0.29w3 + 0.56w4
(31)
dw2/dt ) +0.43w4 + 0.03w3
(32)
dw1/dt ) +0.26w3
(33)
The algorithm also converges to the set of estimated rate constants (0.14, 0, 0.27, 0, 0.35, 0, 0, 0, 0, 0.29) with an error of 8.9 × 10-11. With these sets of parameters, the weight fraction evolution is governed by
Ind. Eng. Chem. Res., Vol. 42, No. 20, 2003 4731
dw4/dt ) -0.99w4
(34)
dw3/dt ) -0.27w3 + 0.525w4
(35)
dw2/dt ) -0.135w2 + 0.27w3
(36)
dw1/dt ) 0.465w4 + 0.135w2
(37)
It can be seen that the evolution equations are hence different for the two sets of kinetic parameters. Both sets of parameters result in the same predictions at t ) 2.5 h, but their predictions at other residence times will be different. We exploit this and compare the predictions at these other residence times with experimental data. This is shown in parts a and b of Figure 1. The first set of rate constants predicts the general trend of the weight fraction evolution better. For instance, the evolution of the distribution of gases and gasoline fraction is better in Figure 1a than in Figure 1b. Hence, the first set of rate constants can be preferred over the second set by investigating the performance at different residence times. Conclusions In this work we have discussed a generalized discrete lumped model for hydrocracking based on the carbon number as well as TBP as characteristic parameters. To make the model robust, we assume only a small subset of reactions out of the total possible reactions to be significant. The significant reactions involved are identified, and their corresponding rate constants are estimated using a combination of GA and SQP. The optimization algorithm proposed was verified using the simulated as well as experimental data published in the literature. Using the carbon number as the basis, the optimization algorithm converges to a unique solution for the set of significant reactions. Using TBP as the basis, the algorithm shows multiple solutions. We can discriminate between these multiple solutions if the reaction rate constants of some of the significant reactions are known a priori or by comparison of the experimental results with the predictions at different residence times of each solution. Acknowledgment The authors thank Chennai Petroleum Corp. Ltd. (CPCL), Chennai, India, for their financial support to carry out this research work at IIT Madras, Chennai, India. Literature Cited (1) Ali, M. A.; Tasumi, T.; Masuda, T. Development of heavy oil hydrocracking catalyst using amorphous silica-alumina and zeolite as catalyst supports. Appl. Catal. A 2002, 233, 77.
(2) Sullivan, R. F.; Egan, C. J.; Langalois, G. E. Hydrocracking of Alkyl Benzene ans Polycyclic Aromatics on Catalyst; Evidence for cyclization of side chains. J. Catal. 1964, 3, 183. (3) Stangeland, B. E.; Kittrell, J. K. Jet Fuel selectivity in hydrocracking. Ind. Eng. Chem. Process Des. Dev. 1972, 11, 16. (4) Zhorov, Yu. M. J.; Panchenkov, G. M.; Tatarintseva, G. M.; Kumin, S. T.; Zen kovskii, S. M. Chemical Scheme and Structure of Mathematical Description of Hydrocracking. Int. Chem. Eng. 1971, 11 (2), 256. (5) Koseoglu, O. R.; Phillips, C. R. Kinetics and Product yield distributions in the CoO-MoO3/Al2 ozone catalyzed hydrocracking of Athabasca Bitumen. Fuel 1988, 67, 1411. (6) Callegas, M. A.; Martinez, M. T. Hydrocracking of a Maya Residue: Kineties and products yield distribution. Ind. Eng. Chem. Res. 1999, 38, 3285. (7) Korre, S. C.; Klein, M. T. Hydrocracking of Polyaromatics hydrocarbons, Development of rate law through inhibition studies. Ind. Eng. Chem. Res. 1997, 36, 2041. (8) Stangeland, B. E. Kinetic model for prediction of hydrocracker yields. Ind. Eng. Chem. Process Des. Dev. 1974, 13 (1), 72. (9) Krishna, R.; Saxena, A. K. Use of axial dispersion model for kinetic description of hydrocracking. Chem. Eng. Sci. 1989, 44, 703. (10) Browarik, D.; Kehlen, H. Hydrocracking process of n -alkanes by continuous kinetics. Chem. Eng. Sci. 1994, 49 (6), 923. (11) Laxminarasimhan, C. S.; Verma, R. P.; Ramachandran, P. A. Continuous lumping model for simulation of hydrocracking. AIChE J. 1996, 42 (9), 2645. (12) Martens, G. G.; Marin, G. B. Kinetics for hydrocracking based on structural classes: Model development and Application. AIChE J. 2001, 47 (7), 1607. (13) Alhumaizi, K. I.; Akhmedov, V. M.; Al-Zahrani, S. M.; Alkhowaiter, S. H. Low-temperature hydrocracking of n-haptane over Ni-supported catalyst: Study of global kinetics. Appl. Catal. A 2001, 219, 131. (14) Kehlen, H.; Ratzsch, M. T.; Bergmann, O. J. Continuous kinetics of first order degradation reaction in poly disperse mixture. Chem. Eng. Sci. 1988, 43 (3), 609. (15) McCoy, B. J.; Wang, M. Continuous-mixture fragmentation kinetics: Particle size reduction and molecular cracking. Chem. Eng. Sci. 1994, 49 (22), 3773. (16) McCoy, B. J.; Madras, G. Degradation kinetics of polymer in solution: Dynamics of molecular weight distributions. AIChE J. 1997, 43 (3), 802. (17) Goldberg, D. E. Genetic algorithm in Search, optimization and machine learning; Addison-Wesley: Reading, MA, 1989. (18) Deb, K. Optimization for Engineering Design: Algorithms and Examples; Prentice Hall: New Delhi, India, 1995. (19) Bagchi, T. P. Multiobjective Scheduling by Genetic algorithm; Kluwer Academic Publisher: Boston, MA, 1999. (20) Costa, L.; Oliveira, P. Evolutionary algorithms approach to the solution of mixed integer nonlinear programming problems. Comput. Chem. Eng. 2001, 25, 257. (21) Bennett, R. N; Bourne, K. H. Hydrocracking for middle distillatesA study of Process Reactions and Corresponding Product Yields and Qualities. ACS Symposium on Advances in Distillate and Residual Oil Technology; American Chemical Society: Washington, DC, 1972.
Received for review December 30, 2002 Revised manuscript received June 27, 2003 Accepted June 27, 2003 IE021057S
