OPTIMIZATION OF KINETIC LUMPING MODEL PARAMETERS TO IMPROVE PRODUCTS QUALITY IN THE HYDROCRACKING PROCESS

In this study a typical continuous lumping model with five parameters has been used for kinetic modeling of thermal and catalytic hydrocracking. Model parameters have been optimized according to experimental product distributions using a particle swarm optimization (PSO) algorithm. Experimental data from the hydrocracker setup have been employed to validate the proposed model. In this setup hydrogen and vacuum gasoil feed were introduced from the top of a vertical reactor and, after passing through a catalyst bed, the liquid and gas products were separated and analyzed. Temperature of the reactor was adjusted in the range of 440-470oC for thermal hydrocracking, and 410-430oC for catalytic hydrocracking. Liquid hourly space velocities (LHSV) were in the range of 0.5-1.5 feed flow rate per catalyst volume in both sets of experiments. Results of optimization showed that the parameters were only temperature dependent. The comparison between model results and experimental data indicates that the model is capable of predicting product yield with maximum errors of 0.986 and 0.041 for RMSE and AARE values, respectively.


INTRODUCTION
Hydrocracking is a catalytic conversion process for producing lighter petroleum fractions such as diesel, kerosene, naphtha and gases from highboiling constituent hydrocarbons in petroleum such as vacuum gasoil (VGO) feedstock, in the presence of hydrogen-rich gas.Normally, the hydrocracking process is carried out in a high pressure trickle bed reactor (Kumar and Balasubramanian, 2009).
Various reviews have been reported on hydrocracking technology in the literature.Different aspects of hydrocracking have been discussed, such as types of hydrocracking, process optimization, catalysis, catalyst acidity, pore diffusion, catalyst poisons and reactor type (Ramakanta et al., 2015;Angeles et al., 2014;Ahmed et al., 2013;Valavarasu et al., 2003;Qader and Hill, 1969;Sadighi and Ahmad, 2013;Matos and Guirardello 2002).These authors also emphasized that the available literature on the kinetics of hydrocracking is scarce, and that kinetic aspects and structural modeling needs to be studied in detail.
In the hydrocracking process, a kinetic model which can predict the performance of the process within the operating severity is very important.Moreover, this model could be extended beyond the experimental *Corresponding author.E-mail address: shokris@ripi.irdomain, for investigating the performance behavior ( M oghadassi et al., 2011; Rana et al., 2007; Alhajree et  al., 2011, Zhou et al., 2011).
Various kinetic models have been proposed by several researchers (Froment, 2005;Botchway et al., 2004;Browning et al., 2016;Ancheyta et al., 2005;Elizalde et al., 2010;Puron et al., 2014) to study the product distribution in hydrocracking.A commonly used approach is to consider the mixture in terms of selected lumps, which can be specified in terms of structural characteristics, such as boiling point.This approach is attractive for the kinetic modeling of complex mixtures because of its simplicity (Ancheyta et al., 1999).Depending on the mixture characteristic function form; there are two approaches for lumping models, which are continuous and discrete approaches (Canan and Arkun, 2012).In the literature, different lumping models have been reported in order to analyze the kinetics of hydrocracking reactions (Govindhakannan and Riggs, 2007;Khorasheh and Zainali, 2001;Adam et al., 2012;Stangeland, 1974;Chou and Ho, 1988;Stangeland and Kittrell, 1972;Becker et al., 2015;Elizalde et al., 2016).Becker et al. (2016) compared a continuous lumping model with a single events microkinetic model.They concluded that the continuous lumping model provides better results for the macroscopic effluent characteristics.Also, the single events model is far more complex and computationally expensive than the continuous lumping model.Furthermore, Elizalde et al. (2009) reported that the continuous kinetic lumping model versus discrete lumped model is able to predict the whole distillation curve of hydrocracked products.Moreover, the prediction of product composition can be performed with different boiling point ranges.In this work, a typical continuous lumping model with five parameters was employed.The main advantages of the lumping technique are its easy computational implementation and the small amount of data required for parameter estimation.
One of the major challenges in continuous lumping modeling is optimization of parameters (Elizalde and Ancheyta, 2011).In this study, the model parameters have been optimized according to experimental product distributions by a PSO algorithm.All of the experimental results have been obtained under limited operating conditions in the hydrocracker setup.
The present work has used the continuous lumping concept for kinetic modeling of hydrocracking of a typical refinery vacuum gasoil.The considered five parameter model used TBP and reactivity (k) of the species as continuous variables for hydrocarbon mixture characterization.Next, Optimization of parameters was performed using a PSO algorithm as a novel optimization method.Finally, optimized parameters were used to predict product yields at various severities.

Continuous Lumping Concept (Continuous Lumping Model Formulation)
The continuous lumping concept was first proposed by De Donder (1931).This concept is used when the mixtures include many different components and their difference is partially detectable.Under these conditions, a continuous variable, such as boiling point, molecular weight or other specifications can be used to describe mixture composition.In order to apply a continuous model, a continuous distribution function for the reactant mixture based on one of their properties needs to be defined.During the reaction, the distribution function changes.The aim of modeling is to investigate the distribution function variation and to obtain its final values.The continuous kinetic model is applicable for hydrocracking feed and products.
In this study, a typical proposed model based on continuous lumping for hydrocracking of petroleum fractions has been applied (Laxminarasimhan et al., 1996).The model considered the dimensionless boiling point (θ) and reactivity (k) as continuous variables for components identification.The dimensionless boiling point for component i is defined as follows: (1) where TBP(l) and TBP(h) indicate the lowest and highest possible boiling point in the reactant mixture, respectively.Usually experimental results of hydrocracking have been expressed in the form of cumulative weight percent (cum wt%) of components versus true boiling point(TBP).On the other hand, material balance and reactor performance equations have been written in the form of concentration (C) or non-cumulative weight percent versus reactivity (k) of components.Thus, for solving material balance equations and coordinate transformation, the continuous model considered a simple power law form relation between k i and θ i as below: (2) Here k max is the reactivity of the species with the highest boiling point, corresponding to θ =1, whereas α is a parameter of model.Both k max and α are positive constants.Also this equation for θ = 0 yields the value of k equal to zero.In fact, the cracking rate constant for the lightest species is equal to zero.Chou et al. (1988) described the coordinate transformation from the i axis (discrete mixture index) to the k axis (continuous mixture variable).They used a component type distribution function D(k) as the correction term for scale similarity along the axes.
(3) D(k) was considered as a Jacobian in coordinate transformation from i to k as follows: (4) Assuming that the total number of components in the mixture is equal to N and (N→ ∞), then di/ dθ can be replaced by N. Thus, the component type distribution function D(k) can be defined as follows: (5) The mass balance equation as a function of reactivity of species i (k i ) and residence time t (or alternatively the inverse of space velocity, 1/LHSV), can be derived as: (6) In the above equation i=1 indicates the lightest component and i=N represents the heaviest component in the mixture.Furthermore, it is assumed that the model is one dimensional and homogeneous and neglects mass transfer resistances in boundary layer and pores of catalyst.The left side of the equation is the accumulation rate of species k i .The first term on the right side represents the consumption reaction rate of k i species during cracking to lower boiling point species, which assumed first order.The second term shows the generation rate of k i species from heavier boiling point species with reactivity k i < K < k max .In the second term p(k i ,K) is the yield distribution function, that represents k i species formation from K species cracking reactions.The function p (k i , K) must have the following conditions: (1) The value of p (k i , K) for k i = K or k i > K should be zero because the component with reactivity ki cannot be generated from itself or lighter components during hydrocracking.
(2) The function p(k i , K) must satisfy the following mass balance equation: (3) Function p (k i , K) should have a finite nonzero value when k i = 0. Experimental findings show that when a component with reactivity K cracks, the lightest component in the mixture (k i =0) is produced in traces.( 4) The value of p (k i , K) must always be positive.
(5) A skewed Gaussian type distribution function such as Equation ( 8) can be considered to best depict the yield distributions and satisfy the above mentioned conditions. (8) Terms A and B and S 0 of Equation ( 8) have been formulated as follows: The parameters a 0 , a 1 and δ are specific for each system and are used for model tuning.

Model Solution
Assuming that the model parameters including a 0 , a 1 , δ, k max , α are known, then the solution method for Equation ( 6) is as follows: For i=N or the component with maximum reactivity k N (k max ), the integral term of Equation ( 6) will be equal to zero because it is not generated during cracking reactions.Therefore, this equation is converted to the following form: Thus the concentration of the heaviest component at time t, C(kN,t), is obtained from its value at time t-δt by the following equation: For components i=1...N-1, the following implicit (backward Euler) finite difference equation was considered for numerical solution: (14) Considering an approximation, replacing concentration values at time t with the values at time t -δt in the integral term, the above equation is changed to Equation ( 15).This equation provides the possibility of calculating the concentration values of species i at time t from their values at time t -δt.Using small time steps δt in the calculation procedure could minimize this approximation error. (15) The calculation procedure starts from residence time t=0 where the component concentrations are known and equal to the feed.Using Equation (2) C versus θ data are converted to C versus k form.After a δt time step, the component concentrations at new time t are calculated from Equation ( 13) for i=N and Equation (15) for i=1...N-1.Since the functionality of C(K,t) versus K in Equation ( 15) is unknown, therefore the integral term has no analytical solution.For node i=N-1 numerical integration from the trapezoidal rule was considered because the integration range included two nodes i=N-1 and i=N.For i=1...N-2 the more accurate cubic spline method for integration was considered because the integration range included three nodes or more.Cubic spline uses a third order polynomial for estimation of the internal function of the integral term.
The calculation procedure continued up to the reactor ultimate residence time with δt step times.Finally, all component concentrations in the reactor outlet will be determined.Due to the high calculation volume, a code was carried out using MATLAB 7.10 simulation software.A personal computer equipped with Intel (R) Core (TM) 2 CPU (3.0 GHz) and 3.25 GB of RAM was used.
In the above procedure for product concentration prediction, it was assumed that the model parameters are known.But usually existing data are in the form of feed and product concentrations and model parameters are unknown.Therefore, another program was developed for parameter optimization, so that the predicted product concentration had minimum difference with experimental values.Thus, the computer program was developed in two modes: the prediction mode and the optimization mode.The developed program was validated with data reported by Elizalde et al. (2009).Then it was employed for the parameter estimation of the experimental data.
In the present work, a PSO algorithm has been applied for the parameter optimization.This method has many advantages including: fast convergence, finding the global optimum in the presence of many local optima, and simple programming (Shokri et al., 2014;Santos et al., 2015).Figure 1 shows a flowchart for the optimization mode of the program.

PSO method
The PSO method is a population based search optimization method that was originally proposed by Kennedy and Eberhart (1995).It is developed from swarm intelligence and is based on bird flocking and fish schooling to a promising position to achieve precise objectives in a multidimensional space.Particles fly around the parameter space with velocities which are dynamically adjusted according to their historical behavior and the particles have a tendency to fly towards the global solution over the course of the search process.In PSO, each potential solution of the problem is called a particle and the population of particles is called a swarm.
In the standard PSO, each particle i has a position X i and velocity V i that change according the following equality: (16) In this equality, m is the iteration number; W represents an inertial weight used as a trade-off between global and local exploration capabilities of the swarm.Small values of W result in more rapid convergence usually on a suboptimal position, while a too large value may prevent divergence.P besti and, G besti represent the position with the 'best' objective value found so far by particle i and the entire population respectively.C 1 and C 2 represent two positive acceleration coefficients regulating the speed of the particle's flying with respect to the best global and local positions respectively, R 1 (m) and R 2 (m) represent random numbers uniformly distributed in the closed interval of [0, 1].

,
, , , , ... This algorithm terminates when no significant improvement in the objective function is observed over a number of iterations.There are different ways in which improvement can be measured.For example, if the average change in particle positions is small, the swarm can be considered to have converged.Alternatively, if the average particle velocity over a number of iterations is approximately zero, only small position updates are made, and the search can be terminated.

Experimental Setup
In this study, a hydrocracking setup was utilized to obtain experimental data sets.The schematic diagram of the system is depicted in Figure 2. The core of the set up is the isothermal tubular reactor having an internal diameter of 2.0 cm and total length of 2.0 m.The reactor material was stainless steel.Thermocouples for temperature indication were inserted in the stainless steel thermowell at the center of reactor.Reaction temperature was supplied by an electrical heating jacket in the four zones.
A commercial NiMo supported on γ-alumina catalyst was used.In order to make the flow regime similar between the setup reactor and the industrial reactor, the ratio between the catalyst diameter and the reactor diameter should be approximately equal.For this purpose, the catalyst pellets were crushed to 0.8 mm.According to this Figure for the process, the feed and high-pressure hydrogen are mixed before entering the reactor section.In the operating conditions of the reactor, a trickle bed flow regime is governing.Inside the reactor, heavier hydrocarbons are converted to lighter liquid and gas hydrocarbons in the presence of hydrogen and catalyst at high pressure.Moreover, organic sulfur and nitrogen are transformed into H 2 S and NH 3 respectively.Reactor product is sent to the cooler.There it is cooled to ambient temperature.Thus, condensable hydrocarbon fractions are separated as liquid.
Then, a two-phase stream is sent into a highpressure separator (HPS) for separation into hydrogenrich gas and liquid hydrocarbon.Liquid hydrocarbon from the HPS is routed to the low pressure separator (LPS).The remaining dissolved gases are separated from the liquid in the LPS.Utilizing an Agilent gas chromatograph liquid products were analyzed by the simulated distillation (ASTM D2887) method.Analysis results are in form of true boiling point (TBP) versus distilled cumulative weight percent.The properties of the vacuum gasoil feed are presented in Table 1.
The experiments were performed at different space velocities of feed and different temperatures for thermal and catalytic cracking.In thermal cracking experiments, for hydrodynamic similarity 0.8 mm quartz pellets were loaded inside the reactor.
Table 2 shows the operating conditions of the experiments.The volume ratio of hydrogen to liquid hydrocarbon feed at the reactor inlet was approximately 1200:1.Moreover, the reactor pressure was maintained at 110 bars in all experiments.
Figures 3(a) and 3(b) give TBP versus cumulative weight percent for feed and product of catalytic and thermal cracking at 410ºC and 440ºC, typically.As can be seen in these Figures with increasing residence time or decreasing LHSV, the curves move down.
Cumulative weight percent (Cum wt%) versus θ are shown in Figures 4(a) and 4(b) for catalytic and thermal cracking at 410ºC and 440ºC, respectively.

RESULTS AND DISCUSSION
The continuous lumping model was applied to the results obtained from catalytic and thermal cracking experiments.Optimization of model parameters were performed for all operating conditions individually, using a PSO algorithm.Tables 3 and 4 show results of the model parameter optimization.As observed in these tables, the parameters are temperature dependent and the LHSV has no effect on them.
For evaluation of the model performance, in any temperature condition parameters optimized according to experimental results at LHSV=1(hr -1 ) were considered as the optimum set.Then by using this optimum set as input in the prediction mode of the computer program, product concentrations were calculated in other severities including LHSV=0.5 and LHSV=1.5.Typically, Figures 5(a) and 5(b) show predicted and experimental results for catalytic and thermal cracking at 430ºC and 470ºC, respectively.
In Figures 5(a) and 5(b), it is observed that the model gives acceptable product predictions at LHSV= 0.5 and 1.5, which indicates it suitable performance.Furthermore, the model results show good agreement with experimental results at LHSV=1.Consequently, the optimization algorithm has good performance.
In this study, the prediction performance of the model was evaluated based on two criteria: Root Mean Square error (RMSE) and Average Absolute Relative Error (AARE).These criteria are tabulated in Table 5.In this table y i shows the experimental data  6 for all experiments.
The results reveal that the continuous lumping model indicates good prediction for LHSV= 1.5 and 0.5, using optimized parameters at LHSV=1.7 and 8 compare product fraction yields (gasoline, diesel and residue) predicted by the model with those of the experimental data at different severities for catalytic and thermal cracking at 430ºC and 470ºC, respectively.The TBP range for gasoline is considered to be 30-150ºC, and 150-370ºC for diesel.The fraction of 370ºC + is classified as residue.Moreover, the normalized TBP (θ) range for these fractions is given in the tables.As observed, the model shows the best prediction for the diesel fraction.

CONCLUSIONS
A continuous lumping kinetic model was applied for the hydrocracking of vacuum gas oil feedstock in an experimental setup.The model parameters were optimized based on the experimental data for LHSV=1(hr -1 ) utilizing a PSO method.Then, the model was used to predict the product distribution for other severities.In order to evaluate the performance of the proposed model, a wide range of experimental data was taken from a hydrocracker setup.Some statistical

Figure 1 .
Figure 1.Flowchart for the optimization mode of the program.

Figure 5 .
Figure 5. Predicted and experimental results: (a) catalytic cracking at 430ºC; (b) thermal cracking at 470ºC.and ŷ i shows the model predicted data.A comparison of RMSE and AARE values for the model prediction is summarized in Table6for all experiments.The results reveal that the continuous lumping model indicates good prediction for LHSV= 1.5 and 0.5, using optimized parameters at LHSV=1.

Table 1 .
Characteristics of feed.

Table 2 .
Operating conditions of the experiments.

Table 3 .
Optimum values of parameters for catalytic cracking experiments.

Table 4 .
Optimum values of parameters for thermal cracking experiments.

Table 5 .
Model performance criteria.

Table 6 .
Comparison of the performance of continuous lumping model.

Table 7 .
Fraction yields prediction for catalytic cracking at 430°C.

Table 8 .
Fraction yields prediction for thermal cracking at 470°C.Exp=experimental criteria (AARE, RMSE) were used to evaluate the prediction performance of models.The comparison between the continuous model results and the experimental data shows that this model is capable of predicting the kinetic behavior well with an acceptable accuracy (maximum RMSE, 0.986 and maximum AARE, 0.041) for complex mixtures.Therefore, it was observed that the model can predict product yields with a highly acceptable degree of accuracy.