Parameter identification of Bouc-Wen model for Magnetorheological (MR) fluid Damper by a Novel Genetic Algorithm

In this research, novel genetic algorithm (nGA) is proposed for Bouc-Wen modle parameters esstimation for magnetorheological (MR) fluid dampers. The optimization efficiency is improved by modifying the crossover and mutation steps of a GA. In the crossover stage, the probability of reproducing offspring from the same parent (same mother and father chromosome) is done to be zero, which may happen in the standard GA, and the probability of a chromosome to be selected for mating is based on error rank weighting of the chromosomes. Additional fitness evaluation of chromosomes will take place in between the crossover and mutation steps to save the best chromosome found so far, which is not implemented in the standard genetic algorithm (GA). The model is validated by comparing its simulation output force (Fsim) with experimentally generated MR damper force (Fexp). The mean absolute error, standard deviation and number of generations for convergence are taken as a criterias for performance evaluation. With these ctriterias, the proposed novel GA outperform better than the other researches. The accuracy is improved by 46.67% compared to standard GA. The proposed novel GA for Bouc-Wen model parameter identification can be used for any MR damper control system with better accuracy.


Introduction
Recently Magnetorheological (MR) damper is one of the best candidates of semi-active nonlinear vibration damping system that can be applied to vibration damping of civil structures, 1,2 automotive vehicles, railway vehicles and Hyperloop capsule train suspension systems 3-6 for its worthwhile features including mechanical simplicity, large force capability, low power requirement, rapid response and robustness. It is also fail-safe, if the control system fails to control the semi-active damper act like passive damper. Hence, MR damper is one of the best research areas for semi-active vibration control. The MR damper fluid is a non-Newtonian, stable rheological properties fluid which consists of magnetic particles of high concentration (26-42% by volume) in a non-magnetic fluid carrier. 7 Most commonly applied magnetic particles are ferrous oxide Fe 2 O 3 and Fe 3 O 4 . The rheological property of MR fluid depends on the concentration of magnetic particles, their shape and size, properties of carrier fluid, applied current to the fluid (magnetic field intensity), temperature and other factors. 7,8 In the absence of an applied current the MR fluid is treated as normal passive suspension system fluid. In this state the magnetic particles are distributed randomly in the fluid as shown in Figure 1(a) and the damper acts as passive damper. When there is an input current to the system, the magnetic particles will form a chain along the flux path in the fluid as shown in Figure 1(b). Such chain arrangement resists and restricts fluid movement as a result of yield stress development in the fluid. 9 The main purpose of vehicle suspension system is to provide a good dynamic performance, such as running stability, ride quality and safety. To fulfill these dynamic performances, the performance and characteristics of suspension system should vary depending on road condition or track irregularity. The conventional passive damper can no longer meet the requirement for high speed railway vehicles. Hence, the semi-active MR damper is introduced to automotive industry as well as railway vehicle suspension system. [10][11][12] In the MR damper suspension system, the inputs are stroke (displacement of MR damper piston) and current from the controller, and output is a continuously varying MR damper force, hence stiffness and damping effect of the MR damper can be continuously adjusted by varying the input current. Based on the ongoing road condition or track irregularity, the controller will manage the amount and characteristics of damping force by controlling the current provided to the MR damper as shown in Figure 2. So that, the dynamic performance of a vehicle will be improved. Interestingly, an automotive vehicles or railway vehicles or capsule trains have an already built-in system to provide the required current.
The general construction of MR damper is shown in Figure 3 above. Structurally, MR damper is similar to the conventional passive hydraulic damper. However, there is an electric circuit in the piston assembly of MR damper to provide current to the MR fluid.
Although the MR damper is the best candidate in vibration damping application, its disadvantage lies in its non-linear force/displacement and hysteretic force/velocity characteristic complexity. In designing controller MR damper, modeling of the actuator is required to be as simple as possible and as much accurate as possible which may be challenging in the case of employing the MR damper. Hence, there should be a trade-off between model accuracy and its complexity. 13,4 The modeling of hysteresis has been studied by different researchers.

Modeling of MR dampers
There are many MR damper modeling approaches proposed in the literatures. In this study, these modeling approaches are divided in to two broad categories, parametric and non-parametric modeling.

Parametric modeling
In parametric modeling approach, the system is idealized to an arrangement of spring and viscos damper which is similar to passive damper. In this approach, the model is represented by mathematical function, and the coefficients of this function can be determined by using an optimization technique, that is, the parameters value are changed until the model's output force is closely match to experimental output force. Some of the models obtained by this approach are Bingham, Viscoelastic-plastic, Bouc-Wen phenomenological, Modified-Bouc-Wen model (Spencer Modeling) and others. 14-18 These models range from simple dryfriction to complicated differential equation representations. Bingham model is simple parametric model which is based on the principle of coulomb friction (dry friction) as shown in Figure 4.
The Coulomb friction element f c is placed parallel to the dashpot c o so for non-zero piston velocity_ x, the damper force F can be expressed as equation (1). 14,15 Viscoelastic-plastic model, which is extension of a Bingham model is proposed in Gomato and Filisko 18 to identify the behavior of electrorheology materials and Spencer adopted this approach this modeling into MR damper. 14 Although the above two models are the simplest MR damper parametric model, they don't provide the force-velocity hysteretic loop and the accuracy is low. Bouc and Wen proposed a Bouc-Wen hysteresis model which possesses an appealing mathematical simplicity and has the ability to represent a large range of hysteresis behavior. 19,20 In the Bouc-Wen model, the mathematical function consists of differential equation to consider the non-linear hysteretic properties of force/ velocity relationships. 14,15 Although Bouc-Wen model is more complex than the Bingham and Bingham-plastic model, it improves the accuracy significantly.
The mechanical idealization of Bouc-Wen model of MR damper is shown in Figure 5, the restoring force F can be predicted by equation (2), Where the evaluation variable z is described by another differential equation (3), Where F: Damping force; c o ; k o ; f o ; a; d; b; g and n are the parameters of shape of Bouc -Wen model to be identified. Though the Bouc-Wen model is extremely versatile and can exhibit a wide variety of hysteretic behavior, the estimated result of the model may not match close to the experimental data in the regions where acceleration and velocity has opposite sign and the magnitude of the velocity is small. 14,15 In order to improve the accuracy at this region, Spencer proposed a modified Bouc-Wen model by introducing another internal displacement 14,15 and modifying the Bouc-Wen model structure by adding dashpot and spring as shown in Figure 6.
The damping force F is calculated as equation (4).
In order to account for the additional internal degree of freedom introduced in the Bouc-Wen model, an additional first order differential equation is required. However, modified Bouc-Wen modeling improves accuracy, the MR damper will be more complex for MR damper control design, and expensive owning to an addition of displacement sensor. The improvement of error found by using modified Bouc-Wen model is insignificant when compared to the MR damper complexity and expensiveness.

Non-parametric modeling
Non-parametric models are entirely based on the performance of a specific MR fluid device. Commonly an elevated amount of experimental data, obtained by observing the magnetorheological response to different excitations under varying operating conditions, is used to predict the device's response to random excitations. 15 Neural-network model and fuzzy model, which are a type of nonparametric modeling, were proposed to model MR damper. [21][22][23][24] Even though both are applicable in modeling of MR damper they require tuning or training the networks in order to derive a feasible model. 16 There are also many other nonparametric modellings proposed by different researchers, i.e., curve fitting modeling, 25,26 black box modeling 27 and nonparametric approach. 28 These groups are more flexible but the physical relationship between modal parameters and hysteresis phenomena may not be explicitly maintained. 16 As a result of all the above-mentioned advantages and drawbacks of both parametric and non-parametric modeling, Bouc-Wen parametric modeling of MR damper has been used extensively to simulate hysteresis loop. Hence, the Bouc-Wen model for MR damper is used in this article.

Parameter identification of Bouc-Wen model for MR Damper
Evolutionary computing techniques such as genetic algorithm, particle swarm optimization, and simulated annealing are widely used for parameter identification of Bouc-Wen model of MR fluid dampers owing to they don't require derivative information, may encode the variables so that the optimization is done with the encoded variables and work with numerically generated data, experimental data or analytical function.
A modified parametric algebraic model is proposed to predict the behavior of hysteresis of prototyped MR damper, and the performance of the model is compared to a Bingham plastic flow model, and a modified Bouc-Wen model through error analysis. 29 In Kwok et al., 28 a hyperbolic tangent function is proposed to describe a Bouc-Wen model hysteresis component to simplify the model, and the parameters are identified using particle swarm optimization. In Kwok et al., 16 a computational efficient genetic algorithm is proposed to identify of Bouc-Wen model parameters of MR damper. Although computational time is minimum and the model is not complex in both the above two approaches, the accuracy is not somehow good. A Bouc-Wen model with adaptive charging system optimization is proposed in Talatahari et al., 30 and a satisfactory result was achieved. A Bouc-Wen-Baber-Noori (BWBM) modeling with the force-lag phenomenon is proposed to describe the distorted behavior of self-constructed MR damper in Peng et al. 13 In Zhu et al., 31 the normalized Bouc-Wen model, which is a variant of the Bouc-Wen model with fewer parameters, is proposed. Although in the normalized Bouc-Wen model parameters to be identified is decreased by one, the outcome error is larger than the outcome error of this paper. The error comparison is done in section 6. In Giuclea et al., 32 a method of finding the modified Bouc-Wen model parameters is proposed by using a genetic algorithm (GA). In Charalampakis and Koumousis, 33 a new method of parameter identification of the Bouc-Wen model is proposed. It is based on a hybrid evolutionary algorithm that utilizes selected stochastic operators, and problemspecific information. In Xiaomin et al., 34 the authors proposed modified genetic algorithm and a simplified Bouc-Wen model for MR damper. They did a sensitivity analysis and concluded that c o ; k o ; f o ; andn parameters significantly affect the shape of the hysteresis loop at different input current, so these parameters are considered to be functions of the excitation current, and the other four parameters are taken as a constant. In Yang et al., 12 the author proposed particle swarm optimization for parameter identification of Bouc-Wen model of an MR damper. In this research, the author regarded the initial damper displacement which contribute the force offset as zero and the exponent term of the Bouc-Wen model, n = 2. So, the parameters to be optimized is decreased to six variables, which will have bad effect on the optimization efficiency. Regarding to Particle swarm optimization, it is similar to genetic algorithm in some ways but one may perform better over the other based on the given data for optimization.
Although many research works have been done for the improvement of magnetorheological damper dynamic control performance, they have some limitations either in the accuracy or in their complexity to design the control system. In this study, a Bouc-Wen model with novel genetic algorithm parameter identification is proposed for MR damper dynamics analysis. The novel genetic algorithm parameter identification method proposed in this paper has the following advantages: (1) All parameters are assumed to affect the shape of the hysteresis loop, though the degree is different. So, the error of the model will be minimized. (2) To improve computation efficiency of the genetic algorithm, the selection, crossover, and mutation stage are modified according to the type of experiment data used for modeling MR damper. (3) The error, which is found with this proposed novel GA, is lower than the one found by the standard genetic algorithm, efficient genetic algorithm proposed in Kwok et al. 16 and normalized Bouc-Wen model. 31 The computational time is also minimum compared to those three models. (4) All the parameters are found to be linearly related to the excitation current. So, MR damper controller design will be easy and more accurate.
The paper is organized in the following section. Section (4) Experiment setup and data extraction; Section (5) Describing the noble genetic algorithm working procedure and find out the parameters of the model to the experimental data; Section (6) Find out the error of simulation result of this model and compare it with the errors found by standard genetic algorithm and efficient genetic algorithm proposed by Kwok et al. 16 and Zhu H et al. 31 ; Section (7) Conclusion.

Experiment setup and data extraction
A MR damper testing setup is arranged as shown in Figure 7, and the experiment is performed by applying a sinusoidal excitation of magnitude 50-250 mm at a low frequency of 1-2 Hz to the MR damper. 35 The damper forces, which are generated under the application of different currents (i.e., 0-1.2A), are measured by the load cell. Force, displacement, and velocity readings are recorded for parameter identification. Figure 8 represents the nonlinearity relationship of force/displacement and the hysteretic relationship of force/velocity of the experiment data.

Parameter identification by noble genetic algorithms
Genetic Algorithms (GAs) are search based algorithms based on the concepts of natural selection and genetics.
GAs are subsets of a much larger branch of computations known as Evolutionary Computations. A GA allows a population composed of many individuals to evolve under specified selection rules to a state that   maximizes the ''fitness'' (i.e., minimize the cost function). 36 Figure 9 shows the general computational layout of novel GA in this paper. Activities performed in each stage are discussed below.

Population initialization
The GA begins by defining a group of chromosomes known as a population. A chromosome C var is an array of variable values to be optimized. In this paper, there are eight Bouc-Wen model parameters to be identified as shown in equations (2) and (3). Therefore, every chromosome has eight variables (genes) of real numbers. The number of chromosomes in this paper is chosen as N pop = 100, hence population size of the GA is 100 chromosomes.

Fitness function calculation and survival selection
A Fitness function is used to evaluate optimality of each chromosome (a set of variables or parameters) in each iteration. Survival selection is the process of selecting the fitter chromosomes and discarding the chromosomes with the highest error. The design of fitness function is also problem dependent, and the fitness function for this paper is designed from the average sum of normalized square error between the simulated and experimental data. 16 Where n is the number of data points and _ x i 2 À1; 1 ½ is the normalized velocity, given by Superscripts sim and exp denote simulation and experimental forces,Df = f max À f min is the difference between the maximum and minimum forces of the experimental force data, and l is the number of experimental data. Each simulation/experiment data is indexed by subscript i. The experimental force f exp i ð Þwill be taken from the experimental data and simulation (predicted) force f sim i À Á is generated from the Bouc-Wen MR damper model using 4th order Rung-Kutta routine with chromosomes (a set of parameters) C var, j obtained from the proposed novel GA (using equations (2) and (3)). So, the error e j of every chromosome will be calculated using the fitness function (equation (7)). Note that the smaller error indicates a fitter chromosome and a better matching between the estimated and the experimental force. After fitness evaluation,N pop errors and associated chromosomes are ranked from lowest error to highest error so that each chromosome is better than the next one. Then, only the first N keep chromosomes are selected to continue, while the rest are discarded. The selection rate, X rate , is the fraction of N pop that survives for the next step of mating. The number of chromosomes that are kept at each generation is Survival selection takes place at each generation or iteration of the algorithm.

Cross over, mutation and elitism steps
Crossover and mutation steps are the most important steps in GA that directs the search for solutions by exploration and exploitation and are used to refine the best solution found so far. In cross over step, parent selection and mating operations will take place at this stage to replace the discarded chromosomes in the survival selection stage. In the parent selection operation, mother and father chromosomes are selected from the survived chromosomes. There are a variety of parent selection methods such as roulette wheel, stochastic reminder, rank weighting, tournament selection, and others. In the rank weighting selection method, the probabilities assigned to the chromosome in the mating pool are inversely proportional to their error. A chromosome with the lowest error has the greatest probability of mating, while the chromosome with the highest error has the lowest probability of mating. Rank weighting is also problem independent and finds the probability from the rank, k, of the chromosomes. Due to the above advantageous features, the rank weighting method is used for parent selection in this paper.
Where: P k is the probability of k th chromosome to be selected. Cumulative probability will then be Therefore, the probability of selection will take place in the following fashion. Take a random number r such that r 2 0; 1 ½ Then the m th chromosome will be selected as a parent (mother and father) chromosome for the reproduction of offspring. Parent selection will take place until the amount of mother and father chromosomes is equal to Mating is the creation of one or more offspring from the parents selected in the parent selection process. In the mating stage, each mother and father chromosomes are paired in random fashion and mated to get an offspring as shown below.
Where; ma and pa are randomly selected mother and father chromosomes respectively and r is any random number such that r 2 0; 1 ½ . In standard GA, there is a probability that ma and pa chromosomes may be the same. Hence, the reproduced offspring is also the same. So at least four chromosomes in the next generation will be the same. This situation will reduce the efficiency of standard GA.
Hence in the proposed novel GA, it is prohibited to select the same chromosome as a mother and father chromosome.
In the former research work, 16 the number of crossovers at each iteration is defined probabilistically, which will influence the efficiency of the GA. To alleviate this problem, a predefined number of crossovers will take place at each iteration in this paper.
A mutation is the second way a GA explores solution space. It introduces trials, not in the original population, and keeps the GA from converging too fast before sampling the entire solution space, hence the algorithm is not trapped in a local minimum. Increasing the number of mutations increases the algorithm's freedom to search outside the current region of variable space. A problem-dependent mutation rate, Y rate , a fraction of ðN pop Ã C vr Þ;is decided in this paper. Then, N pop Ã C vr Ã Y rate variables (genes), Var, of a population are selected randomly, and a Gaussian random number is added to them to mutate the chromosomes.

Var i ð Þ Var i ð Þ + rl
Where r 2 0; 1 ½ is uniform random number and rl;N 0; s 2 ð Þ is a normally distributed random number with s as the standard deviation of the initial population. 16 As shown in Figure 9, in this paper, an additional fitness evaluation of chromosomes will take place in between the crossover and mutation stages to save the best chromosome found after crossover from being mutated, which is not implemented in standard GA and efficient GA. 16 It increases the efficiency the GA.
Additionally, in this paper, selection for mutation takes place at a variable level, and not at a chromosome level. But in Kwok et al., 16 the whole variable values of a chromosome will be mutated if a chromosome is selected for mutation, which may have some impact on the efficiency of the GA.
As a result of the application of elitism, the best candidate chromosome obtained during GA iteration is copied to the next offspring without modification by crossover or mutation. Therefore, the population keeps at least one copy of the fittest chromosome found so far, and the smallest error cannot be increased over iterations, that is Where, e q min is the minimum error at the q th iteration. The current running average minimum error is calculated as The minimum error and running average minimum error vs. generations are plotted as shown in Figure 10.

Termination
There are many possible ways to terminate the GA based on the problem domain, such as termination after reaching the maximum limit of generation, termination after reaching the predefined minimum error of the fitness function, termination if there is no improvement in the best fitness values (smallest error values) through a predetermined number of consecutive iterations, and termination when the chance of achieving significant changes of fitness function error in the next generation is excessively low. The method used by Kwok et al. 16 for the decision of termination of GA iteration is also used here in this paper.
By using the principle of maximum entropy, the minimum error is represented by exponential distribution as shown in Figure 11, where the distribution is given as a histogram of 10 bins. The exponential distribution can be given as Pe = l À1 exp l À1 e i min À Á ; i = 1:::q ð17Þ Where l = e q min is the running average minimum error and e i min the minimum error of the i th generation To make the error exponentially distributed a sufficient minimum sample shall be available, and in this paper q min .15 is taken as the minimum iteration for termination.
Considering the behavior of the error's cumulative probability distribution the GA iteration is proposed to terminate when the minimum error reduction greater than 20% of the maximum error is not expected in further iterations. Hence, the termination criteria in this paper is defined as Finally, the GA will terminate at generation q such that Where L is the logical and operator

Result and discussion
Experimental data, which are extracted from the experiment under different combinations of inputs, are fed to the proposed novel GA, efficient GA proposed on article, 16 normalized Bouc-Wen model proposed on article, 31 and standard GA to find out model parameters. Then, the simulation forces generated by each method are compared to the experimentally obtained force to evaluate the performance and accuracy of each of the above proposed methods.
In the normalized Bouc-Wen model which is stated in Zhu et al., 31 the general equation is stated in equation (20).
The main objective of the proposed method in Zhu et al. 31 is decreasing the number of parameters to be identified to seven, and using three steps to find out the optimum value of each parameter. In the first step the value of k 1 is determined considering two wave T-periodic displacement signals. In the second step, the optimum values of c 1 ; a and f o are determined by least square method. In the third step the values of r; s and n are determined by standard genetic algorithm. The detail optimization work can be referred in Zhu et al. 31 Error and computational time comparison The reconstructed force/velocity hysteresis from the above four methods of Bouc-Wen model simulation results and from experimental data for randomly selected input currents of 0.2A and 0.8A are compered graphically as shown in Figure 12.  Figure 11. Distribution of minimum error and accumulative distribution.
The comparison between simulation force (F sim ), found by Bouc-Wen model with the above four parameter identification methods, and experimentally recorded force (F exp ) allows us to compare the accuracy of each method. The error comparison is done by equations (21) and (22). The errors of standard GA, the proposed efficient GA of paper, 16 efficient parameter identifications of normalized Bouc-Wen model proposed on paper, 31 and the proposed novel GA of this paper are listed in Figure 13 below. From the given bar graph the performance of the novel GA proposed in this paper is better than the other methods proposed on papers 16,31 and standard GA owing to its error value and number of generation (computation time).
Where MAE: mean absolute error N: number of samples of experimental data where : In addition to accuracy performance, the proposed novel GA converges faster. It converges faster than the other methods as shown in Figure 10 and 13c.

Relationship between identified parameters and input current
Parameters are identified by using different sets of experimental data and current input. It is found that every parameters value vary as the current is changed as shown in Table 1. So, the relationship between current and parameter values is found using Lagrange's polynomial interpolation based on supplied current.

Conclusion
This paper has presented Bouc-Wen model for magnetorheological fluid damper modeling. The parameters of the model are identified using different sets of experimental data with the implementation of proposed novel GA. The efficiency of the proposed novel GA is increased by using Rank weighting approach in the selection of chromosomes for mating and prohibiting same chromosomes (mother and father) to be mate to produce an offspring. An additional best chromosomes selection takes place after crossover at each iteration to save the best one from mutation. Mutation takes place in a fashion different from the common one. Setting the initial boundaries of a population has an indispensable role to converge the algorithms faster. So, in this work the initial boundary is set after two or three trails.    The proposed novel GA optimization performance is evaluated with respect to the experimental data and found that the error produced by this novel GA is smaller than the error produced by standard GA, efficient GA proposed in Kwok et al. 16 and efficient normalized Bouc-Wen model proposed by Zhu et al. 31 as shown in Table  2 and the computational time is better than the other three as shown graphically in Figure 13.

Declaration of conflicting interests
The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding
The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This paper has been supported by Core Technology Development of Subsonic Capsule Train research project, Kora Railroad Research Institute, Korea as part of grant No. PK2001A1A.