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

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.

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 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.

Figure 1.

m-k-c 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}x+k_{h}z$$E1

$$\dot{z}=A\dot{x}-\beta \left|\dot{x}\right|{\left|z\right|}^{n-1}z-\gamma \dot{x}{\left|z\right|}^n$$E2

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, $\dot{x}=\mathit{dx}/\mathit{dt}$ and $\dot{z}=\mathit{dz}/\mathit{dt}$; $k_e=\mathit{\alpha k}$ and $k_h=\left(1-\alpha \right)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].

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{x}\approx \frac{x_{k+1}-x_k}{∆t}$$E3

Where $x_{k+1}$, $x_k$ 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.

Figure 2.

Euler backward approximation.

Now, from Eq. (1), it is possible to define the differential equation $\dot{R}=k_e\dot{x}+k_h\dot{z}$, therefore applying the discrete approximation (3) to $\dot{R},\dot{x}$ and $\dot{z}$ the following equation is applied.

$$R_{k+1}-R_k=k_e\left(x_{k+1}-x_k\right)+k_h\left(z_{k+1}-z_k\right)$$E4

$$z_{k+1}-z_k=A\left(x_{k+1}-x_k\right)-\beta \left|x_{k+1}-x_k\right|{\left|z_k\right|}^{n-1}{z}_k-\gamma \left(x_{k+1}-x_k\right){\left|z_k\right|}^n$$E5

This model can be solved in programming if the deformation $∆x=x_{k+1}-x_k$ is fixed or $∆R=R_{k+1}-R_k$ is fixed, is clear that if $∆x$ is fixed, so $∆R$ is easy to determine. In the sample, k = 0 is clear that $x_0=0$, $R_0=0$ and $z_0=0$, now the next values can be calculated using the next Pseudocode as shown in Table 1.

1. $∆x=x_{k+1}-x_k$ is fixed and given.

2. $A,\beta ,\gamma ,n,\alpha$ and k is given.

3. In k = 0, $z_k=0$, $R_k=0$.

4. Determine $z_{k+1}$ using (5).

5. Determine $R_{k+1}$ using (4).

6. Set $z_k=z_{k+1}$, $R_k=R_{k+1}$

7. Go to (4) to calculate the values of the next sample.

Table 1.

Pseudocode to solve the BW discrete model.

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 \left\{y_1,y_2,y_3,\dots ,y_m\right\}$$E6

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 \left\{x_1,x_2,x_3,\dots ,x_n\right\}$$E7

Decision values ($\stackrel{\sim }{Y}$). These are the values of Decision variables (Y).

$$\stackrel{\sim }{Y}=\left[{\stackrel{\sim }{y}}_1{\stackrel{\sim }{y}}_2\cdots {\stackrel{\sim }{y}}_m\right]$$E8

Design values ($\stackrel{\sim }{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 p-samples required, therefore $\stackrel{\sim }{X}$ is a matrix of $p\times n$ dimension

$$\stackrel{\sim }{X}=\left[\begin{array}{c}\stackrel{\sim }{X}\left(1\right)\\ \stackrel{\sim }{X}\left(2\right)\\ ⋮\\ \stackrel{\sim }{X}\left(p\right)\end{array}\right]=\left[\begin{array}{cccc}{\stackrel{\sim }{x}}_1\left(1\right)& {\stackrel{\sim }{x}}_2\left(1\right)& \cdots & {\stackrel{\sim }{x}}_n\left(1\right)\\ {\stackrel{\sim }{x}}_1\left(2\right)& {\stackrel{\sim }{x}}_2\left(2\right)& \cdots & {\stackrel{\sim }{x}}_n\left(2\right)\\ ⋮& ⋮& \ddots & ⋮\\ {\stackrel{\sim }{x}}_1\left(p\right)& {\stackrel{\sim }{x}}_2\left(p\right)& \cdots & {\stackrel{\sim }{x}}_n\left(p\right)\end{array}\right]$$E9

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 ($\overline{e}$). It is a selected value such that the fitness function evaluation, allows to affirm that the $\stackrel{\sim }{Y}$ values are correct:

$$\sum _{i=1}^p⌊f\left(\stackrel{\sim }{X}\left(i\right),\stackrel{\sim }{Y}\right)⌋<\overline{e}$$E10

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

$$P=\left[\begin{array}{c}\stackrel{\sim }{Y}\left(1\right)\\ \stackrel{\sim }{Y}\left(2\right)\\ ⋮\\ \stackrel{\sim }{Y}\left(q\right)\end{array}\right]=\left[\begin{array}{cccc}{\stackrel{\sim }{y}}_1\left(1\right)& {\stackrel{\sim }{y}}_2\left(1\right)& \cdots & {\stackrel{\sim }{y}}_m\left(1\right)\\ {\stackrel{\sim }{y}}_1\left(2\right)& {\stackrel{\sim }{y}}_2\left(2\right)& \cdots & {\stackrel{\sim }{y}}_m\left(2\right)\\ ⋮& ⋮& \ddots & ⋮\\ {\stackrel{\sim }{y}}_1\left(q\right)& {\stackrel{\sim }{y}}_2\left(q\right)& & {\stackrel{\sim }{y}}_m\left(q\right)\end{array}\right]$$E11

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\left(j\right)=\sum _{i=1}^p⌊f\left(\stackrel{\sim }{X}\left(i\right),\stackrel{\sim }{Y}\left(J\right)\right)⌋$$E12

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}^{\mathit{Ts}}{\stackrel{\sim }{\mathit{y}}}_{\mathit{k}}\left(i\right)}{\mathit{Ts}}{\sigma }_k=\sqrt{\frac{{\sum }_{i=1}^{\mathit{Ts}}{\left({\stackrel{\sim }{\mathit{y}}}_{\mathit{k}}\left(i\right)-{\mu }_k\right)}^2}{\mathit{Ts}}}$$E13

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

$${\stackrel{\sim }{\mathit{y}}}_{\mathit{k}}\left(i\right)=\left\{\begin{array}{c}{\mu }_k+N\left(0,1\right){\sigma }_k\text{if}U\left(.\right)>0.5\\ {\stackrel{\sim }{\mathit{y}}}_{\mathit{k}}\left(i\right)+N\left(0,1\right){\sigma }_k\text{if}U\left(.\right)\le 0.5\end{array}\right.$$E14

Where $N\left(0,1\right)={\sum }_{i=1}^{12}U\left(.\right)-6$, and $U\left(.\right)$ is a random value between 0 and 1.

The pseudocode of the Evonorm algorithm is shown in Table 2.

1. Generation of initial population (11).

2. Evaluation of initial population.

3. Selection of the best $\mathit{Ts}<q$ individuals.

4. Mutation and generation of new population (13) and (14).

5. If a criterion (10) satisfied, then end else go to step 2.

Table 2.

Pseudocode of Evonorm.

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.

Figure 3.

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

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

Figure 4.

Experimental setup for the static monotonic and cyclic tests.

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

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 $\alpha ,\beta ,\gamma ,n,A,k$ are the decision variables:

$$Y\equiv \left\{\alpha ,\beta ,\gamma ,n,A,k\right\}$$E15

$$X\equiv \left\{x_{k+1},x_k,R_{k+1},R_k\right\}$$E16

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::

$$X_v=\left[\begin{array}{cccc}-1.53& -1.51& -10.37& -10.32\\ -1.17& -1.28& -4.61& -5.93\\ -1.04& -1.17& -3.30& -4.61\\ -0.20& -0.39& 2.66& 1.62\\ 0.70& 0.53& 7.32& 6.57\\ 1.36& 1.27& 0.51& 9.21\\ 1.52& 1.64& 10.04& 10.19\end{array}\right]$$E17

Now, it is necessary to determine the fitness function $F\left(\mathit{Xd},\mathit{Yd}\right)$in the algorithm, this function included the design variables $\mathit{Xd}$ and de decision variables $\mathit{Yd}$. Therefore, the fitness function is the approximated model discrete (6).

$$F\left(X,Y\right)=R_{k+1}-R_k-k_e{∆}_X{k}_h\left(A{∆}_X+\beta \left|{∆}_X\right|{\left|z_k\right|}^{n-1}{z}_k\gamma {+∆}_X{\left|z_k\right|}^n\right)$$E18

The major objective will be to minimize the function $F\left(\mathit{Xd},\mathit{Yd}\right)\cong 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:

$$\mathit{\alpha ϵ}\left[0,1\right],\mathit{\beta ϵ}\left[\mathrm{0.1,0.9}\right],\mathit{\gamma ϵ}\left[\mathrm{0.1,0.9}\right],$$E19

$$\mathit{Aϵ}\left[1,10\right],\mathit{nϵ}\left(1,2\right],\mathit{kϵ}\left[10,1000\right]$$E20

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.

α	A	N	β	γ	k	Error
0.40	1.69	1.1	0.81	0.72	7.78	0.21
0.43	1.35	1.1	0.78	0.79	8.53	0.23
0.53	1.98	1.1	0.86	0.64	6.96	0.23
0.46	1.50	1.1	1.12	0.37	7.92	0.25

Table 3.

Evonorm algorithm results.

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 6.

Error graph.

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.

Figure 7.

Comparison between real graph and loops obtained from evolutionary algorithm.

Figure 8.

Graph of actual data vs estimated data

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; 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\left(.\right)>0.4,\text{if}U\left(.\right)\le 0.4$.