Energy & Fuels 2008, 22, 867–870
867
Logistic Regression Model for Isoconversional Kinetic Analysis of Cellulose Pyrolysis Junmeng Cai,* Ronghou Liu,* and Chen Sun Biomass Energy Engineering Research Centre, School of Agriculture and Biology, Shanghai Jiao Tong UniVersity, 800 Dongchuan Road, Shanghai 200240, P.R. China ReceiVed NoVember 6, 2007. ReVised Manuscript ReceiVed December 14, 2007
The possibility of applying the logistic regression model for fitting the kinetic conversion data of cellulose pyrolysis has been presented. It has been obtained that the pyrolysis conversion curves at different heating rates can be successfully described by the sixth degree polynomial logistic regression model. In this way, the generalized function can replace the original kinetic conversion curve, which is very useful to overcome the noise problems when analyzing the experimental data. On the basis of the analytical derivative of the known logistic regression model, the activation energies at the different values of the conversion degree have been obtained by means of the Friedman’s method.
1. Introduction The use of biomass as an energy source is of great importance, as it constitutes part of an alternative solution to the replacement of fossil fuel.1 The pyrolysis of biomass has received renewed attention because of the possibility of converting biomass into useful energetic products or into valuable chemicals.2 The information about biomass pyrolysis kinetics is necessary to accurately predict reactor behavior as well as to optimize the process of biomass pyrolysis. As cellulose is a main constituent of biomass, its pyrolysis kinetics has been the subject of extensive research.3–7 The thermogravimetric analysis (TGA) is one of the most commonly used techniques to study cellulose pyrolysis kinetics.8 Thermally stimulated solid-state reactions are commonly studied using a linear heating program. In nonisothermal TGA, the mass of the sample is continuously measured while the sample is subjected in a controlled atmosphere and the temperature is ramped at a constant heating rate. The differential isoconversional methods are usually used for performing the kinetic analysis of the experimental data by * To whom correspondence should be addressed. E-mail: jmcai@ sjtu.edu.cn (J.C.) or
[email protected] (R.L.). (1) Cai, J.; Liu, R.; Deng, C. An assessment of biomass resources availability in Shanghai: 2005 analysis. Renewable Sustainable Energy ReV. [Online early access]. DOI: 10.1016/j.rser.2007.04.003. Published Online: 2007. http://dx.doi.org/10.1016/j.rser.2007.04.003. (2) Bridgwater, A. V. Biomass fast pyrolysis. Therm. Sci. 2004, 8, 21– 49. (3) Völker, S.; Rieckmann, Th. Thermokinetic investigation of cellulose pyrolysis–impact of initial and final mass on kinetic results. J. Anal. Appl. Pyrolysis 2002, 62, 165–177. (4) Grønli, M.; Antal, M. J.; Várhegyi, G. A round-robin study of cellulose pyrolysis kinetics by thermogravimetry. Ind. Eng. Chem. Res. 1999, 38, 2238–2244. (5) Antal, M. J.; Várhegyi, G.; Jakab, E. Cellulose pyrolysis kinetics: revisited. Ind. Eng. Chem. Res. 1998, 37, 1267–1275. (6) Kim, S.; Eom, Y. Estimation of kinetic triplet of cellulose pyrolysis reaction from isothermal kinetic results. Korean J. Chem. Eng. 2006, 23, 409–414. (7) Capart, R.; Khezami, L.; Burnham, A. K. Assessment of various kinetic models for the pyrolysis of a microgranular cellulose. Thermochim. Acta 2004, 417, 79–89. (8) Antal, M. J.; Varhegyi, G. Cellulose pyrolysis kinetics: the current state of knowledge. Ind. Eng. Chem. Res. 1995, 34, 703–717.
thermogravimetry. To apply those differential isoconversional methods, one has to determine the values of the temperature and the conversion rate at a given degree of conversion. The values of the temperature at a given degree of conversion are usually estimated from the experimental data via a nonlinear interpolation. The corresponding values of the conversion rate are usually estimated by numerical differentiation of the experimental data, which are very sensitive to the noise in the data. Small errors in the experimental conversion data are magnified in the obtained differential conversion rates. To overcome the above problems, an alternative solution is the use of a fitting function for the experimental data. Several methods were proposed to do this research. Adma{evic´ et al.9 investigated the possibility of applying the normalized Weibull distribution function to fit the conversion data of nonisothermal dehydration of equilibrium poly(acrylic acid) hydrogel. Naya et al.10 proposed a method to find the parametric logistic regression model that fits the curves resulting from the decomposition of polymers during a TGA experiment. In the papers of Barbadillo et al.,11 Cao et al.,12 and Naya et al.,13 the method was presented for fitting an overall TGA curve of polymer degradation by the logistic mixture model, which could be used to separate the overlapping reactions involved in polymer degradation and to overcome the noise problems when analyzing the experimental data. Cai and Liu14 proposed a method to model the nonisothermal kinetic conversion data of (9) Adma{evic´, B.; Jankovic´, B.; Kolar-Anic´, L.; Minic´, D. Normalized Weibull distribution function for modelling the kinetics of non-isothermal dehydration of equilibrium swollen poly (acrylic acid) hydrogel. Chem. Eng. J. 2007, 130, 11–17. (10) Naya, S.; Cao, R.; Artiaga, R. Local polynomial estimation of TGA derivatives using logistic regression for pilot bandwidth selection. Thermochim. Acta 2003, 406, 177–183. (11) Barbadillo, F.; Fuentes, A.; Naya, S.; Cao, R.; Mier, J. L.; Artiaga, R. Evaluating the logistic mixture model on real and simulated TG cures. J. Therm. Anal. Calorim. 2007, 87, 223–227. (12) Cao, R.; Naya, S.; Artiaga, R.; García, A.; Varela, A. Logistic approach to polymer degradation in dynamic TGA. Polym. Degrad. Stab. 2004, 85, 667–674. (13) Naya, S.; Cao, R.; Ullibari, I. L.; Artiaga, R.; Barbadillo, F.; García, A. Logistic mixture model versus Arrhenius for kinetic study of material degradation by dynamic thermogravimetric analysis. J. Chemom. 2006, 20, 158–163.
10.1021/ef7006672 CCC: $40.75 2008 American Chemical Society Published on Web 02/07/2008
868 Energy & Fuels, Vol. 22, No. 2, 2008
Cai et al. Table 2. Values of R2 and ASE between the Experimental Data and the Logistic Regression Model Prediction
Figure 1. Overlay of the experimental data (dots) and the logistic regression model prediction (lines). Table 1. Values of T0 and Tf for Cellulose Pyrolysis Determined by TGA at Different Heating Rates β (K
min-1)
T0 (K)
Tf (K)
5 25 50
573.15 582.42 587.02
666.46 695.05 707.80
solid-state reactions by one or by the linear combination of a few Weibull distribution functions. The aim of this paper is to show the possibility of applying the nth-degree polynomial logistic regression model for fitting the kinetic conversion data of cellulose pyrolysis and to determine the activation energies dependent on the conversion degree based on the analytical data of the fitting. 2. Experimental Data The experimental data of cellulose pyrolysis kinetics has been obtained from the literature.15,16 The experiments were carried out on a thermogravimetric analyzer. The atmosphere used was nitrogen with a flow rate of 60 mL min-1. In all the experiments the material used was cellulose (Whatman No. 6 paper). Experiments in dynamic conditions were performed at heating rates of 5, 25, and 50 K min-1. The degree of conversion is expressed as R )
w0 - w w0 - wf
(1)
where w0, w, and wf refer to the initial, actual, and final mass fraction of the sample.
3. Mathematical Consideration 3.1. Logistic Regression Model. In this paper, the following expression of the logistic regression model17 has been proposed to fit the nonisothermal kinetic conversion data of cellulose pyrolysis: R(T) ) 1 - w
ePk(T - T0) 1 + ePk(T - T0)
(2)
where R is the conversion degree, T is the temperature, T0 is the initial temperature, w is the indeterminate constant, and Pk(T (14) Cai, J.; Liu, R. Weibull mixture model for modeling nonisothermal kinetics of thermally stimulated solid-state reactions: application to simulated and real kinetic conversion data. J. Phys. Chem. B 2007, 111, 10681–10686. (15) Conesa, J. A.; Caballero, J. A.; Marcilla, A.; Font, R. Analysis of different kinetic models in dynamic pyrolysis of cellulose. Thermochim. Acta 1995, 254, 175–192. (16) Conesa, J. A.; Caballero, J. A.; Reyes-Labarta, J. A. Artificial neural network for modeling thermal decompositions. J. Anal. Appl. Pyrolysis 2004, 71, 343–352. (17) Hosmer, D. W.; Lemeshow, S. Applied Logistic Regression; John Wiley & Sons Inc.: New York, 2000.
β (K min-1)
R2
ASE
5 25 50
0.9999769 0.9999547 0.9999453
1.13232 × 10-4 7.35586 × 10-5 7.92661 × 10-5
- T0) is the kth degree polynomial, Pk(T - T0) ) p0 + p1(T - T0) + p2(T - T0)2 + ... + pk(T - T0)k . The proposed model is a modified form of the logistic regression model presented in the paper of Naya et al.10 For fitting the kinetic conversion data to a logistic regression model, some estimation of the parameters in eq 2 is needed. We considered a nonlinear regression model: Re,j ) Rc,j + εj, j ) 1, 2, ... , nd
(3)
where Re,j is the experimental value of the kinetic conversion degree in the data point j, Rc,j is the value calculated from eq 2 in the corresponding point, nd is the number of data points, and εj are the errors, assumed to have normal distribution, with zero mean and constant variance. The parameters of the model have nd nd been estimated by minimizing ∑ j)1 εj2()∑ j)1 (Re,j - Rc,j)2). There are a number of algorithms for minimizing nd ∑ j)1 (Re,j - Rc,j)2 enabling us to find the best values of the parameters. In this study, the Levenberg–Marquardt method has been used for the calculation of the parameter values that nd minimize ∑ j)1 (Re,j - Rc,j)2 . Many methods use the gradient of the function to be minimized, while the Levenberg–Marquardt method uses a Jacobian instead. Detailed information about this algorithm can be found in the literature.18 For performing the Levenberg–Marquardt method, either general purposed mathematical software or a computer program developed in any programming language can be used. In this work, the Origin software developed by OriginLab Corporation has been employed. The software provides a data analysis and graphing workspace for scientists and engineers.19 Details about the implementation of the Levenberg–Marquardt algorithm in the Origin software can be found in its online help file.20 One of the problems that appears when using this fit is the choice of some starting points for the different parameters of the logistic regression model. Since this method is not easy and requires previous expertise, we propose the following method to overcome the problem. Supposing w ) w0, using the transformed form of eq 2, given by the relation
[
ln
]
1 - R(T) ) Pk(T - T0) w0 + R(T) - 1
(4)
we determine the coefficients of Pk(T - T0) by least-squares from the transformed data points: the horizontal coordinates do change (the temperature) while the new vertical coordinates are ln[(1 - R(T))/(w0 + R(T) - 1)] , where R(T) is the experimental value. These calculations can be easily performed by means of statistical software. 3.2. Friedman’s Method for Determination of the Activation Energy Dependent on the Conversion Degree. To obtain the dependence of the activation energy on the conversion (18) Watson, G. A. A Levenberg-Marquardt method for estimating polygonal regions. J. Comput. Appl. Math. 2007, 208, 331–340. (19) http://www.originlab.com/ (accessed Nov. 1, 2007). (20) http://www.originlab.com/pdfs/GettingStarted.PDF (accessed Nov. 1, 2007).
Kinetic Analysis of Cellulose Pyrolysis
Energy & Fuels, Vol. 22, No. 2, 2008 869
Figure 2. Differential conversion curves of cellulose pyrolysis at different heating rates, based on the logistic regression model. Table 3. Values of Tmax, rmax, and (dr/dT)max for Cellulose Pyrolysis at Different Heating Rates β (K min-1)
Tmax (K)
Rmax
(dR/dT)max (K-1)
5 25 50
635.14 660.35 670.45
0.63536 0.62964 0.59877
0.0255789 0.0242022 0.0227913
degree, some isoconversional methods are usually used.21–24 The differential isoconversional method suggested by Friedman25 is extensively used in the literature, as indicated by the more than 500 citations that we have found for the papers where the Friedman’s method has been reported (information on the number of citations has been obtained from the ISI Web of Science database). For nonisothermal conditions, when the temperature varies with time with a constant heating rate, the kinetics of a solidstate reaction are usually described by the following expression:26,27 β
dR ) Ae-E⁄RTf(R) dT
(5)
where β is the heating rate, A is the frequency factor, E is the activation energy, R is the universal gas constant, and f(R) is the differential conversion function associated with a certain reaction mechanism. Friedman’s method is based on eq 5 in logarithmic form and leads to ER,i
[ ( dRdT ) ] ) ln[A f(R) ] - RT
ln βi
R,i
R
i
(6)
R,i
where i and R designate a given value of the heating rate and the degree of conversion, respectively. For a given value of the degree of conversion, the plot of ln[βi(dR/dT)R,i] versus 1/TR,i should be a straight line whose slope can be used to evaluate the activation energy. Through Friedman’s method, we can get (21) Burnham, A. K. Computational aspects of kinetic analysis.: Part D: The ICTAC kinetics project–multi-thermal–history model-fitting methods and their relation to isoconversional methods. Thermochim. Acta 2000, 355, 165–170. (22) Simon, P. Isoconversional methods: Fundamentals, meaning and application. J. Therm. Anal. Calorim. 2004, 76, 123–132. (23) Pratap, A.; Lilly, T.; Rao, S. ; Lad, K. N.; Dhurandhar, H. D. Isoconversional vs. model fitting methods: A case study of crystallization kinetics of a Fe-based metallic glass. J. Therm. Anal. Calorim. 2007, 89, 399–405. (24) Ortega, A. Isoconversional method in CRTA. Thermochim. Acta 1997, 298, 161–164. (25) Friedman, H. L. Kinetics of thermal degradation of char-forming plastics from thermogravimetry. Application to a phenolic plastic. J. Polym. Sci., Part C: Polym. Symp. 1964, 1, 183–195. (26) Vyazovkin, S.; Dollimore, D. Linear and nonlinear procedures in isoconversional computations of the activation energy of nonisothermal reactions in solids. J. Chem. Inf. Comput. Sci. 1996, 36, 42–45. (27) Cai, J.; He, F.; Yi, W.; Yao, F. A new formula approximating the Arrhenius integral to perform the nonisothermal kinetics. Chem. Eng. J. 2006, 124, 15–18.
Figure 3. Dependence of ln[βi(dR/dT)R,i] vs (1000/(RTR,i)) at the different values of R, for cellulose pyrolysis (solid lines are linear fitting corresponding to different R).
Figure 4. Values of ER and ln[ARf(R)] at various values of R for cellulose pyrolysis.
different values of the activation energies for a process depending on the conversion degree. 4. Results and Discussion The experimental conversion data at different heating rates for cellulose pyrolysis are given in Figure 1. Values of the initial temperature (T0) and the final temperature (Tf) of the sample at various heating rates are presented in Table 1. The logistic regression model presented has been applied to the experimental kinetic conversion data of cellulose pyrolysis. The best fitting has been obtained with the sixth degree polynomial logistic regression model. The curves calculated by the logistic regression model are also plotted in Figure 1 for comparison. Table 2 lists the values of R2 (where R is the correlation coefficient) and ASE (average squared error) between the experimental data and the logistic regression model prediction. It can be observed that the logistic regression model matches the kinetic conversion data very well (for all heating rates, R2 is higher than 0.99994). Since the temperature dependence of the conversion degree is well described by the logistic regression model, the corresponding conversion rate can be given by the first derivative of eq 2 with respect to temperature. When this is done, it overcomes the noise problems when analyzing the experimental data. Figure 2 shows the differential conversion curves calculated by differentiating directly the logistic regression model. Table 3 lists the values of Tmax (the temperature at the maximum conversion rate), Rmax (the conversion degree at the maximum conversion rate), and (dR/dT)max (the maximum conversion rate) for cellulose pyrolysis at different heating rates.
870 Energy & Fuels, Vol. 22, No. 2, 2008
Since it has been shown that the experimental conversion data are well fitted by the sixth degree logistic regression model, the activation energies at different values of R (where R is taken from the mentioned fitting function) can be calculated. For this purpose, Friedman’s method has been applied. The isoconversion lines obtained by Friedman’s method at different conversion degree values are given in Figure 3, where the symbols designate the values of corresponding magnitudes on the left-hand and right-hand sides of the Friedman’s equation (eq 6), at different heating rates. Figure 4 shows the calculated dependencies of the ER and ln[ARf(R)] values on the degree of conversion for cellulose pyrolysis. From Figure 4, it can be seen that the activation energies change from 220 to 200 kJ mol-1 in the range of the conversion degree R ) 0.05-0.95. The activation energy values are similar to the 209.2 kJ mol-1 value presented in the paper of Conesa et al.16 Detailed knowledge about the phenomenon of the variations in the activation energy with the degree of conversion for cellulose pyrolysis is out of the scope of this work, and interested readers can refer to the paper of Vyazovkin.28
Cai et al.
5. Conclusions A method is proposed to find the sixth degree polynomial logistic regression model that fits the kinetic conversion data of the nonisothermal pyrolysis of cellulose. Since the fit is very good, the logistic regression model can replace the original conversion curves of cellulose pyrolysis. In this way, it is possible to obtain the conversion rate values by differentiating directly the conversion curve predicted by the logistic regression model. Making use of the analytical derivative of the fitting, the activation energies at the different values of the degree of conversion have been obtained by means of Friedman’s method, a classical differential isoconversion method. Acknowledgment. This research was supported by the National Natural Science Foundation of China, under Project No. 50776059. EF7006672
(28) Vyazovkin, S. On the phenomenon of variable activation energy for condensed phase reactions. New J. Chem. 2000, 24, 913–917.
