Engineering Applications

## **Chapter 4**

## Programming an Evolutionary Algorithm for the Estimation of Non-Linear Damping Vibration Parameters

*Carlos A. Lara and Cesar Guerra*

## **Abstract**

The use of genetic algorithms (GAs) has branched out into various disciplines such as mechanical engineering, providing solutions in cases where some models do not have a mathematical solution. In the field of mechanical vibrations, there are empirical nonlinear models that seek to represent the physical behavior of certain elements, such as in aeronautical applications, where stiffness and damping in structures can present hysteresis and can be represented by means of the Bouc-Wen (BW) model. This model includes constants that define the behavior of non-linear stiffness and damping, which are difficult to obtain since they are empirical models. This work presents the results of programming a GA to estimate the BW model constants for wire rope springs, commonly used as vibration isolators that have nonlinear stiffness and damping resulting in hysteresis behavior.

**Keywords:** genetic algorithms, Bouc-Wen model, non-linear stiffness, damping, rope springs

### **1. Introduction**

The study of vibration control has taken great importance and necessity in recent years. The emergence of new isolating and vibration-dissipating elements has prompted new studies about their behavior in hysteresis damping. Damping is difficult to model this is often due to more than one phenomenon, for example, a combination of viscous damping, internal damping, dry friction, viscoelastic effects, etc. [1]. There are various of insulator configurations, among which are: metal springs, i.e., helical or leaf springs, and viscoelastic elements such as rubber, neoprene, silicon, air springs, etc. Steel cable springs are characterized by their high energy storage and dissipation capacity, based on dry friction [2].

In the study of hysteresis damping, there is a problem in describing its non-linear behavior [1]. One solution is the model of BW [3], from which heuristic algorithm techniques have been proposed, which seek to estimate the parameters, factors, error reduction, among others, of which it is possible to establish an estimated

model of the system, and from it, design control and/or estimation strategies [4]. Given the non-linear nature of the model, it has been approached by different methods, including the following: Gauss–Newton, modified Gauss–Newton, Least squares, Simplex, Levenberg–Marquardt, extended Kalman filter, reduced gradient methods, genetic algorithms (GAs), real-coded GAs, Differential Evolution, adaptive laws, etc.

It is well known that the classical linearized analysis of the dynamical systems can lead to results that are reasonably accurate only when the minimum (rest position) force and the displacements are of such magnitude that the relative change in force during the motion is small. In practice, however, very often some or all of these assumptions are violated, so that in many dynamical systems nonlinear phenomena may completely alter intuitively expected behavior and can drastically change their dynamical responses [5].

In mechanical vibrations systems, the nonlinear phenomenon can be presented principally in the springs elements and/or the dampers models, significant results have also been obtained to represent these phenomena.

Recently, the use of new insulator mechanical in several systems, for example in Aeronautics, has prompted research to design new non-linear model representatives of this elements, where a memory-dependent, multivalued relation between force and deformation, i.e., hysteresis, is often observed in structural materials and elements, such as reinforced concrete, steel, base isolators, dampers, and soil profiles.

Many mathematical models have been developed to efficiently describe such behavior for use in time history and random vibration analyses. One of the most popular is the BW class of hysteresis which is used to describe the hysteretic behaviors of structures in nonlinear dynamic and stochastic analyses.

BW model is used to describe the non-linear behavior of the stiffness and damping of an element, where the restoring force becomes highly nonlinear showing significant hysteresis. The hereditary nature of this nonlinear restoring force indicates that the force cannot be described as a function of the instantaneous displacement and velocity. Accordingly, many hysteretic restoring force models were developed to include the time-dependent nature using a set of differential equations.

BW model is a semi-empirical model that contains several parameters and is one of the most used hysteretic models, and it was introduced by Robert Bouc [6] and extended by Yi-Kwei Wen [7] who demonstrated its versatility by producing a variety of hysteretic patterns.

Being a semi-empirical model, the BW model contains semi-empirical parameters which should be esteemed using several mathematical and empirical strategies, such as Gauss–Newton [8], modified Gauss–Newton [9], Least squares [10], Simplex [11], Levenberg–Marquardt [12], extended Kalman filters [13] among others.

Recently, the use of GAs for the estimation of BW parameters has been used, for example, Kwok et al. [14] used a GAs to estimate the parameters of the BW model with characteristics of non-symmetrical hysteresis; Wang et al. [15] used a novel differential evolution algorithm for estimation of parameters of asymmetric hysteresis loops.

Meanwhile, Charalampakis et al. [16] presented a new identification method that determines the parameters of Bouc-Web hysteresis based on a hybrid evolutionary algorithm which utilizes selected stochastic operators.

In most cases, the problem lies in that the BW model can be easily solved because it combines an algebraic equation with a differential equation, in addition, it is found that there are redundant parameters [17]. One solution to deal with

*Programming an Evolutionary Algorithm for the Estimation of Non-Linear Damping Vibration… DOI: http://dx.doi.org/10.5772/intechopen.113070*

this problem, users of the BW model often fix some parameters to arbitrary values, while other users eliminate the redundant parameters via a process of normalization.

In this work, a discrete approximation of the BW model is proposed to facilitate the estimation of BW parameters using an efficient evolutionary algorithm called "Evonorm". The programming of algorithm helped model the behavior of physically loaded/unloaded springs in an experimental setup.

### **2. BW discrete model approximated**

Starting from a mass-spring-damper vibratory system, where the interest is based in accordance with the type of damping present (hysteresis or a combination of phenomena) as shown in **Figure 1**, the application of BW for the study of its behavior in terms of the damping present in the system.

Mass (m) represents the inertia in *kg* and Damper (c) is the viscous damping in *N\*s/m*; both elements are considered as lineal elements of the system. Now, the spring (k) represents the stiffness lineal in N/m, but in our study, this element contains both the non-linear stiffness and non-linear damping.

Actuality, no any mathematical representation exists for this non-lineal phenomena, but an empirical representation known as BW model and are described in the following.

The BW model was proposed initially by Bouc early in 1967 and subsequently generalized by Wen in 1976. Its typical equations are expressed as follows:

$$R = k\_e \varkappa + k\_h z \tag{1}$$

$$\dot{\boldsymbol{z}} = \mathbf{A}\dot{\boldsymbol{x}} - \beta |\dot{\boldsymbol{x}}||\boldsymbol{z}|^{n-1}\boldsymbol{z} - \gamma \dot{\boldsymbol{x}}|\boldsymbol{z}|^{n} \tag{2}$$

Where *R* is the restoring force, *x=x*(*t*) is the deformation of spring, *z = z(t)* is known as hysteresis displacement and not physically measured, *x*\_ ¼ *dx=dt* and *z*\_ ¼ *dz=dt*; *ke* ¼ *αk* and *kh* ¼ ð Þ 1 � *α k*; where *k* is the lineal stiffness coefficient.

The α parameter is the ratio of post-yield to pre-yield stiffness, *A, n, β, γ* are parameters that control the hysteresis shape.

Analyzing Eqs. (1) and (2), it is found that there are redundant parameters in BW model. That is why to estimate the control parameters in the BW model, have been defined with alternative models or modified; for example: "The normalized BW model" presented in [18], "A multi-objective optimization algorithm for Bouc–Wen– Baber–Noori model [19].

**Figure 1.** *m-k-c system.*

In this work, the first step in order to solve the BW equation consists of eliminating the derivative function using a numerical approximation, in this case discrete approximation of the first-order of the original model called "Euler-backward discretization" is used as:

$$
\dot{\omega} \approx \frac{\varkappa\_{k+1} - \varkappa\_k}{\Delta t} \tag{3}
$$

Where *xk*þ1, *xk* are the two-sample data of the x(t) and Δt is the time sample of the sample data width as shown in **Figure 2**.

Now, from Eq. (1), it is possible to define the differential equation *<sup>R</sup>*\_ <sup>¼</sup> *kex*\_ <sup>þ</sup> *khz*\_, therefore applying the discrete approximation (3) to *R*\_ , *x*\_ and *z*\_ the following equation is applied.

$$R\_{k+1} - R\_k = k\_e(\varkappa\_{k+1} - \varkappa\_k) + k\_h(z\_{k+1} - z\_k) \tag{4}$$

$$z\_{k+1} - z\_k = A(\boldsymbol{\omega}\_{k+1} - \boldsymbol{\omega}\_k) - \beta |\boldsymbol{\omega}\_{k+1} - \boldsymbol{\omega}\_k| |\boldsymbol{z}\_k|^{n-1} z\_k - \gamma (\boldsymbol{\omega}\_{k+1} - \boldsymbol{\omega}\_k) |\boldsymbol{z}\_k|^n \tag{5}$$

This model can be solved in programming if the deformation Δ*x* ¼ *xk*þ<sup>1</sup> � *xk* is fixed or Δ*R* ¼ *Rk*þ<sup>1</sup> � *Rk* is fixed, is clear that if Δ*x* is fixed, so Δ*R* is easy to determine. In the sample, *k=0* is clear that *x*<sup>0</sup> ¼ 0, *R*<sup>0</sup> ¼ 0 and *z*<sup>0</sup> ¼ 0, now the next values can be calculated using the next Pseudocode as shown in **Table 1**.


#### **Table 1.**

*Pseudocode to solve the BW discrete model.*

**Figure 2.** *Euler backward approximation.*

*Programming an Evolutionary Algorithm for the Estimation of Non-Linear Damping Vibration… DOI: http://dx.doi.org/10.5772/intechopen.113070*

### **3. Evolutionary algorithm model**

Evonorm [20] is an evolutionary algorithm used in this work to estimate the parameters of BW model. Next definitions are required to understand the algorithm.

#### **3.1 Variables definition**

*Decision variables (Y)*. These are the *m-*variables of the system with unknown values. The determination of these values will be the target of the algorithm.

$$Y \equiv \{ \mathcal{y}\_1, \mathcal{y}\_2, \mathcal{y}\_3, \dots, \mathcal{y}\_m \}\tag{6}$$

*Design variables (X)*. These are the *n*-variables of the system with known values, these values set that will allow determinate the decision values (*Y*) using the algorithm.

$$X \equiv \{ \mathbf{x}\_1, \mathbf{x}\_2, \mathbf{x}\_3, \dots, \mathbf{x}\_n \} \tag{7}$$

Decision values (*Y*~). These are the values of Decision variables (*Y*).

$$
\tilde{Y} = \begin{bmatrix} \tilde{y}\_1 \ \tilde{y}\_2 \ \cdots \ \tilde{y}\_m \end{bmatrix} \tag{8}
$$

Design values (*X*~). These are the values of Design variables that will be used to determine the decision values in the algorithm. The efficiency of the algorithm depends on selecting the *<sup>p</sup>*-samples required, therefore *<sup>X</sup>*<sup>~</sup> is a matrix of *<sup>p</sup>* � *<sup>n</sup>* dimension

$$
\bar{X} = \begin{bmatrix}
\bar{X}(\mathbf{1}) \\
\bar{X}(\mathbf{2}) \\
\vdots \\
\bar{X}(p)
\end{bmatrix} = \begin{bmatrix}
\ddot{\mathbf{x}}\_1(\mathbf{1}) & \ddot{\mathbf{x}}\_2(\mathbf{1}) & \cdots & \ddot{\mathbf{x}}\_n(\mathbf{1}) \\
\ddot{\mathbf{x}}\_1(\mathbf{2}) & \ddot{\mathbf{x}}\_2(\mathbf{2}) & \cdots & \ddot{\mathbf{x}}\_n(\mathbf{2}) \\
\vdots & \vdots & \ddots & \vdots \\
\ddot{\mathbf{x}}\_1(p) & \ddot{\mathbf{x}}\_2(p) & \cdots & \ddot{\mathbf{x}}\_n(p)
\end{bmatrix} \tag{9}
$$

*Fitness function.* It is a function or heuristic algorithm required to evaluate if the values of decision variables and design variables are correct, namely.

*Minimal error (e)*. It is a selected value such that the fitness function evaluation, allows to affirm that the *Y*~ values are correct:

$$\sum\_{i=1}^{p} |f\left(\tilde{X}(i), \tilde{Y}\right)| < \overline{e} \tag{10}$$

*Population (P)*. It is a *q*-data set (individuals) of the Decision values candidates to solve (10), therefore *P* is a matrix of *q* � *m* dimension

$$P = \begin{bmatrix} \check{Y}(\mathbf{1}) \\ \check{Y}(\mathbf{2}) \\ \vdots \\ \check{Y}(q) \end{bmatrix} = \begin{bmatrix} \check{y}\_1(\mathbf{1}) & \check{y}\_2(\mathbf{1}) & \cdots & \check{y}\_m(\mathbf{1}) \\ \check{y}\_1(\mathbf{2}) & \check{y}\_2(\mathbf{2}) & \cdots & \check{y}\_m(\mathbf{2}) \\ \vdots & \vdots & \ddots & \vdots \\ \check{y}\_1(q) & \check{y}\_2(q) & \cdots & \check{y}\_m(q) \end{bmatrix} \tag{11}$$

#### **3.2 Selection (Ts)**

Since each row of the *P*-matrix in (9) represents a possible solution of the condition (10), it is necessary determinate the value of the error for each row of the population matrix (11),

$$E(j) = \sum\_{i=1}^{p} |f\left(\tilde{X}(i), \tilde{Y}(j)\right)| \tag{12}$$

Now, for each j-individual in the matrix population, it is necessary to evaluate the contribution it makes to the solution (12), for which the rows of the P-matrix (11), must be ordered in ascending order. Now, must be selected the fits *Ts*-individuals (rows) to mutate and crossover to generate the new Population.

#### **3.3 Mutation and crossover**

In order to avoid the algorithm being trapped in local optima, Evonorm uses random variables with normal distribution. The normal distribution function is a random variable and describes many random phenomena that occur in everyday life. It simulated the normal distribution function with two parameters, the first is the mean and it is a numeric measure of the central tendency of the random variable. The second parameter is the standard deviation, and it is a measure of the dispersion of a variable around the mean. A normal distribution function can be used to represent a set of possible values of a decision variable, so it is necessary to use a set of parameters (mean and standard deviation) of the normal distribution function per decision variable (18).

Therefore, the mutation is generated using each *k*-column of the (11) and to *Ts*-individuals select as follows:

$$\mu\_k = \frac{\sum\_{i=1}^{T\mathfrak{s}} \tilde{\mathbf{y}}\_k(i)}{T\mathfrak{s}} \sigma\_k = \sqrt{\frac{\sum\_{i=1}^{T\mathfrak{s}} \left(\tilde{\mathbf{y}}\_k(i) - \mu\_k\right)^2}{T\mathfrak{s}}} \tag{13}$$

At the same time, the new population is generated as follows:

$$\check{\mathcal{Y}}\_{k}(i) = \begin{cases} \mu\_{k} + N(0, \mathbf{1})\sigma\_{k} \text{ if } U(.) > \mathbf{0.5} \\ \check{\mathcal{Y}}\_{k}(i) + N(0, \mathbf{1})\sigma\_{k} \text{ if } U(.) \le \mathbf{0.5} \end{cases} \tag{14}$$

Where *<sup>N</sup>*ð Þ¼ 0, 1 <sup>P</sup><sup>12</sup> *<sup>i</sup>*¼<sup>1</sup>*U*ðÞ�*:* 6, and *<sup>U</sup>*ð Þ*:* is a random value between 0 and 1. The pseudocode of the Evonorm algorithm is shown in **Table 2**.


*Programming an Evolutionary Algorithm for the Estimation of Non-Linear Damping Vibration… DOI: http://dx.doi.org/10.5772/intechopen.113070*

## **4. Experiment result**

## **4.1 Experimental data**

The experimental data were obtained from tests made to a set of four wire rope isolators (WRI) in parallel (like the one shown in **Figure 3**), to have stability during the test.

The mechanical tests were carried out in the universal machine Shimadzu (**Figure 4**) of 10 KN, with controlled displacement for compression-decompression load.

The values of Force-deformation in the wire ropes were plotted in **Figure 5**, showing the hysteretic behavior of the WRI.

*Picture and physical characteristics of the wire rope isolator (WRI) used in the experimental procedure.*

**Figure 4.** *Experimental setup for the static monotonic and cyclic tests.*

**Figure 5.** *Data result from experimentation.*

From the data that are plotted in **Figure 5**, representative samples were taken, this is to input variables for the evolutionary algorithm, which performs the calculation of distribution function parameters to generate new individuals.

#### **4.2 Algorithm programming**

In order to apply the Evonorm algorithm to estimate the BW parameters model, which is necessary to define the decision variables and design variables of the system. In the BW model the data set of hysteresis loop are the design variables, Meanwhile the parameters *α*, *β*, *γ*, *n*, *A*, *k* are the decision variables:

$$Y \equiv \{a, \beta, \gamma, n, A, k\} \tag{15}$$

$$X \equiv \{ \mathfrak{x}\_{k+1}, \mathfrak{x}\_k, R\_{k+1}, R\_k \} \tag{16}$$

The values of the design variables can be obtained from the result experiment with the values in **Figure 5**. In this case, the values that appear below are selected::


Now, it is necessary to determine the fitness function *F Xd* ð Þ , *Yd* in the algorithm, this function included the design variables *Xd* and de decision variables *Yd*. Therefore, the fitness function is the approximated model discrete (6).

$$F(X,Y) = R\_{k+1} - R\_k - k\_t \Delta\_X k\_h \left( A \Delta\_X + \beta |\Delta\_X| |\mathbf{z}\_k|^{n-1} \mathbf{z}\_k \gamma + \Delta\_X |\mathbf{z}\_k|^n \right) \tag{18}$$

The major objective will be to minimize the function *F Xd* ð Þffi , *Yd* 0.

Next values was used in the algorithm programming: Number of individuals: p = 50, number of selected individuals: *Ts = 25* and the numbers of iterations: (iters) *Nr* = 100.

On the other hand, it is necessary establish limits in the Design variables, this provides to algorithm's heuristic find the optimal values and in minor iterations.

The limits values of the Design variables are:

$$a\epsilon[0,1], \beta\epsilon[0.1,0.9], \gamma\epsilon[0.1,0.9],\tag{19}$$

$$A\epsilon[1,10], n\epsilon(1,2], k\epsilon[10,1000] \tag{20}$$

The selection of these values and the number of samples is important to make the algorithm run more efficient, so they must be strategic and minority. Strategic because, for this case, the middle of the loop was selected, taking as samples: the ends, the point at the intersection with the vertical, as well as two intermediate points at the ends of the intersection with the vertical; and minorities so as not to increase the computational cost of the algorithm. The result obtained is shown in **Table 3** for different runs.

*Programming an Evolutionary Algorithm for the Estimation of Non-Linear Damping Vibration… DOI: http://dx.doi.org/10.5772/intechopen.113070*


#### **Table 3.**

*Evonorm algorithm results.*

#### **Figure 6.** *Error graph.*

**Table 3** shows the results of the Evonorm algorithm with the values that define the shape of the hysteresis loop.

The convergence of the error showing the efficiency of the algorithm is also shown. **Figure 6** shows the graph of the error percentage against the number of iterations performed by the algorithm.

**Figure 8.** *Graph of actual data vs estimated data*

On the other hand, the real results obtained from experimental tests (represented by the red graph) and the results with the BW parameters obtained from the evolutionary algorithm (represented by the blue graph) were graphed (**Figure 7**). The data that was chosen to feed the BW model were those that were obtained with the minimum percentage of error. **Figure 8** shows the comparison of the real values versus those produced by the evolutionary algorithm.

### **5. Conclusion and future work**

The evolutionary algorithm can obtain adjusted BW parameters that can be fed into the model. The experimental hysteresis loops were fitted, with the parameter from Evonorm, in only 4 runs of the evolutionary algorithm.

The convergence of the minimum value of the error of the different permutations was achieved in 25 iterations, which is acceptable and can be improved by feeding more amount of experimental data to the evolutionary algorithm.

Of the 4 runs that were carried out with the Evonorm evolutionary algorithm, very similar ranges of errors were obtained. With the best option, very similar graphs, of the real experimentation and of the algorithm, were also obtained, and for that reason can be confident with the values that the algorithm delivers.

From **Figure 7**, the parameters provided by the EVONORM evolutionary algorithm, which are fed to the BW model, generate a graph that is very close to the graph that is made with the values obtained from experimental tests. Given the above, performing the selection of the input variables to the evolutionary algorithm correctly and strategically makes the estimated output values that are fed to the model guaranteed precision to the hysteresis loop compared to that produced by experimental tests.

It is clear that the correct selection of limits values in Evonorm is a condition to ease the convergence, these values can be obtained using the knowledge of the expert; *Programming an Evolutionary Algorithm for the Estimation of Non-Linear Damping Vibration… DOI: http://dx.doi.org/10.5772/intechopen.113070*

but is possible to use a neural network, so that after the training, estimate the values of the limits, which a future work.

Finally, concerning the error value, in **Figure 6** this remains constant after 25 iterations and the ideality is that this value declines gradually, some changes in the limits values can tackle this, furthermore, it is possible to improve tunning the percentage values in Eq. (14), for example, if *U*ð Þ*:* >0*:*4, if *U*ð Þ*:* ≤0*:*4.

## **Author details**

Carlos A. Lara\* and Cesar Guerra Universidad Autónoma de Nuevo León, San Nicolás de los Garza, México

\*Address all correspondence to: carlos.laraoc@gmail.com

© 2023 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

## **References**

[1] Rao SS. Mechanical Vibrations. 6th ed. Pearson; 2016. p. 1152

[2] Ledezma D, Tapia P, Ferguson N, Brennan M, Tang B. Recent advances in shock vibration isolation: An overview and future possibilities. ASME. Applied Mechanics Reviews. 2019;**71**(6):060802. DOI: 10.1115/1.4044190

[3] de Silva CW. Vibration Damping, Control, and Design. Press; 2019. p. 634

[4] Hassani V, Tjahjowidodo T, Do TN. A survey on hysteresis modeling, identification and control. Mechanical Systems and Signal Processing. 2014;**49** (1–2):209-233. DOI: 10.1016/J.YMSSP. 2014.04.012

[5] Warminski J, Lenci S, Carttmell P, Rega G, Wiercigroch M. Nonlinear Dynamic Phenomena, Mechanical Engineering Systems in Solids Mechanics and Its Applications. 2012. p. 181

[6] Bouc R. Forced vibration of mechanical systems with hysteresis. In: Proceedings of the Fourth Conference on Nonlinear Oscillation. Prague, Czechoslovakia; 1967. p. 315

[7] Wen YK. Method for random vibration of hysteretic systems. Journal of Engineering Mechanics. American Society of Civil Engineers. 1976;**102**(2): 249-263

[8] Yar M, Hammond J. Parameter estimation for hysteretic systems. Journal of Sound and Vibration. 1987; **117**(1):161-172. DOI: 10.1016/0022-460X (87)90442-1

[9] Kunnath S, Mander J, Fang L. Parameter identification for degrading and pinched hysteretic structural concrete systems. Engineering

Structures. 1997;**19**(3):224-232. DOI: 10.1016/S0141-0296(96)00058-2

[10] Sues R, Mau S, Wen Y. System identification of degrading hysteretic restoring forces. Journal of Engineering Mechanics. 1988;**114**(5):833-846. DOI: 10.1061/(ASCE)0733-9399(1988) 114:5(833)

[11] Zhang H, Foliente G, Yang Y, Ma F. Parameter identification of inelastic structures under dynamic loads. Earthquake Engineering and Structural Dynamics. 2002;**31**:1113-1130. DOI: 10.1002/eqe.151

[12] Ni Y, Ko J, Wong C. Identification of non-linear hysteretic isolators from periodic vibration tests. Journal of Sound and Vibration. 1998;**217**(4):737-756. DOI: 10.1006/jsvi.1998.1804

[13] Lin J, Zhang Y. Nonlinear structural identification using extended Kalman filter. Computers and Structures. 1994; **52**(4):757-764. DOI: 10.1016/0045-7949 (94)90357-3

[14] Kwok N, Ha Q. Bouc-Wen model parameter identification for a MR fluid damper using computationally efficient FA. ISA Transaction. 2007;**46**:167-179. DOI: 10.1016/j.isatra.2006.08.005

[15] Wang G, Chen G, Bai F. Modeling and identification of asymmetric Bouc-Wen hysteresis for piezoelectric actuator via a novel differential evolution algorithm. Sensor and Actuators. 2015; **235**:105-118. DOI: 10.1016/j. sna.2015.09.043

[16] Charalampakis A, Koumousis V. Identification of Bluc-Wen hysteresis system by a hybrid evolutionary algorithm. Journal of Sound and

*Programming an Evolutionary Algorithm for the Estimation of Non-Linear Damping Vibration… DOI: http://dx.doi.org/10.5772/intechopen.113070*

Vibration. 2008;**314**:571-185. DOI: 10.1016/j.jsv.2008.01.018

[17] Zhu X, Lu X. Parametric identification of Bouc-Wen model and its application in mild steel damper modeling. Procedia Engineering. 2011; **14**:318-324. DOI: 10.1016/j.proeng. 2011.07.039

[18] Zhu H, Rui X, Yang F, Zhu W, Wei M. An efficient parameters identification method of normalized Bouc-Wen model for MR damper. Journal of Sound and Vibration. 2019; **448**:146-158. DOI: 10.1016/j.jsv.2019. 02.019

[19] Li Z, Noori M, Zhao Y, Wan C, Feng D, Altabey WA. A multi-objective optimization algorithm for Bouc–Wen– Baber–Noori model to identify reinforced concrete columns failing in different modes. Proceedings of the Institution of Mechanical Engineers, Part L: Journal of Materials: Design and Applications. 2021;**235**(9):2165-2182. DOI: 10.1177/14644207211020028

[20] Torres L. Evonorm: Easy and effective implementation of estimation of distribution algorithms. Journal of Research in Computing, Science. 2008; **23**:75-83

## **Chapter 5**

## The Genetic Algorithm and its Application in Calculating the Kinetic Parameters of the Thermoluminescence Curve

*Nguyen Duy Sang*

## **Abstract**

This chapter explores the use of genetic algorithms as a tool for calculating the kinetic parameters of the thermoluminescence curve. Genetic algorithm is a search algorithm inspired by the process of natural selection, and it has proven to be effective in solving optimization problems in various fields. Author used genetic algorithm to estimate the activation energy and frequency factor of the thermoluminescence curve, which are important parameters in determining the dosimetric properties of materials. The results showed that genetic algorithm could accurately estimate the kinetic parameters of the thermoluminescence curve with high precision and efficiency compared to conventional methods. This approach can also handle noisy data and reduce the impact of outliers on the estimation process. Furthermore, author demonstrated that genetic algorithm can be generalized to different types of the thermoluminescence curves, such as those generated by different materials or under different experimental conditions. The proposed method is fast, accurate, and robust, making it useful for researchers in the field of dosimetry who require precise estimations of these parameters.

**Keywords:** genetic algorithms, kinetic parameters, thermoluminescence, materials, dosimetry

### **1. Introduction**

The study of thermoluminescence (TL) has been an important area of research for many years due to its numerous applications in dating, dosimetry, and material characterization. The TL curve is a graph of the light emitted by a material as it is heated, and it provides valuable information about the energy states of the material. In order to extract useful information from the TL curve, it is necessary to analyze its kinetic parameters, such as activation energy and frequency factor. One powerful method for calculating these parameters is the genetic algorithm, which is a

computational technique based on the principles of natural selection and evolution. Genetic algorithm is based on the simulation of genetic processes in living organisms and the principle of natural evolution. The experimentally obtained TL spectra are curves that are fitted according to different kinetic models [1–3]. The algorithm works with a biological population, each of which represents the ability to adapt to the explore space through synchronous combinations of evolutionary processes such as selection, crossover, and mutation [4].

Python provides a powerful platform for TL analysis, allowing for complex data analysis and visualization for accurate interpretation of experimental results. Python can be used for TL analysis by performing data analysis on experimental data obtained from TL measurements.

Here are some steps that can be followed: Import necessary modules: The first step is to import the necessary modules that will be required for data analysis. Some commonly used modules include NumPy, Pandas, Matplotlib, and Seaborn. Load experimental data: Load the experimental data into Python using the Pandas module. The data should be stored in a CSV or TXT file format that can be easily loaded into Python. Data processing: The next step is to clean and process the data. This involves removing any noise or unwanted signals present in the data. This can be done using NumPy's signal processing functions such as filtering and smoothing. Plotting data: Visualizing the processed data is an important step in TL analysis. Using Matplotlib and Seaborn, create different types of plots such as scatter plots, line plots, histograms. Curve fitting: The next step is to fit a curve to the data. This helps in determining the kinetic parameters of the material being studied. Python provides various curve fitting algorithms such as least squares fitting and nonlinear regression. Analyzing results: Finally, interpret the results obtained from curve fitting and apply them to the material being studied.

### **2. Genetic algorithm approach to the TL curves**

#### **2.1 Theoretical basis**

Genetic algorithm (GA) is based on the simulation of genetic processes in living organisms and the principle of natural evolution. The algorithm works with a biological population, each of which represents the ability to adapt to the explore space through synchronous combinations of evolutionary processes such as selection, crossover, and mutation [4]. The GA tuning process includes the following steps: (i) Setting initial chromosome and its encoding; (ii) Calculate GOK model distributions for each variable from individuals of a population and evaluate the fit of the fitness function; (iii) Select randomly parents and go to the next step; (iv) Crossover and mutation, go to the next step; (v) Select randomly individuals and go to the next step; (vi) Accept the results if there is better fitness value than the worst explore in the population and go to the next step, reject the worst explore and return step (iii); (vii) If the number of pre-determined steps (stopping condition) is reached and go to the next step; and (viii) Print results for best explore in population and GA finish. The block diagram of the GA is shown in **Figure 1**.

Explanation of **Figure 1**: **Initial population**: create an initial population of candidate solutions randomly. Each solution represents a potential solution to the problem being solved. Selection: The selection process is typically based on a **fitness function** that evaluates each individual's performance. **Crossover:** Two selected individuals are *The Genetic Algorithm and its Application in Calculating the Kinetic Parameters… DOI: http://dx.doi.org/10.5772/intechopen.112198*

**Figure 1.** *Flowchart of the genetic algorithm.*

combined to produce offspring with characteristics from both parents. The crossover point is chosen randomly, and the resulting offspring replace their parents in the population. **Mutation**: To maintain diversity in the population, some individuals undergo random mutation, which introduces new characteristics not present in the parent population. **Evaluation**: The fitness function is used to evaluate the performance of each individual in the population, including the newly generated offspring. **Termination**: The algorithm terminates when either a termination criterion is met (e.g., a maximum number of generations) or the best solution has been found.

```
Code in Python:
```
The Python file FSG\_GA.py at GitHub [5].

#### **2.2 The thermoluminescence kinetic equation**

The most commonly used models for analyzing such signals are the first order kinetics (FOK), second order kinetics (SOK), the general order kinetics (GOK) [6, 7].

The empirical equation describing the first-order TL curve has the form:

$$\mathbf{I}(\mathbf{T}) = \mathbf{I}\_{\mathbf{m}} \exp\left[\mathbf{1} + \frac{\mathbf{E}}{\mathbf{k}\mathbf{T}} \frac{\mathbf{T} - \mathbf{T}\_{\mathbf{m}}}{\mathbf{T}\_{\mathbf{m}}} - \frac{\mathbf{T}^2}{\mathbf{T}\_{\mathbf{m}}^2} \exp\left(\frac{\mathbf{E}}{\mathbf{k}\mathbf{T}} \frac{\mathbf{T} - \mathbf{T}\_{\mathbf{m}}}{\mathbf{T}\_{\mathbf{m}}}\right) \left(\mathbf{1} - \frac{2\mathbf{k}\mathbf{T}}{\mathbf{E}}\right) - \frac{2\mathbf{k}\mathbf{T}\_{\mathbf{m}}}{\mathbf{E}}\right] \tag{1}$$

Code in Python:

Graph of first-order kinetic TL curve:

#### **Figure 2.**

*First-order glow peaks of TL curve.*

Explanation of **Figure 2**: Graph of the first-order kinetic model of the TL spectrum with a peak and symmetry.

The empirical equation describing the second-order TL curve has the form:

$$\mathbf{I}(\mathbf{T}) = 4\mathbf{I}\_{\mathrm{m}} \exp\left(\frac{\mathbf{E}}{\mathrm{k}\mathbf{T}} \,\frac{\mathbf{T} - \mathbf{T}\_{\mathrm{m}}}{\mathbf{T}\_{\mathrm{m}}}\right) \left[\frac{\mathbf{T}^{2}}{\mathbf{T}\_{\mathrm{m}}^{2}} \left(\mathbf{1} - \frac{2\mathbf{k}\mathbf{T}}{\mathbf{E}}\right) \exp\left(\frac{\mathbf{E}}{\mathbf{k}\mathbf{T}} \,\frac{\mathbf{T} - \mathbf{T}\_{\mathrm{m}}}{\mathbf{T}\_{\mathrm{m}}}\right) + \mathbf{1} + \frac{2\mathbf{k}\mathbf{T}\_{\mathrm{m}}}{\mathbf{E}}\right]^{-2} \tag{2}$$

Code in Python:

*The Genetic Algorithm and its Application in Calculating the Kinetic Parameters… DOI: http://dx.doi.org/10.5772/intechopen.112198*

Graph of second-order kinetic TL curve:

**Figure 3.** *Second-order glow peaks of TL curve.*

Explanation of **Figure 3**: A graph of the second-order kinetic model of the TL spectrum with a peak and an asymmetrical shape, sloped forward of the peak.

General-order glow peaks are produced in intermediate cases (neither of firstorder nor of second-order). The four parameters describing a glow peak are Im, E, Tm, and b.

The empirical equation describing the general-order TL curve has the form:

$$\mathbf{J(T)} = \mathbf{I\_m} \mathbf{b}^{\frac{\mathbf{b}}{\mathbf{b}-1}} \exp\left(\frac{\mathbf{E}}{\mathbf{kT}} \frac{\mathbf{T} - \mathbf{T\_m}}{\mathbf{T\_m}}\right) \left[ (\mathbf{b} - \mathbf{1}) \left(\mathbf{1} - \frac{2\mathbf{k}\mathbf{T}}{\mathbf{E}}\right) \frac{\mathbf{T}^2}{\mathbf{T\_m}^2} \exp\left(\frac{\mathbf{E}}{\mathbf{k}\mathbf{T}} \frac{\mathbf{T} - \mathbf{T\_m}}{\mathbf{T\_m}}\right) + \mathbf{1} \quad \text{(3)}\right]$$

$$+ \frac{2\mathbf{k}\mathbf{T\_m}(\mathbf{b} - \mathbf{1})}{\mathbf{E}}$$

Code in Python:

Graph of general-order kinetic TL curve:

**Figure 4.** *General-order glow peaks of TL curve.*

Explanation of **Figure 4**: The graph of the general-order kinetic model of the TL spectrum has the form of an intermediate peak of a first- and second-order kinetics model, with a slightly sloping front peak.

In addition to the three kinetic models above, the TL spectrum is also matched according to the mixing model and the one-trap one-center recombination model. Each model has its own advantages and disadvantages in calculating TL trap parameters.

Details and full Python source code for TL kinetic models can be found on GitHub [5].

#### **2.3 Fitting TL curves to estimate the energy value**

### *2.3.1 Straight line fitting*

The initial rise (IR) method introduced by Garlick and Gibson [8] is used to estimate the trapping energy value of the TL curve. The IR method works as follows: A sample is irradiated with a known dose of ionizing radiation. The sample is then heated at a constant rate, and the emitted light is measured as a function of temperature. This is called a TL curve. The TL curve is divided into equal temperature intervals, and the total integrated light output for each interval is calculated. The integrated light output is plotted against the natural logarithm of the heating rate for each interval. The slope of the resulting straight line is proportional to the activation energy required to release the trapped charges in the sample. This method is based on the principle that the intensity of TL increases

*The Genetic Algorithm and its Application in Calculating the Kinetic Parameters… DOI: http://dx.doi.org/10.5772/intechopen.112198*

initially with the temperature. The activation energy can be calculated using the Arrhenius Eq. (4)

$$\mathbf{I}\approx\exp\left(-\frac{\mathbf{E}}{\mathbf{k}\mathbf{T}}\right)\tag{4}$$

where E is the activation energy, k is the Boltzmann constant, and slope is the slope of the straight line. The IR method is repeated at different doses of ionizing radiation, and the activation energy is plotted against the dose. The trapping energy value can be estimated from the intercept of this plot with the x-axis (dose axis). In summary, the IR method involves measuring the total integrated light output as a function of temperature for a sample irradiated with a known dose of ionizing radiation [8–11]. An example of an IR region of a glow peak is shown in **Figure 5**.

Explanation of **Figure 5**: The low-temperature peak tail in this region increases up to a critical temperature TC which is less than Tm. The values of E from the IR remain true for some critical values of temperature up to TC, corresponding to a TL intensity IC smaller than about 10% of the maximum TL intensity Im [12].

Code in Python:

**Figure 5.** *Diagram of initial rise (IR).*

Linear spectral fitting graph:

**Figure 6.** *Fitting line-fit of the TL1 to estimate E.*

Explanation of **Figure 6**: On the left is the TL1 sample matched by GA and applying the IR method, the kinetic parameters are also calculated, E = 0.72 eV.

#### *2.3.2 Gaussian peak spectral fitting*

The peak shape (PS) method is generally called as Chen's [13] method, which is used to determine the kinetic parameters of the main glow peak of the TL materials. This

*The Genetic Algorithm and its Application in Calculating the Kinetic Parameters… DOI: http://dx.doi.org/10.5772/intechopen.112198*

**Figure 7.** *Diagram of peak shape (PS).*

method is mainly based on the temperatures Tm, T1, and T2, which are the peak temperatures, the temperatures at half of the maximum intensity, on the ascending and descending parts of the peak, respectively. Calculation of the activation energy by PS method is shown in **Figure 7**. The expression deduced by Chen [13] and valid for any kinetics is given by Eq. (9), where α stands for τ, δ, and ω, in which the low-temperature half width τ = Tm - T1, the high-temperature half width δ = T2 -Tm, and the total half intensity width ω = T2 - T1. The values of c<sup>α</sup> and b<sup>α</sup> are summarized as defined in Eq. (6).

$$\mathbf{E}\_{\alpha} = \mathbf{c}\_{\alpha} \left( \frac{\mathbf{k} \mathbf{T}\_{\mathbf{m}}^{2}}{\alpha} \right) - \mathbf{b}\_{\alpha} (2 \mathbf{k} \mathbf{T}\_{\mathbf{m}}) \tag{5}$$

$$\begin{cases} \mathbf{c}\_{\mathbf{t}} = \mathbf{1}, \mathbf{51} + \mathbf{3}, \mathbf{0} \left( \boldsymbol{\mu}\_{\mathbf{g}} - \mathbf{0}, 42 \right) \\\\ \mathbf{b}\_{\mathbf{t}} = \mathbf{1}, \mathbf{51} + 4, 2 \left( \boldsymbol{\mu}\_{\mathbf{g}} - \mathbf{0}, 42 \right) \\\\ \mathbf{c}\_{\boldsymbol{\delta}} = \mathbf{0}, \mathbf{97} \boldsymbol{\delta} + \mathbf{7}, \mathbf{3} \left( \boldsymbol{\mu}\_{\mathbf{g}} - \mathbf{0}, 42 \right) \\\\ \mathbf{b}\_{\boldsymbol{\delta}} = \mathbf{0} \\\\ \mathbf{c}\_{\alpha} = 2, \mathbf{52} + \mathbf{10}, 2 \left( \boldsymbol{\mu}\_{\mathbf{g}} - \mathbf{0}, 42 \right) \\\\ \mathbf{b}\_{\alpha} = \mathbf{1} \end{cases} \tag{6}$$

where μ<sup>g</sup> is the so-called geometrical shape or symmetry factor that determines the order of the kinetics. The order of the kinetics depends on the glow PS. The value of μ<sup>g</sup> for first- and second-order kinetics is 0.42 and 0.52, respectively. The symmetry factor in GOK model can be evaluated from the following Eq. (7). The TL glow peaks corresponding to second-order kinetics are characterized by an almost symmetrical shape, whereas first-order peaks are asymmetrical [6].

$$
\mu\_{\rm g} = \frac{\delta}{\alpha} = (\mathbf{T\_2} - \mathbf{T\_m})/(\mathbf{T\_2} - \mathbf{T\_1})\tag{7}
$$

#### **Figure 8.** *The TL spectrum of the R4 sample is matched with four peaks.*

Explanation of **Figure 7**: describes how to fit the Gaussian peak spectrum and the PS method to calculate the trap energy. The calculation results depend on the T1, T2, and T values of the TL peak.

#### Code in Python:

TL curves from R package tgcd [14] with 22 TL curves are measured from different materials provided by George Kitis. Among them is the TL curve denoted as R4, which is used to fit the Gaussian peak shape and calculate the kinematic parameters. The results are shown in **Figure 8**.

Details and full Python source code for calculating the activation energy can be found on GitHub [5].

Explanation of **Figure 8**: spectrum of R4 consisting of four peaks matched and peaked by genetic algorithm. The kinetic parameters of the spectrum are calculated reasonably, and the obtained FOM coefficient is very small.

#### **2.4 Calculation of the frequency factor**

The intensity I(t) of the TL signal is measured at time t after the start of the experiment. This magnitude TL I(t) is proportional to the derivative �dn/dt,

*The Genetic Algorithm and its Application in Calculating the Kinetic Parameters… DOI: http://dx.doi.org/10.5772/intechopen.112198*

#### **Figure 9.**

*The spectrum of K2GdF5:Tb curve and results of calculating E and s values.*

depending on the measurement conditions. In experimental research, experimentalists pay much attention to the frequency factor of the TL signal. Kitis et al. [7] obtain the following analytic equation for s with the GOK model:

$$\mathbf{s} = \left(\frac{\beta \mathbf{E}}{\mathbf{k} \mathbf{T}\_{\mathrm{m}}}\right) \exp\left(\frac{\mathbf{E}}{\mathbf{k} \mathbf{T}\_{\mathrm{m}}}\right) \left[\mathbf{1} + (\mathbf{b} - \mathbf{1}) \frac{2 \mathbf{k} \mathbf{T}\_{\mathrm{m}}}{\mathbf{E}}\right]^{-1} \tag{8}$$

Thus, kinematic parameters such as E and s of the TL curve will be estimated according to Eqs. (5) and (8). Each peak coordinate of the TL curve including two parameters of temperature and TL intensity will be recorded after each mouse click on the screen containing the TL curve. The kinematic order of the GOK model is also selected and changed from 1 to 2 until the FOM values of the TL curve reach the minimum condition.

Code in Python:

K2GdF5 materials doped with concentrations of 10%Tb are widely used in radiation dosimetry and materials science [15]. The spectrum of K2GdF5:Tb curve and results of calculating E and s values are shown in **Figure 9**.

Explanation of **Figure 9**: depicts the result of spectral matching of sample K2GdF5: Tb with four peaks. The calculation results of sample K2GdF5:Tb in terms of E and s values are also calculated with high accuracy.

Details and full Python source code for calculating the frequency factor can be found on GitHub [5].

### **3. Conclusion**

The application of GA to TL curve analysis provides an efficient and effective method for the estimation of kinetic parameters such as activation energy, frequency factor, and order of reaction. The algorithm can explore a large search space of candidate solutions and converge toward a solution that optimizes the fitness function. The use of GA in TL curve analysis has significant implications for the field of geochronology and archeological dating, as it provides a powerful tool for the accurate and precise determination of the age of materials.

## **Author details**

Nguyen Duy Sang Can Tho University, Can Tho, Vietnam

\*Address all correspondence to: ndsang@ctu.edu.vn

© 2023 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

*The Genetic Algorithm and its Application in Calculating the Kinetic Parameters… DOI: http://dx.doi.org/10.5772/intechopen.112198*

## **References**

[1] Sang ND, Hung NV, Hung TV, Hien NQ. Using the computerized glow curve deconvolution method and the R package tgcd to determination of thermoluminescence kinetic parameters of chilli powder samples by GOK model and OTOR one. Journal of Nuclear Instruments Methods B. 2017;**394**:113-120

[2] Sang ND. Estimate half-life of thermoluminescence signals according to the different models by using python. Journal of Taibah University for Science. 2021;**15**(1):599-608

[3] Sang ND, Thi HHQ. Using the genetic algorithm to detect kinetic parameters of thermoluminescence glow curves. Nuclear Instruments and Methods in Physics Research Section B: Beam Interactions with Materials and Atoms. 2022;**517**:33-42

[4] Coley DA. An Introduction to Genetic Algorithms for Scientists and Engineers. World Scientific; 1999. p. 244

[5] Sang ND. Github. 2023. Available from: https://github.com/sangduynguye n/ga\_book

[6] Pagonis V, Kitis G, Furetta C. Numerical and Practical Exercises in Thermoluminescence. United States of America: Springer; 2006

[7] Kitis G, Gomez-Ros JM, Tuyn JWN. Thermoluminescence glow-curve deconvolution functions for first, second and general orders of kinetics. Journal of Physics D: Applied Physics. 1998;**31**(19): 2636-2641

[8] Garlick GFJ, Gibson AF. The electron trap mechanism of luminescence in Sulphide and silicate phosphors. Proceedings of the Physical Society. 1948;**60**:574-590

[9] Sang ND. Study of the effect of gamma-irradiation on the activation energy value from the thermoluminescence glow curve. Journal of Taibah University for Science. 2017; **11**(6):1221-1221

[10] Rawat NS et al. Use of initial rise method to analyze a general-order kinetic thermoluminescence glow curve. Nuclear Instruments and Methods in Physics Research Section B: Beam Interactions with Materials and Atoms. 2009;**267**(20):3475-3479

[11] Correcher V, Garcia-Guinea J. Potential use of the activation energy value calculated from the thermoluminescence glow curves to detect irradiated food. Journal of Radioanalytical and Nuclear Chemistry. 2013;**298**(2):821-825

[12] McKeever SWS. Thermoluminescence of Solids. Cambridge: Cambridge University Press; 1988

[13] Chen R. On the calculation of activation energies and frequency factors from glow curves. Journal of Applied Physics. 1969;**40**:570

[14] Peng J, Dong Z, Han F. Tgcd: An R package for analyzing thermoluminescence glow curves. SoftwareX. 2016;**5**:112-120

[15] Van Tuyen H et al. Spectroscopic studies of K2GdF5:Nd3+ single crystals for incredibly strong NIR emission at 864 nm. Journal of Physics and Chemistry of Solids. 2022;**161**:110454

## *Edited by Yann-Henri Chemin*

In this edition of *Genetic Algorithms - Theory, Design and Programming*, we present a series of scientific contributions that delve into the intricate theoretical foundations and practical nuances of genetic algorithms (GAs). Beyond the academic realm, GAs have demonstrated profound applications in societal decision-making and engineering optimization, showcased through real-world examples and case studies. A dedicated section on programming principles offers a thorough guide for implementing GAs across diverse languages. This edition, tailored for researchers and academics, serves as a testament to the scientific advancements within the field, inviting readers to explore the nuanced journey from theoretical constructs to pragmatic applications in the dynamic landscape of GAs.

*Andries Engelbrecht, Artificial Intelligence Series Editor*

Published in London, UK © 2024 IntechOpen © your\_photo / iStock

Genetic Algorithms - Theory, Design and Programming

IntechOpen Series

Artificial Intelligence, Volume 23

Genetic Algorithms

Theory, Design and Programming

*Edited by Yann-Henri Chemin*