A comparative study of GA, PSO and SCE algorithms for estimating kinetics of biomass pyrolysis

Compared with fossil fuel energy, renewable energy offers many promising advantages, such as improved economic efficiency and reduced pollutant emission, and therefore has drawn much attention in a wide range of applications. Existing studies have shown that renewable energy sources will be the main supply in future development processes^[1−2]. Biomass is one of the most commonly used renewable energy resources, and comprehensive knowledge of pyrolysis reaction mechanisms and kinetics is extremely important to promote the development of biomass pyrolysis technologies. Determination of the unknown kinetic parameters for a specific type of biomass is crucial for accurately describing the pyrolysis process. By contrast to the time-consuming and difficult manual determination of parameters^[3], heuristic algorithms^[4−5] are more effective in solving the problem and are therefore widely used.

The most extensively used heuristic optimization algorithms include Genetic Algorithm (GA)^[6], Hill Climbing algorithm (HC)^[7], Particle Swarm Optimization (PSO)^[8], Shuffled Complex Evolution (SCE)^[9], etc. Li et al.^[10] extracted kinetics of MDF (medium density fiberboard) pyrolysis using GA. Abdelouahed et al.^[11] discussed the calculation method of kinetic parameters of biomass pyrolysis based on GA and concluded that the Kissinger method was the best method to determine kinetics. Ferreiro et al.^[6] studied the effect of biomass type on pyrolysis-related behaviors under different heating conditions using GA. Gong et al.^[7] developed a numerical model for pyrolysis of oriented strand board (OSB) using HC. Xu et al.^[12] extracted kinetics of lignocellulose, wood and rape straw using PSO. Aghbashlo et al.^[13] combined PSO with an adaptive network-based fuzzy inference system (ANFIS) to correlate the prediction of the kinetic constant of lignocellulose pyrolysis. Ding et al.^[14] compared various biomass models using SCE and concluded that a three-component parallel reaction scheme was more appropriate. Purnomo et al.^[15] compared the efficiency and accuracy of five different optimization algorithms in calculating biomass kinetics and showed that SCE was the most accurate. Different optimization algorithms have their own advantages and limitations and are useful in solving different types of problems. When applying these algorithms, rational selection according to the specific problems and appropriate parameter settings are crucial for obtaining reliable outcomes.

GA^[11] is an optimization algorithm based on Darwinian survival of the fittest theory. By simulating the natural evolutionary process, it transforms the problem to be solved into multiple individuals, and each individual encompasses a set of parameters that need to be solved. Through genetic, crossover and mutation operations, new individuals are created, evaluated and selected to gradually approach the optimal solution. GA^[6] has powerful parallel processing capability and is suitable for solving complex, nonlinear, multi-peaked and difficult problems. In addition, it's easy to implement with simple parameter settings. PSO^[8] transforms the encountered problem into a minimization problem by simulating the behavior of a flock of birds when searching for food. The algorithm first initializes random groups of particles that can be considered as possible solutions. By adjusting the position and speed of each particle, the swarm approaches the optimal solution. PSO has a global search capability with the influence of local optimal solutions. However, it should be noted that the search accuracy of PSO depends on the initial state and parameter setting of the population^[16], so reasonable initialization and parameter tuning are required. Additionally, the algorithm may converge prematurely, resulting in the search of only the local optimal solution but not the global optimal solution. SCE^[9,15] combines the concept of complexity science and the ideas of biological evolution. In each iteration, the algorithm updates the current optimal solution by randomly swapping parameters in multiple small systems contained in complex systems. SCE uses the idea of random swapping and can effectively deal with high-dimensional optimization problems with multiple local minima^[17]. However, the parameter setting of this algorithm is relatively complicated and requires reasonable parameter adjustment. In addition, the computational complexity of SCE is high, requiring a large amount of computational and storage resources.

Although all these three prevailing algorithms have been widely used in determining kinetics of biomass pyrolysis, few studies have compared their optimization capabilities, such as computation efficiency and accuracy. To challenge this issue and fill the research gap of existing studies, we launch a comparative study of these three algorithms by extracting kinetics of wood, a representative biomass material. Based on the originally obtained thermogravimetric experimental data, the contribution of each sub-reaction in the pyrolysis process is firstly identified by Gauss multi-peak fitting method, followed by the estimation of kinetics of each separated reaction using several model-free methods. Subsequently, the average values of the multiple sets of solutions are served as the initial solution to determine the search ranges of the three optimization algorithms. Then, the efficiency and accuracy of each algorithm are discussed in detail based on the optimization results.

Nitrogen atmosphere was selected to conduct the pyrolysis tests of beech wood using a NETZSCH STA449F3 thermal analyzer. Three different heating rates of 5, 10 and 20 K/min were used to heat the 5−7 mg wood powders from 290 to 1,070 K. Heating rate plays an important role in collecting proper experimental data. Very slow heating rate, such as 1 K/min, allows the reactions to come closer to equilibrium and there is less thermal lag in the sample and apparatus. Contrarily, high heating rates give faster experiments, which are more representative of the heating rates in fires, but deviate more from equilibrium and result in greater thermal lag. Larger heating rates are suitable for finding a wide range of decomposition, while smaller heating rates show better performance in the separation of individual events. Heating rates used in tests should preferably be in the range of 1−20 K/min as recommended by ICTAC (International Confederation for Thermal Analysis and Calorimetry) Kinetics Committee^[18]. Consequently, these three representative heating rates, 5, 10 and 20 K/min, were used in this study. A ceramic crucible was employed and a small hole was set in the center of the crucible lid to allow the release of volatiles during pyrolysis. All samples were dried in an oven for at least 72 h before tests to minimize the impact of moisture. Given that thermogravimetric experiments are highly reproducible, only three replicate experiments were performed at each heating rate to estimate the experimental uncertainty.

Model-free methods

Utilization of a suitable kinetic parameter calculation method allows a rough estimate of the search range for the subsequent optimization. Model-free methods, including the Kissinger method and isoconversional methods, are introduced in this section and will be utilized to calculate the kinetics of wood. Previous studies showed that the kinetic parameters calculated by the Kissinger method are very close to other isoconversional methods, such as KAS, Tang, DAEM methods.

In a single-step pyrolysis reaction, the degree of conversion of the solid, α, can be expressed as:

where $m$, ${m_0}$ and ${m_\infty }$ are the transient, initial and final masses of the sample, respectively. The reaction kinetic equation is:

where $t$ is time, T is the absolute temperature, $f\left( \alpha \right)$ is the differential form of the reaction model and $\lambda $ is the reaction rate constant. For a first order reaction, $\lambda $ can be expressed as:

where $A$, ${E_a}$ and $R$ refer to the pre-exponential factor, activation energy and ideal gas constant, respectively. With constant heating rate $\beta = dT/dt$ , Eq. (2) can be converted to:

The integral function of the conversion rate, $g\left( \alpha \right)$, can be expressed as:

where ${T_0}$ is the initial temperature.

Kissinger method

The Kissinger^[19] method uses the peak temperature of the reaction peak for plotting, and the equation can be expressed as:

where the subscript $n$ and $p$ refer to the n-th heating rate and peak of MLR curve. When the model is a first-order kinetic model, the value of $f'\left( {{\alpha _p}} \right)$ is −1, Eq. (6) can be further simplified as:

A straight line can be obtained by plotting $ln\left( {{\beta _n}/{T_{p,n}^2}} \right)$ versus $1/{T_{p,n}}$. The slope and intercept of the line can be employed to estimate ${E_a}$ and A, respectively.

KAS method

The formula for the KAS^[20] method is expressed as:

Using this formula, a straight line can be made by plotting $ln\left( {\beta /{T^2}} \right)$ to $1/T$ at different heating rates and arbitrary conversion rate. The slope of this line can be used to derive ${E_a}$. For first order reactions, the reaction mechanism functions are $f\left( \alpha \right) = 1 - \alpha $ and $g\left( \alpha \right) = - ln\left( {1 - \alpha } \right)$, and Eq. (8) becomes:

Tang method

Tang et al.^[21] proposed the Arrhenius temperature integral approximation for Eq. (5):

Substituting the approximation into Eq. (5) and taking logarithms on both sides, the expression of the Tang method is obtained:

At different heating rates and arbitrary conversion rate, plotting $ln\left( {\beta /{T^{1.894661}}} \right)$ versus $1/T$ makes a straight, and ${E_a}$ can be estimated based on the slope. Similarly, for a first order reaction, Eq. (13) can be transformed into:

The pre-exponential factor, $A$, can be calculated from the intercept.

DAEM method

DAEM^[22] is an effective method to study the reaction behavior of complex systems including an infinite number of parallel first order reactions. The simplified DAEM model can be expressed as:

Similarly, ${E_a}$ and $A$ can be calculated based on the slope and intercept of the linearly fitted line of $ln ( {\beta /{T^2}} )$ to $1/T$.

To simulate the measured mass and mass loss rate (MLR) of sample in the thermogravimetric analysis experiments, a 0D numerical model for the pyrolysis of a thermally thin solid is developed. The 0D model implies the temperature gradient inside the condensed phase is neglected compared to the traditional 1D heat transfer model. The general forms of pyrolysis reaction and reaction rate are:

where $Comp$ denotes component, $\theta $ is the stoichiometric coefficient by mass, $w$ and $n$ are reaction rate and the reactant concentration index, $c$ is the transient mass of the reactant after normalization, the subscript $j$ refers to the j-th reaction. The transient mass change rate of each component is calculated as:

where the subscript i denotes i-th component, ${N_j}$ is the total number of reactions. ${\theta _{j,i}}$ is positive or negative when the i-th component serves as a reactant or a product. The total transient mass (m) and MLR can be calculated as:

where ${N_s}$ and ${N_g}$ are the total numbers of solid and gaseous components, respectively. To commence simulation, the initial mass fraction of the starting reactants, the mass stoichiometry coefficients, and the three components of Arrhenius kinetics for each reaction need to be assigned. The initial values of m and MLR for the intermediate and final products are set to be zero.

Optimization algorithms

Three most frequently utilized optimization algorithms, GA, PSO, and SCE, are used to compare their performance by determining kinetics of wood combining the numerical model and experimental results. Detailed information of the three algorithms was introduced in our recent publications^[23,24]. A four-component reaction scheme is applied to describe the pyrolysis of wood, namely the evaporation of water, pyrolysis of cellulose, hemicellulose and lignin. Each reaction includes four unknown parameters, namely ${E_a}$, $A$, stoichiometric coefficient of solid product, and the reaction order. To gain comparison purpose, the initial search ranges, population sizes, iteration numbers, and objective functions of the three algorithms when implementing optimizations are set to be identical.

GA (Genetic algorithm)

GA simulates the process of population evolution, adopting a series of genetic operations such as selection, crossover and mutation for the current population to create a new generation and gradually progress the population to a state close to the optimal solution. In GA, each set of unknown parameters is referred to as an individual, and a combination of tens to thousands of individuals is defined as a population. Offspring, namely all potential solutions, are continuously produced by the overall population through genetic, crossover and mutation operations. An objective function, also known as fitness function, is essential for the assessment process. The objective function utilized in current study takes both m and MLR into account, as recommended by ICTAC Kinetics Committee^[25]:

where ${N_m}$, ${N_{MLR}}$ are the total numbers of experimental data points of m and MLR, the subscript exp and num denote experimental and numerical values, respectively, $\bar m$ and $\overline {MLR} $ are the average values of experimental m and MLR, respectively. A lower value of R² generally represents higher accuracy of the algorithm and better fit between experimental and simulation results. This objective function will also be used in the following PSO and SCE algorithms.

PSO (particle swarm optimization) algorithm

Inspired by the behaviors of bird populations, PSO is developed as an alternative optimization algorithm. PSO algorithm includes velocity and position models, where the velocities of particles are used to update the positions of the particles, and the positions represent the potential solutions in the search ranges. The velocity and position updating processes require each particle to keep in mind the previous optimal position as well as the global optimal position searched by all particles. At the beginning of iteration, the velocities and positions of the particles are randomly assigned according to a specified range. Then they are updated according to the following relationship:

where $i$ and $j$ are the numbers of particles and parameters in the particles, $ x $ and $ v $ represent the position and velocity of the particles, $ {x^{pb}} $, $ {x^{gb}} $ are the local and global best positions of the particles, q₁ and q₂ are two random numbers located within [0,2].

SCE (shuffled complex evolution) algorithm

SCE is an efficient optimization algorithm due to its excellent global search performance and convergence speed, and it is suitable for solving high-dimensional complex nonlinear problems. The main principle of SCE is that each parameter has its own specific search range, within which the fitting value is calculated for each randomly generated individual as well as the ranking. The probability that an individual is selected is:

where ${x_k}$ represents the k-th individual ranked from the lowest to the highest individual fitness value, $n$ is the total number of individuals, $f\left( {{x_k}} \right)$ is the function determining the assignment of individuals, $p\left( {{x_k}} \right)$ is the probability that the k-th individual is selected. Individuals are divided into multiple complexes for evolution based on sorting. During evolution, the better adapted particles would replace the less adapted particles as parents to generate the next generation. Then, all groups are mixed and reordered and the process is repeated until the convergence condition is satisfied.

In this section, we analyze the kinetics of wood pyrolysis based on thermogravimetric experimental data obtained at the three different heating rates (5, 10 and 20 K/min), followed by an analytical discussion of the performance of the three algorithms, GA, PSO and SCE, in terms of both efficiency and accuracy in optimizing the kinetic parameters.

Analyses of thermogravimetric results

Figure 1 shows the measured mass and MLR curves during wood pyrolysis at different heating rates. Both mass and MLR data are normalized by the initial mass ${m_0}$. As expected, higher heating rates shift the mass and MLR curves toward higher temperature range. In Fig. 1b, it can be reasonably inferred that wood pyrolysis encompasses four main decomposition reactions: water evaporation, decomposition of hemicellulose, cellulose, and lignin, corresponding to the first minor peak before 400 K, the asymmetrical left shoulder of the main peak, the main peak, and the long tail, respectively. In order to evaluate the contribution of each reaction to the total MLR and estimate the relevant kinetics, Gauss multi-peak fitting method^[26] is adopted to resolve the MLR curve at each heating rate into multiple elemental curves, as shown in Fig. 2. Each elemental curve corresponds to an elemental pyrolysis reaction. Gauss multi-peak fitting method needs only one parameter when separating overlapped multiple peaks, which renders its simplicity in shape and easy optimization compared with other separation methods, such Weibull, Gauss, Gamma, and Symmetric logistic, as demonstrated by ICTAC Kinetics Committee^[25].

Figure 1. 

Dependencies of (a) mass and (b) MLR of beech wood on temperature at different heating rates.

Figure 2. 

Separated MLR curves using Gauss multi-peak fitting method at different heating rates.

Kinetic analyses

Pyrolysis reaction scheme of beech wood, including four elemental reactions, is listed in Table 1. Based on previous studies^[24], the reaction of water, hemicellulose and cellulose can be described by first order reactions, while lignin is a high-order reaction. Consequently, there are totally 12 parameters, excluding the stoichiometric coefficient of water evaporation, in each optimization run. Therefore, in this section the kinetic analyses are implemented for these four components using the methods introduced previously.

Based on the separated curves in Fig. 2, the Kissinger method is first used to estimated kinetics of each reaction. Plotting ${\text{ln}}( {\beta /{T_p^2}} )$ ~ $1/{T_p}$ and executing linear fitting, shown in Fig. 3, $A$ and ${E_a}$ of each reaction can be estimated based on the slop and intercept, as listed in Table 2. The relatively good linearity of the fitted lines confirms the reliability of the Kissinger method. In addition, $A$ and ${E_a}$ are also estimated by the KAS method (${\text{ln}} ( {\beta /{T^2}} )$ ~ $1/T$), Tang method ($\ln ( {\beta /{T^{1.894661}}} )$ ~ $1/T$), and DAEM method (${\text{ln}}( {\beta /{T^2}} )$ ~ $1/T$). The derived A and ${E_a}$ are listed in Table 2, and the linear fittings at different conversion rates are depicted in Figs 4−6, respectively.

Figure 3. 

Linear fittings of ${\text{ln}}( {\beta /{T_p}^2} )$ vs $1/{T_p}$ in the Kissinger method.

Figure 4. 

Linear fittings of ${\text{ln}}( {\beta /{T^2}} )$ vs $1/T$ plots in the KAS method: (a) water, (b) hemicellulose, (c) cellulose, (d) lignin.

Figure 5. 

Linear fittings of $\ln ( {\beta /{T^{1.894661}}} )$ vs $1/T$ plots in the Tang method: (a) water, (b) hemicellulose, (c) cellulose, (d) lignin.

Figure 6. 

Linear fittings of ${\text{ln}}( {\beta /{T^2}} )$ vs $1/T$ plots in the DAEM method: (a) water, (b) hemicellulose, (c) cellulose, (d) lignin.

The calculated A and ${E_a}$ by the four analytical methods in Table 2 are similar despite some minor deviations. Figure 7 shows the variation trends of calculated ${E_a}$ of the four main components in wood with varying conversion rates. Since in the Kissinger and DAEM methods ${\text{ln}}( {\beta /{T^2}} )$ linearly depends on $1/T$ which is similar to those in the KAS and Tang methods, only the results of the KAS and Tang methods are plotted in Fig. 7. In Fig. 7a & c, ${E_a}$ of water and cellulose decline linearly with increasing α, whereas in Fig. 7b & d opposite variation trend is observed for hemicellulose and lignin. Vyazovkin et al.^[25] suggested a pyrolysis process can be described by a single step reaction only if the difference between the maximum and the minimum values of ${E_a}$ is smaller than 20% of the average value. In Fig. 7a−d, the ratios of $\left( {{E_{a,max}} - {E_{a,min}}} \right)/{E_{a,ave}}$ are 18.6%, 4.01%, 3.38%, 4.15%, implying all these reactions could be regarded as single step reactions.

Figure 7. 

Dependencies of calculated ${E_a}$ on α using the KAS and Tang methods: (a) water, (b) hemicellulose, (c) cellulose, (d) lignin.

Comparison of accuracy of GA, PSO and SCE

To compare the optimization performance of the three focused algorithms, the same optimization settings are used for each algorithm. The average values of the four analytical methods in Table 2 are used in determining the search ranges. More specifically, these average values are employed as the mean values of the search ranges during optimization. The lower bounds of search ranges are set to be 0.1 times of the mean values, while the upper bounds are selected to ensure these mean values are the average values of the search ranges. Figure 8 shows the values of objective function (R²) and computation times (t_com) of the optimization runs using the three algorithms with 200-3000 population sizes and fixed iteration number of 200. R² of GA changes irregularly as the population size is smaller than 1000, but it declines with further increase of population size. This phenomenon indicates the accuracy of GA strongly depends on population size. R² of PSO descents quickly for population size smaller than 800 and increase slightly beyond this range. While the R² of SCE are always very low and changes slightly with varying population size, suggesting the accuracy of SCE is very high and is barely affected by population size. In Fig. 8b, t_com of the three algorithms all increase with population size. The difference is that both t_com and its increasing rate of SCE are much larger than those of GA and PSO, implying the computation efficiencies of GA and PSO are approximately identical to each other and both are higher than that of SCE. As introduced in by Shi et al.^[24], each iteration of SCE involves multiple complex systems, and therefore the computational complexity is higher than the other two algorithms, requiring a larger amount of computational and storage resources. Meanwhile, Table 3 lists t_com and R² of the three algorithms with varying population sizes. Distinctly, the computation efficiencies of the three algorithms are ranked as GA ≈ PSO > SCE, while the accuracies are ranked as SCE > PSO > GA.

Figure 8. 

Objective function values and the computation times of the three algorithms when optimizing kinetics of wood pyrolysis with 200-3000 population sizes and 200 iterations.

Comparison of efficiency of GA, PSO and SCE

Similarly, the convergence efficiencies of the three algorithms are compared using fixed population size of 3000 and 6000 iterations, and the evolutions of R² are portrayed in Fig. 9. GA, PSO, and SCE converge at 1000, 500, and 800 iterations, respectively. Meanwhile, the decreasing rate of R² of PSO is higher than the two others. Consequently, it can be reasonably concluded that the convergence efficiencies of three algorithms can be ranked as PSO > SCE > GA. However, in Fig. 8b, the computation efficiencies are ranked as GA ≈ PSO > SCE. In Fig. 8b and Table 3, SCE consumes much more time than GA and PSO, up to approximately 12 times. Apparently, each algorithm has its inherent merits and limits. Overall, PSO is a more favorable algorithm featuring high accuracy, computation efficiency and convergence efficiency. Nevertheless, particular care should be taken when applying it since it may fall into local optimal solution.

Figure 9. 

Objective function value evolutions of the three algorithms when optimizing kinetics of wood pyrolysis with 3,000 population size and 6,000 iterations.

Parameter validation

In all the optimization runs in above, the minimum values of R² for GA, PSO, and SCE are 0.0826, 0.0747, and 0.0743, respectively. These optimizations are conducted at 10 K/min heating rate and the corresponding results are listed in Table 4. Apparently, with identical optimization setting, only minor deviations exit among different algorithms. In order verify the reliability of the optimized parameters, the experimental results of 5 and 20 K/min heating rates are predicted using the numerical model and the kinetics listed in Table 4, as exhibited in Figs 10−12 corresponding to GA, PSO, and SCE, respectively. Existing studies^[23] revealed that hemicellulose pyrolysis is mainly responsible for the asymmetric shoulder of the MLR curve, cellulose decomposition is related to the main peak, and lignin decomposition is located at a higher temperature range. In Figs 10 & 11, corresponding to GA and PSO algorithms, the numerical curves fit the experimental data well and the relative locations of the subpeaks agree with the literature. Nevertheless, in Fig. 12 the peak temperature of lignin is lower than that of cellulose despite the overall good agreement, conflicting with the existing conclusion. This divergence is presumably caused by the compensation effect among A, ${E_a}$ and reaction order when extracting kinetics from TGA data using optimization algorithms, which is still an unsolved problem^[25].

Figure 10. 

Comparison between experimental and numerical MLRs using optimized parameters of GA at 5, 10 and 20 K/min heating rates.

Figure 11. 

Comparison between experimental and numerical MLRs using optimized parameters of PSO at 5, 10 and 20 K/min heating rates.

Figure 12. 

Comparison between experimental and numerical MLRs using optimized parameters of SCE at 5, 10 and 20 K/min heating rates.

Pyrolysis of beech wood was investigated experimentally at three heating rates of 5, 10, and 20 K/min. Based on the measured MLR curves, the overlapping peaks of wood was first separated by Gauss multi-peak fitting method to identify their contributions. Then, four analytical methods were used to determine ${E_a}$ and A of each reaction. Subsequently, the accuracy, computation efficiency, and convergence efficiency of GA, PSO and SCE algorithms were compared at 10 K/min heating rate. It was found that in terms of optimization accuracy, SCE was the best followed PSO, and then GA. While for computation efficiency, PSO was the best, then GA and SCE. Whereas considering convergence efficiency, the three algorithms were ranked as PSO > SCE > GA. All these indicated each algorithm had its inherent advantages and limits, and PSO featured better overall performance. Furthermore, the reliability of the optimized kinetics was verified by predicting the remaining experimental data at 5 and 10 K/min which were not used during parametrization.

An interesting conclusion attained is that no any single algorithm excels others in all aspects. A potential solution to this issue may be developing more advanced hybrid algorithms which could better balance the accuracy, computation efficiency, convergence efficiency, storage resource, etc. Meanwhile, these heuristic optimization algorithms can also be coupled with some artificial intelligence (AI) algorithms, such as machine learning, deep learning, support vector machines, decision trees, random forest, and metaheuristics. Even though AI has been successfully applied in many engineering applications, few attempts invoking AI have been made to challenge the complex pyrolysis process of biomass. All these need in-depth exploration in future studies.