# Modeling and Parameter Identification of MR Damper considering Excitation Characteristics and Current

Authors: Shuguang Zhang; Wenku Shi; Zhiyong Chen
Journal: Shock and Vibration
 (2021)

---

## Abstract

Smart structures such as damping adjustable dampers made of magnetorheological (MR) fluid can be used to attenuate vibration transmission in vehicle seat suspension. The main research content of this paper is the nonlinearity and hysteresis characteristics of the MR damper. A hysteretic model considering both excitation characteristics and input current is proposed to fit the damper force-velocity curve for the MR damper under different conditions. Multifactor sensitivity analysis based on the neural network method is used to obtain importance parameters of the hyperbolic tangent model. In order to demonstrate the fitting precision of the different models, the shuffled frog-leaping algorithm (SFLA) is employed to identify the parameters of MR damper models. The research results indicate that the modified model can not only describe the nonlinear hysteretic behavior of the MR damper more accurately in fixed conditions, compared with the original model, but also meet the fitting precision under a wide range of magnitudes of control current and excitation conditions (frequency and amplitude). The method of parameter sensitivity analysis and identification can also be used to modify other nonlinear dynamic models.


---

## Body

## 1. Introduction

Vehicle ride comfort plays an important role in vehicle dynamics research, in which the transmission of vibration through seat exerts a crucial influence. Traditional passive seat suspensions are very limited in terms of ride comfort and cannot meet the increasingly stringent comfort requirements. In recent years, the semiactive seat suspension system is becoming more and more popular due to simple structure, low power consumption, good vibration damping behavior, and high reliability [1–6].There are many types of semiactive suspensions such as throttle variable damper suspension, air suspension, electrorheological damper suspension, and magnetorheological damper suspension [7, 8]. Therefore, an MR damper is a fluid damping device, widely used in vibration control of civil engineering structures and mechanical systems because of its wide damping adjustable range, short adjustment time, high reliability, and low cost [9–12].The structure of the magnetorheological damper is similar to the structure of an ordinary vehicle cylinder damper. The medium inside the MR damper is a magnetorheological fluid. The MR fluid can reversibly change from free flowing, linear viscous liquids to non-Newtonian fluid with a controllable yield strength when exposed to a magnetic field [13].However, the magnitude of the generated damper force is related to many factors that involve many nonlinear elements and increase model complexity of the MR damper. The causes of nonlinearity are listed as follows: (1) different magnetoresistivity of electromagnetic circuits; (2) hydrodynamics of MR fluid flowing in piston gap; (3) the chemical properties of MR fluid; and (4) the relationship between the yield stress generated by MR fluid and applied magnetic field strength.Various papers have established MR damper mechanical models by various proposed methods. According to relevant literature, MR damper force models are divided into three main groups [10, 14–16]: pseudostatic models, parametric models, and nonparametric models. In pseudostatic models, the damper force is calculated by the pressure difference between the two ends of the damper piston, structural parameters, MR fluid parameters, and magnetic field strength [17]. This type of model can generally reflect the force-displacement characteristics properly, but the fitting of the force-velocity characteristics is not accurate enough. Pseudostatic models are mainly used for the design of MR dampers, but not suitable for the research of semiactive suspension control. Parametric models and nonparametric models are mainly distinguished according to whether the parameters make sense in physical. Nonparametric models generally consider that the damper force consists of a polynomial with input current and damper velocity, or neural network model [18]. And, in parametric models, the damper force is calculated by physically meaningful parameters. There are many types of parametric models, among which the Bingham model is the most basic one. Although this model can describe the relationship between the damper force and displacement well, it cannot reflect the hysteresis relationship of force-velocity. At present, the Bouc–Wen model is widely used, which can exactly describe the hysteretic characteristics of the MR damper force versus piston velocity by designing the hysteresis variable z. However, the nonlinear differential equations contained in the Bouc–Wen model inevitably increases difficulty to parameter identification. In 1997, Spencer and Dyke [19] proposed a modified Bouc–Wen model, which further improved the fitting accuracy of the force-velocity hysteresis curve. There is a wide range of parametric models: nonlinear hysteretic biviscous model, generalized sigmoid hysteresis model, hyperbolic tangent function-based model, and so on. However, the previous models have shortcomings in terms of model accuracy, complexity, or limitation [14, 16, 20–22].The purpose of this paper is to build a detailed model of the MR damper used for semiactive suspension control. The paper is structured as follows: in Section2, the mechanical performance test of the MR damper was performed to study the characteristics of damper force in different excitation inputs. In Section 3, several classic MR damper models were introduced in detail, and the more suitable model was selected. In Section 4, the sensitivity analysis method of parameters of the hyperbolic tangent model was applied to illustrate the importance of every model parameter, as references for model extensions in subsequent chapters. The shuffled frog-leaping algorithm used for parameter identification is explained in Section 5. Subsequently, the identification results of damper hysteretic model parameters are discussed and an extended hysteresis model is achieved and verified in Section 6. In the end, conclusions are presented in Section 7.


## 2. Experimental Study of the MR Damper

### 2.1. Experimental Setup

Two types of MR dampers, shown in Figure1, were employed for dynamic tests in this study. The strokes of the pistons of these two dampers are 74 mm (long stroke, RD-8041-1) and 55 mm (short stroke, RD-8040-1), respectively. Except the length of stroke, the two dampers have the same structure parameters and MR fluid. The damper consists of two liquid chambers that are separated by a piston with small orifices. The rated value of the input current is 1 A in the electromagnet coil around the piston head.Figure 1
Picture of the MR damper.Figure2 shows the experimental setup including actuator, computer closed-loop control electro-hydraulic servo system, force sensor, displacement sensor, fixture, and constant current power supply. The damper was driven by a Schenck actuator, and the position sensor and force transducer were used to collect the measured displacement and damper force of the two dampers.Figure 2
Experimental setup.Sine vibrations, with the frequency range from 0.5 Hz to 9 Hz and vibration amplitude of 5 mm, 10 mm, 15 mm, and 20 mm, were provided by the actuator. And, the current inputs ranged from 0 to 1 A with 0.2 A intervals. The specification of excitations is shown in Table1.Table 1
Specification of excitation.
Amplitude in (mm)Frequency51015200.5 Hz√√√√1 Hz√√√√1.5 Hz√√√√2 Hz√√√√2.5 Hz√√√√3 Hz√√√√3.5 Hz√√√√4 Hz√√√√4.5 Hz√√√—5 Hz√√√—6 Hz√√——7 Hz√√——9 Hz√——


### 2.2. Experimental Results

In the following analysis, only the results of the RD-8040-1 damper were discussed, in which the mechanical behaviors of the two dampers were similar, as was shown in the experiment. The responses of the MR damper applied a 20 mm 2 Hz sinusoid excitation which are shown in Figure3, at six constant current values, 0 A, 0.2 A, 0.4 A, 0.6 A, 0.8 A, and 1 A. Figures 3(a) and 3(b) illustrate the force-displacement (f-d) and the force-velocity (f-v) characteristics, respectively. It is obvious that the f-d curve is approximately rounded rectangle. There are evident bumps near the zero-displacement area, which is caused by the elasticity of the damper. The effect of changing current value is readily observed, and the damper force curve rises monotonically as the current value increases. The increase of current value is found to be nonlinear in nature. Under 20 mm 2 Hz sinusoid excitation, the maximum damper force increases by 450.71N when current value increases from 0.2 A to 0.4 A. Similarly, for increase from 0.4 A to 0.6 A and from 0.6 A to 0.8 A, the maximum damper force increases by 297.02 N and 324.06 N severally. When the current value increases from 0.8 A to 1 A, the increment of maximum damping force is 64.52N. It can be also observed that the magnetic field intensity generated by the electromagnetic coil is close to the magnetic saturation intensity once the current value increases above 0.8 A, reflected in which the increment of the damper force of the MR damper becomes smaller.Figure 3
Mechanical properties of different currents, frequency = 2 Hz, and amplitude = 20 mm. (a) Damper force-displacement. (b) Damper force-velocity.
(a)(b)The damper force-velocity curve is symmetrical and nonlinear, with apparent hysteresis characteristics. When the velocity of the piston increases, the damper force increases along the lower branch curve. The hysteresis curve of damper force moves counterclockwise with damper velocity. When the input current value of the coil increases, the hysteresis loop of the damper force-velocity of the MR damper gradually increases.Figure4(a) shows the force-displacement characteristics with 0.8 A current applied under 5 mm sine vibration and various frequency (1 Hz, 2 Hz, 3 Hz, 4 Hz, 5 Hz, 6 Hz). It can be observed that the maximum damper force increases when increasing frequency, keeping the current and amplitude invariable. As the frequency grows higher, the area enclosed by the damper force-displacement loops increases, and the energy dissipation in each vibration cycling increases. As is shown in Figure 4(b), the damper force at the same velocity is approximately equal and shows the Bingham plastic behavior of the MR damper when the velocity and acceleration are in the same direction. However, the damper force in the same velocity decreases when frequency grows higher, in case the velocity and acceleration are reversed. This phenomenon of the MR damper is because the MR fluid presents viscoelasticity at low velocity and viscoplasticity at a high velocity. Furthermore, the size of the hysteresis loop mainly depends on the inertia effect of damper. The area of hysteresis loops increases due to the large inertia resulting from the high acceleration as the velocity direction changes.Figure 4
Mechanical properties of different frequencies, current = 0.8 A, and amplitude = 5 mm. (a) Damper force-displacement. (b) Damper force-velocity.
(a)(b)In Figures5(a) and 5(b), the area enclosed by the damper force-displacement curve increases as the excitation amplitude increases under the same input current and frequency, and the area of force-velocity hysteresis loop becomes larger as well.Figure 5
Mechanical properties of different amplitudes, current = 0.8 A, and frequency = 2 Hz. (a) Damper force-displacement. (b) Damper force-velocity.
(a)(b)Based on the characterizations of the MR damper above, it can be observed that the MR damper behaves strong nonlinearity. All of the excitation amplitudes, frequencies, and input currents have impacts on the mechanical characteristics of the MR damper.It is difficult to accurately describe the MR damper model by conventional methods. The establishment of a comprehensive dynamic model of the MR damper has an important influence on the control of the semiactive seat suspension system. In this study, a detailed hysteretic model was established to describe the mechanical characteristics of the MR damper for the purpose of vehicle suspension control.


## 2.1. Experimental Setup

Two types of MR dampers, shown in Figure1, were employed for dynamic tests in this study. The strokes of the pistons of these two dampers are 74 mm (long stroke, RD-8041-1) and 55 mm (short stroke, RD-8040-1), respectively. Except the length of stroke, the two dampers have the same structure parameters and MR fluid. The damper consists of two liquid chambers that are separated by a piston with small orifices. The rated value of the input current is 1 A in the electromagnet coil around the piston head.Figure 1
Picture of the MR damper.Figure2 shows the experimental setup including actuator, computer closed-loop control electro-hydraulic servo system, force sensor, displacement sensor, fixture, and constant current power supply. The damper was driven by a Schenck actuator, and the position sensor and force transducer were used to collect the measured displacement and damper force of the two dampers.Figure 2
Experimental setup.Sine vibrations, with the frequency range from 0.5 Hz to 9 Hz and vibration amplitude of 5 mm, 10 mm, 15 mm, and 20 mm, were provided by the actuator. And, the current inputs ranged from 0 to 1 A with 0.2 A intervals. The specification of excitations is shown in Table1.Table 1
Specification of excitation.
Amplitude in (mm)Frequency51015200.5 Hz√√√√1 Hz√√√√1.5 Hz√√√√2 Hz√√√√2.5 Hz√√√√3 Hz√√√√3.5 Hz√√√√4 Hz√√√√4.5 Hz√√√—5 Hz√√√—6 Hz√√——7 Hz√√——9 Hz√——


## 2.2. Experimental Results

In the following analysis, only the results of the RD-8040-1 damper were discussed, in which the mechanical behaviors of the two dampers were similar, as was shown in the experiment. The responses of the MR damper applied a 20 mm 2 Hz sinusoid excitation which are shown in Figure3, at six constant current values, 0 A, 0.2 A, 0.4 A, 0.6 A, 0.8 A, and 1 A. Figures 3(a) and 3(b) illustrate the force-displacement (f-d) and the force-velocity (f-v) characteristics, respectively. It is obvious that the f-d curve is approximately rounded rectangle. There are evident bumps near the zero-displacement area, which is caused by the elasticity of the damper. The effect of changing current value is readily observed, and the damper force curve rises monotonically as the current value increases. The increase of current value is found to be nonlinear in nature. Under 20 mm 2 Hz sinusoid excitation, the maximum damper force increases by 450.71N when current value increases from 0.2 A to 0.4 A. Similarly, for increase from 0.4 A to 0.6 A and from 0.6 A to 0.8 A, the maximum damper force increases by 297.02 N and 324.06 N severally. When the current value increases from 0.8 A to 1 A, the increment of maximum damping force is 64.52N. It can be also observed that the magnetic field intensity generated by the electromagnetic coil is close to the magnetic saturation intensity once the current value increases above 0.8 A, reflected in which the increment of the damper force of the MR damper becomes smaller.Figure 3
Mechanical properties of different currents, frequency = 2 Hz, and amplitude = 20 mm. (a) Damper force-displacement. (b) Damper force-velocity.
(a)(b)The damper force-velocity curve is symmetrical and nonlinear, with apparent hysteresis characteristics. When the velocity of the piston increases, the damper force increases along the lower branch curve. The hysteresis curve of damper force moves counterclockwise with damper velocity. When the input current value of the coil increases, the hysteresis loop of the damper force-velocity of the MR damper gradually increases.Figure4(a) shows the force-displacement characteristics with 0.8 A current applied under 5 mm sine vibration and various frequency (1 Hz, 2 Hz, 3 Hz, 4 Hz, 5 Hz, 6 Hz). It can be observed that the maximum damper force increases when increasing frequency, keeping the current and amplitude invariable. As the frequency grows higher, the area enclosed by the damper force-displacement loops increases, and the energy dissipation in each vibration cycling increases. As is shown in Figure 4(b), the damper force at the same velocity is approximately equal and shows the Bingham plastic behavior of the MR damper when the velocity and acceleration are in the same direction. However, the damper force in the same velocity decreases when frequency grows higher, in case the velocity and acceleration are reversed. This phenomenon of the MR damper is because the MR fluid presents viscoelasticity at low velocity and viscoplasticity at a high velocity. Furthermore, the size of the hysteresis loop mainly depends on the inertia effect of damper. The area of hysteresis loops increases due to the large inertia resulting from the high acceleration as the velocity direction changes.Figure 4
Mechanical properties of different frequencies, current = 0.8 A, and amplitude = 5 mm. (a) Damper force-displacement. (b) Damper force-velocity.
(a)(b)In Figures5(a) and 5(b), the area enclosed by the damper force-displacement curve increases as the excitation amplitude increases under the same input current and frequency, and the area of force-velocity hysteresis loop becomes larger as well.Figure 5
Mechanical properties of different amplitudes, current = 0.8 A, and frequency = 2 Hz. (a) Damper force-displacement. (b) Damper force-velocity.
(a)(b)Based on the characterizations of the MR damper above, it can be observed that the MR damper behaves strong nonlinearity. All of the excitation amplitudes, frequencies, and input currents have impacts on the mechanical characteristics of the MR damper.It is difficult to accurately describe the MR damper model by conventional methods. The establishment of a comprehensive dynamic model of the MR damper has an important influence on the control of the semiactive seat suspension system. In this study, a detailed hysteretic model was established to describe the mechanical characteristics of the MR damper for the purpose of vehicle suspension control.


## 3. The Mechanical Models of the MR Damper

As a controllable fluid, MR fluid is a Newtonian fluid that can freely flow in the initial state and turns into non-Newtonian fluid when exposed to a magnetic field. The change of MR fluid viscous is reversible and reflected by a significant variation in the damper output force.A lot of models have been proposed to describe the nonlinear behavior of the damper force versus velocity of MR dampers. There are a couple of models that were commonly employed in the previous work such as the Bingham model, Bouc–Wen model, nonlinear hysteretic biviscous model, generalized sigmoid hysteresis model, hyperbolic tangent function-based model, and so on [23–29].Bouc–Wen model is widely used in the modeling of hysteretic systems (Figure6). It can exhibit various hysteretic behaviors by changing its model parameters and has very good adaptability. The expression is shown as follows:(1)F=c0x˙+k0x−x0+αz,z˙=−γx˙zzn−1−βx˙zn+Ax˙,where, α, β, γ, and n are the model parameters and z is a hysteresis variable. The bias force generated by the accumulator inside the MR damper is converted to the initial offset x0 of the spring k0. The Bouc–Wen model has good fitting accuracy for the damper force-velocity characteristics of MR dampers. However, it makes great difficulty to the parameter identification and may affect the robustness of the control system, due to the existence of the differential term z˙.Figure 6
Bouc–Wen model.Research shows that the Bouc–Wen model can better predict the force-displacement and force-velocity characteristics of MR dampers, compared with the Bingham model. But, in the low-speed area, the model is not accurate enough to fit the velocity hysteresis characteristics.Spencer proposed a modified Bouc–Wen model in order to more exactly predict the damper response in low velocity. The damper force is given as(2)c0y˙=αz+k1x−xin+c1x˙−x˙in,z˙=−γx˙−x˙inzzn−1−βx˙−x˙inzn+Ax˙−x˙in,F=c0x˙in+k0x−x0,where c0 and c1 represent the viscous damping at low velocity and high velocity, respectively; k1 represents stiffness at high speed; and xin is the internal variate of displacement, as shown in Figure 7. Other parameters are consistent with the Bouc–Wen model. Although the accuracy of the model is improved, the complexity of the model also increases significantly. The newly brought differential intermediate variable in this model brings greater challenges in parameter identification and the robustness of control system.Figure 7
Modified Bouc–Wen model.Kwok and Ha [30] proposed a damper model in 2006. They used a hyperbolic tangent function to represent the hysteresis element z (see Figure 8). The viscosity and stiffness elements are expressed by linear functions c0x˙ and k0x. Due to lack of differential equations, the calculation efficiency of this model is promoted, and parameter identification becomes simple to operate. The model is expressed as follows:(3)z=tanhβx˙+δsignx,F=c0x˙+k0x+αz+F0.Figure 8
Hyperbolic tangent model.Through calculation and analysis, it can be found that the hyperbolic tangent model can accurately fit the velocity-damping force hysteresis curve of the MR dampers. Figure9 shows the comparison of the calculation results of the two models. The fitting accuracy of the hyperbolic tangent model is higher than the Bouc–Wen model especially in the velocity transition zone. The output damping force error RMS of the two models compared to the test damping force are 30.88N and 49.49N, respectively.Figure 9
The calculation results of different models.By comparing different hysteretic models, it can be found that the hyperbolic tangent model is simple in structure and does not contain differential terms. Meanwhile, it can fit the velocity-damping force hysteretic characteristics of the MR damper precisely. Therefore, the hyperbolic tangent model is selected for extended research in this paper.


## 4. Parameter Sensitivity Analysis

Sensitivity analysis of the mathematical model parameters is not only critical to model validation but also can provide guidance for research of parameters analysis. The analysis of hyperbolic tangent hysteresis model shows that it can predict response for fixed current input and vibration excitation. However, the model parameters under different currents and excitation conditions are different. In order to establish a general model suitable for different working conditions, a hysteresis model including current, amplitude, and frequency characteristics is needed. The method of sensitivity analysis can distinguish the important parameters from unimportant parameters in the mathematical model [31]. Due to the large number of parameters in the hyperbolic tangent model, the sensitivity analysis about the contribution of the parameters to the MR damper model output is carried out to filter out the parameters with poor sensitivity. The relationships between the sensitive parameters and the excitation characteristics or input current are established, so as to obtain the extensional MR damper model that is suitable for different excitation conditions.In sensitivity analysis, the sensitivity coefficient can be expressed as follows:(4)S=1N∑i=1NFi′−Fi20.5,where Fi and Fi′represent the model output force values at the ith sampling spot of parameters before and after changing; N = total number of sampling spot; and S = sensitivity of parameters. The sensitivity of the parameter increases as the value of S becomes larger.The main parameters that have relatively significant influence on the sensitivity coefficient of damper force were obtained via adjusting each parameter within its allowable range, so as to study the degree of influence of these parameters on the accuracy of the model. Furthermore, in subsequent model extensions, model design was based on the ranking of parameter sensitivity, focusing on parameters with larger sensitivity coefficients, and establishing the relationships between model parameters and excitation characteristics or input current. In this way, the number of model parameters is minimized as much as possible on the premise of ensuring model accuracy. Therefore, the complexity of the model is reduced, and the computational efficiency of the model is improved.There are plenty of sensitivity techniques in previous research, and the results of sensitivity ranking may exist slightly disparity for different methods [32]. In this paper, multiparameter sensitivity analysis based on BP neural network was put in use.The neural network (NN) models are well known for strong self-learning adaptability, high robustness, and excellent fault tolerance. It can reflect the implicit feature rules by the connection weights between neurons and solve many nonlinear problems including complex information, unclear rules, and noise pollution. The neural network models can evaluate the sensitivity values of model parameters in the situation that multiple parameters change simultaneously in complex models [33].In this paper, the BP neural network was employed to quantify the sensitivity of parameters of the damper force model. The weights and thresholds of the input layer to the hidden layer and of the hidden layer to the output layer were extracted through training the BP neural network with sample database. The BP neural network model and weights represented the mapping relationship between each parameter and output. The method of Garson (equation (5)) was applied to compute the sensitivity values of model parameters.(5)RI=∑j=1mwijwjk/∑i=1nwij∑i=1n∑j=1mwijwjk/∑i=1nwij,where Wij = the connection weight from input unit to hidden unit and Wjk = the connection weight of hidden unit to output unit. The sine excitation at amplitude 10 mm and frequency 4 Hz and input current 0.8 A were loaded on the MR damper to analyse parameters sensitivity.The main steps of sensitivity analysis method based on BP neural network are as follows:Step 1. Structure design of neural network: a three-layer neural network was created, including input layer, hidden layer, and output layer, to calculate the sensitivity values of each model parameter. There are six units in input layer: α, β, δ, F0, c0, and k0, respectively. The sensitivity values of parameters were calculated by equation (4). The number of units in the hidden layer was calculated by equation (6)(6)n1=n+m+a.where,n, n1, and m represent separately the number of input unit, hidden unit, and output unit; a = any constant between 1 and 10. After calculation, the number of units in the hidden layer was determined to be 10. The structure of neural network is illustrated in Figure 10.Step 2. Selection of training samples for neural network: six model parameters, α, β, δ, F0, c0, and k0, were selected as the training data, whose basic values, α = 1025.79, β = 0.01843, δ = 0.9715, F0 = −104.55N, c0 = 1.83, and k0 = 0.0025, were calculated under 10 mm, 4 Hz, and 0.8 A excitation condition. And, five levels for each parameter, 0, ±30%, and ±50% varied from basic values. By combination, 15625 (56) samples were obtained. Therfore, 15,000 samples were selected randomly as training data, and 625 samples from the remaining samples were picked out to verify the trained neural network model.Step 3. Analysis of the neural network model: by analysing the calculation results of the neural network, conclusions can be made that there are more than 90% samples whose errors are within 5%. Obviously, the mapping established through the neural network can get accurate relationship between model parameters and damper force output. The weights (Wij and Wjk) are listed in Table 2.Figure 10
The neural network structure.Table3 lists the sensitivity coefficients of each model parameter calculated by the methods of Garson (equation (14)). It is noteworthy that the parameter a has the greatest impact on the damper force output, with the sensitivity value of 39.43%. And, β, c0, and d rank behind, with values of 30.68%, 14.76%, and 14.71%, respectively. The remaining two parameters have only slight effect on the damper force output.Table 2
The weight of each neuron.
WijWjk10.7636770.773383−0.32296−0.001670.301491.02E−05−0.015021.2831210.2197590.129862−0.025440.205786−0.001621.28763−0.7358−0.141360.100185−1.36419−0.20055−0.00179−0.247240.934392−0.144−0.13751−0.021740.5765720.0017140.23885−4.48097−0.57274−0.00747−0.0086−1.428840.002235−0.10466−0.166821.768808−0.8945−0.00960.29989−0.000442.092271.833262−0.5778−0.07458−0.00660.476606−0.00023−47.63258−0.82338−0.690840.3101030.000858−0.309018.14E−0549.280090.869601−0.314370.3436720.0095630.285540.000773−8.873610−0.22652−1.6121.6863860.019954−0.192150.0090690.1586Table 3
The sensitivity results from the neural network model.
ParameterΑβδSensitivity (%)39.432530.683414.71123ParameterF0c0k0Sensitivity (%)0.394714.76000.0180The calculation result of sensitivity analysis show that the effects of parameterα, β, δ, and c on the model output are significant ones; meanwhile, the parameters F0 and k0 can be set to constants due to the low sensitivity.


## 5. The Shuffled Frog-Leaping Algorithm (SFLA)

To describe the force-velocity hysteretic characteristics of the MR damper, the parameters in the MHT model need to be identified.Parameter identification is a key issue in this research due to the complexity and nonlinearity of the model. In order to accurately determine the parameters, researchers have successively proposed a series of identification methods aimed at nonlinear models, such as genetic algorithm (GA), particle filter algorithm (PA), and least square method (LSC) [34].In this study, the processes of parameter identification utilizing shuffled frog-leaping algorithm and the experiment data fitting were carried out [35, 36]. The shuffled frog-leaping algorithm was firstly proposed in 2003 by Lansey and Eusuff [37] to solve the combinatorial optimization problem. This algorithm combines the advantages of particle swarm algorithm and memetic algorithm and is a new metaheuristic group evolution algorithm. It has been successfully applied in many fields especially in solving multiobjective complex function optimization problems. Figure 11 shows the parameter identification process of the model.Figure 11
The algorithm process of SFLA. (a) The flowchart of SFLA. (b) Algorithm steps.
(a)(b)Although there are some resemblances between SFLA and genetic algorithm (both are updated iteratively through fitness sorting), the SLFA has a faster spread speed. Memes in SFLA can flexibly use different mechanisms to spread information from one member of the memeplex to another, while genes in GA can only pass from parents to offspring. The evolution in genetic algorithm is limited by the number of offspring that a single parent can produce.Steps of the shuffled frog-leaping algorithm are as follows:Step 1. Initialize parameters and population: set t = 1; determine the population size F and randomly generate the initial population X (t); determine the number of subpopulations m and the number of individuals in each subpopulation n; F = m∗n.Step 2. Rank frogs: calculate the fitness values of each individual frog in the current population X (t).Step 3. Partial update iteration: t = t + 1. Sort the entire frog group according to their fitness and perform a local search in each subpopulation to update the worst frog.(1)
During the inner-group update process of the subpopulation, the algorithm only adjusts the worst frog’s position by computing(7)SW=MININTrandFPB−FPW,SWmax,for positive step=MAXINTrandFPB−FPW,−SWmax,for negative stepNL=Lw+SW.where SW and NL are the regeneration step size and the updated position of the worst frog, respectively; SWmax = the maximum step size of frog movement specifically; FPB = the best solution in perform; FPw = the worst solution in subpopulation.If the calculation results in step 3(1) meet the iteration requirements, move to the step 4; Otherwise, go to step 3(2).(2)
ReplaceFPBin equation (7) with the global optimal value FPXbest as follows:(8)SW=MININTrandFPbest−FPW,SWmax,for positive step=MAXINTrandFPbest−FPW,−SWmax,for negative stepNL=Lw+SW.(3)
If the two substeps above cannot meet the iteration requirements, the frog that is in the worst position in the subpopulation will jump to a random position (rp) in the solution space.(9)NL=rp.Step 4. The iteration ends after all the subpopulations have completed the partial search. Mix all the frogs together and regroup according to the updated fitness value.Step 5. Determine whether the current result meets the convergence demand. If not, go to step 2 and resort and divide the memeplex to perform the local search process. If it does, the process of algorithm ends.


## 6. Parameter Identification and Model Validation

### 6.1. Parameter Identification

The dynamic performance of the MR damper is not only related to input current but also influenced by the characteristics of the excitation (amplitude and frequency) according to the research in previous sections. The identified parameters extracted from a fixed working condition cannot accurately fit for other excitation conditions. A general model that is suitable for different excitation conditions should be built by establishing the relationship between model parameters and loading conditions including input current, vibration frequency, and amplitude.Although several studies have made effort to introduce the applied load to the damper model, most of these papers only study the influence of input current on the model output, without considering excitation characteristics [38]. Cheng and Chen [39] proposed a model containing amplitude and frequency.Noteworthy, the amplitude and frequency of vibration excitations cannot be measured in real time, especially under random excitation. Therefore, it is necessary to explore the relationship between excitation characteristics and measurable parameters during modeling. Wang and Ma [21] proposed a relational expression between movement states (displacement, velocity, and acceleration) of an MR damper and excitation characteristics. The maximum velocity (vm) of the damper has a linear relationship with the sine excitation amplitude and frequency. vm can thus be expressed by(10)vm=c0∗a⋅f,where c0 is a constant coefficient and a and f represent the amplitude and frequency, respectively.To verify the equation, the maximum velocity values under a series of working conditions were sorted out and then observed the relationship between thevm and amplitude or frequency (see Figure 12). It is evident that the vm amplitude curve and vm frequency curve both show linear relationships. As for vm variation under fixed amplitude (5 mm, 10 mm, 15 mm, and 20 mm), it increases with the increase of frequency straightly. Similarly, the linear increase relationship of vm versus amplitude is observed under each fixed frequency (1–5 Hz). Figure 13 demonstrates that the value of vm does not depend on the input of current approximately under each constant amplitude and frequency.Figure 12
The relationship between maximum velocity and excitation characteristic.Figure 13
The effect of current on maximum velocity under different excitations.Because the variablevm cannot be directly measured in real time either, and the damper instantaneous displacement and acceleration are connected with vm as follows:(11)vm=x˙2−x¨⋅x,where x˙ is the velocity of the damper andx¨ and x represent instantaneous displacement and acceleration, respectively. Therefore, the excitation frequency and amplitude are expressed by measurable dynamics variables of the damper (when the vehicle is operating in real time) by taking advantage of equations (10) and (11).Here, we firstly identified the parameters of the damper model from experiment data under a wide range of exciting conditions. The shuffled frog-leaping algorithm was utilized to obtain the accurate model parameters. And then, the correlation rules between each parameter and the excitation characteristics (current andvm) were excavated in purpose of building the general damper model.The relationships between parameterα and the intensity of the input current or variable vm are shown in Figure 14. It is clearly observed that the parameter a can be linearly represented with the input current (see Figure 14(a)). The quadratic polynomial has good fitting effect on the a-vm curve in Figure 14(b), described as(12)α=α1I+α2⋅α3vm2+α4vm+α5.Figure 14
Identification results of parameterα.
(a)(b)The relationships between parameterß and excitation characteristics or control current are shown in Figure 15. The results show that ß decreases nonlinearly as the excitation velocity amplitude increases. It can be found that the accuracy of the power function fitting is adequate. In the meanwhile, ß almost remains unchanged with the variation of control current. Therefore, ß can be expressed as follows:(13)β=β1⋅vmβ2.Figure 15
Identification results of parameterβ.
(a)(b)Noticeably, in equations (12) and (13), the coefficients α1, α2, α3, α4, α5, β1, and β2 are constants.Figures16(a) and 16(b) reveals that parameter δ depends on only the current value, rather than vibration characteristic vm. For fixed vibration input (vm), the parameter δ increases as the control current increases. A first-order polynomial can approximate the relationship between δ and current precisely.(14)δ=δ1I+δ2.Figure 16
Identification results of parametersδ and c0.
(a)(b)(c)(d)Observations from Figures16(c) and 16(d) are that the parameter c0 linearly rises as current increases, while the effect of vm on parameter c0 can be neglected. The parameter c0 is obtained as follows:(15)c0=c01I+c02.The variations of other parameters,k and F0, have nothing to do with the input current and excitation characteristics (Figure 17), and the sensitivity values of these parameters to the model output are quite low without exception. Therefore, it can be regarded as constants.Figure 17
Identification results of parametersF0 and k.
(a)(b)Plugged the calculated results in the hyperbolic tangent model, equations are shown as follows:(16)F=αz+F0+c0x˙+k0x,z=tanhβx˙+δsignx,vm=x˙2−x¨⋅x,α=α1I+α2⋅α3vm2+α4vm+α5,β=β1⋅vmβ2,δ=δ1I+δ2,c0=c01I+c02.The number of parameters to be identified increases from 6 (α, β, δ, F0, c0, and k0) to 13 (a1, a2, a3, a4, a5, ß1, ß2, d1, d2, F0, c01, c02, and k0). The parameters were reidentified under every input current and vibration condition, using the shuffled frog-leaping algorithm. The recalculated parameters are shown in Table 4.Table 4
Summary of identification results.
ParameterValueParameterValueα11.4478δ20.4354α20.1233F0−88.7903α3−0.0002c011.4962α40.0608c010.9767α5702.0492k00.0208β12.9106β2−0.8621δ10.5472The hyperbolic tangent model only considering current variable was established through the same process as above:(17)z=tanh0.09276x˙+0.361signx,F=cx˙+0.1897x+αz+F0,α=938.1199I+58.0515,δ=0.6160I+0.4273,c0=2.3285I+0.9799.Figure18 shows velocity-damping force characteristics of different models under 5 mm, 2 Hz sine excitation. The black curve represents the damping force measured by the MR damper mechanical test, the red curve represents the hyperbolic tangent model only considering the input currents, and the blue curve represents the hyperbolic tangent model considering both excitations and currents.Figure 18
Damper force calculation results, hyperbolic tangent model, at 2 Hz excitation with amplitude of 5 mm. (a) Input current = 0.4 A. (b) Input current = 0.6 A. (c) Input current = 0.8 A. (d) Input current = 1 A.
(a)(b)(c)(d)We can notice that the damper force calculated with the hyperbolic tangent model which only considers current variable has a large margin of error whether in low or high velocity area. The reason leading to great error is that the hysteretic model parameters do not include the elements related to excitation characteristics, while the parametersα, β, and c0 are significantly influenced by the excitation frequency and amplitude as well.The damper force calculated by the detailed model including input currents and excitation conditions is more in line with the experimental curve under different working conditions. The output damping force error RMS of the two models compared to the test damping force is 83.7974 N and 65.3787 N, respectively, under 0.4 A current input. The error RMS of the two models is 205.77 N and 137.75 N under 1 A current input. The error RMS of modified model are reduced by 21.98%, 24.77%, 29.33%, and 33.06% compared with the hyperbolic tangent model only considering current. Therefore, it can be found that the improvement effect of the model is better as the input current increases.


### 6.2. Model Validation

A two-degree-of-freedom seat suspension model is established to verify the damping control of the proposed model applied to the semiactive suspension. The model structure is shown in Figure19.Figure 19
Semiactive seat suspension system model.The motion differential equation of the seat suspension model shown in Figure19 is shown as follows:(18)Mpz¨p+kpzp−zs+cpz˙p−z˙s=0,Msz¨s+kszs−z0+Fd−kpzp−zs−cpz˙p−z˙s=0,Msz¨s+kszs−z0+cpz˙s−z˙0−kpzp−zs−cpz˙p−z˙s=0.The parameter meanings in equation (18) are listed in Table 5.Table 5
Parameter meaning of seat suspension.
ParameterRepresentationValuesMpPassenger weight50 kgMsSeat weight17 kgkpCushion stiffness19556 N/mksSuspension stiffness78000 N/mcpCushion damping1865 Nm/scsPassive suspension damping942 Nm/sFdMR damper force—zpPassenger displacement—zsSeat displacement—z0Excitation displacement—The sky-hook control strategy is used for semiactive control simulation. The damping coefficient of sky-hook control was selected as 2000 Nm/s. The control current of MR damper was calculated by the lookup table method. The database of lookup table was obtained from the established forward model. Figure20 shows the simulation results under random excitation. The included frequency range is 1–10 Hz. The acceleration of the seat upper surface is reduced when the semiactive control is applied. When the modified MR damper model is applied to the control system, the acceleration is further reduced, as the black curve in the figure shows.Figure 20
The acceleration contrast of seat upper surface.In order to see the vibration reduction effect more clearly, the vibration transmissibility under different frequencies is shown in Figure21. Semiactive suspension of the MR damper can effectively reduce the transmission of vibration, especially at frequencies sensitive to the human body (4–8 Hz). At the resonance peak frequency, the original tanh model reduces the transmissibility by 0.79, and the modified model further reduces the transmissibility by 0.18.Figure 21
Comparison of vibration transmissibility of different models.Figure22 shows the comparison of the control current of the two models. In order to clearly see the difference, the sampling time of the current is set as 0.1 s. It can be seen that the output current of the two models is different. According to Figures 20–22, it can be determined that the control current calculated by the modified model is more accurate. The accurate model can further improve the vibration reduction effect of the whole semiactive suspension control system.Figure 22
Control current for two kinds of model.


## 6.1. Parameter Identification

The dynamic performance of the MR damper is not only related to input current but also influenced by the characteristics of the excitation (amplitude and frequency) according to the research in previous sections. The identified parameters extracted from a fixed working condition cannot accurately fit for other excitation conditions. A general model that is suitable for different excitation conditions should be built by establishing the relationship between model parameters and loading conditions including input current, vibration frequency, and amplitude.Although several studies have made effort to introduce the applied load to the damper model, most of these papers only study the influence of input current on the model output, without considering excitation characteristics [38]. Cheng and Chen [39] proposed a model containing amplitude and frequency.Noteworthy, the amplitude and frequency of vibration excitations cannot be measured in real time, especially under random excitation. Therefore, it is necessary to explore the relationship between excitation characteristics and measurable parameters during modeling. Wang and Ma [21] proposed a relational expression between movement states (displacement, velocity, and acceleration) of an MR damper and excitation characteristics. The maximum velocity (vm) of the damper has a linear relationship with the sine excitation amplitude and frequency. vm can thus be expressed by(10)vm=c0∗a⋅f,where c0 is a constant coefficient and a and f represent the amplitude and frequency, respectively.To verify the equation, the maximum velocity values under a series of working conditions were sorted out and then observed the relationship between thevm and amplitude or frequency (see Figure 12). It is evident that the vm amplitude curve and vm frequency curve both show linear relationships. As for vm variation under fixed amplitude (5 mm, 10 mm, 15 mm, and 20 mm), it increases with the increase of frequency straightly. Similarly, the linear increase relationship of vm versus amplitude is observed under each fixed frequency (1–5 Hz). Figure 13 demonstrates that the value of vm does not depend on the input of current approximately under each constant amplitude and frequency.Figure 12
The relationship between maximum velocity and excitation characteristic.Figure 13
The effect of current on maximum velocity under different excitations.Because the variablevm cannot be directly measured in real time either, and the damper instantaneous displacement and acceleration are connected with vm as follows:(11)vm=x˙2−x¨⋅x,where x˙ is the velocity of the damper andx¨ and x represent instantaneous displacement and acceleration, respectively. Therefore, the excitation frequency and amplitude are expressed by measurable dynamics variables of the damper (when the vehicle is operating in real time) by taking advantage of equations (10) and (11).Here, we firstly identified the parameters of the damper model from experiment data under a wide range of exciting conditions. The shuffled frog-leaping algorithm was utilized to obtain the accurate model parameters. And then, the correlation rules between each parameter and the excitation characteristics (current andvm) were excavated in purpose of building the general damper model.The relationships between parameterα and the intensity of the input current or variable vm are shown in Figure 14. It is clearly observed that the parameter a can be linearly represented with the input current (see Figure 14(a)). The quadratic polynomial has good fitting effect on the a-vm curve in Figure 14(b), described as(12)α=α1I+α2⋅α3vm2+α4vm+α5.Figure 14
Identification results of parameterα.
(a)(b)The relationships between parameterß and excitation characteristics or control current are shown in Figure 15. The results show that ß decreases nonlinearly as the excitation velocity amplitude increases. It can be found that the accuracy of the power function fitting is adequate. In the meanwhile, ß almost remains unchanged with the variation of control current. Therefore, ß can be expressed as follows:(13)β=β1⋅vmβ2.Figure 15
Identification results of parameterβ.
(a)(b)Noticeably, in equations (12) and (13), the coefficients α1, α2, α3, α4, α5, β1, and β2 are constants.Figures16(a) and 16(b) reveals that parameter δ depends on only the current value, rather than vibration characteristic vm. For fixed vibration input (vm), the parameter δ increases as the control current increases. A first-order polynomial can approximate the relationship between δ and current precisely.(14)δ=δ1I+δ2.Figure 16
Identification results of parametersδ and c0.
(a)(b)(c)(d)Observations from Figures16(c) and 16(d) are that the parameter c0 linearly rises as current increases, while the effect of vm on parameter c0 can be neglected. The parameter c0 is obtained as follows:(15)c0=c01I+c02.The variations of other parameters,k and F0, have nothing to do with the input current and excitation characteristics (Figure 17), and the sensitivity values of these parameters to the model output are quite low without exception. Therefore, it can be regarded as constants.Figure 17
Identification results of parametersF0 and k.
(a)(b)Plugged the calculated results in the hyperbolic tangent model, equations are shown as follows:(16)F=αz+F0+c0x˙+k0x,z=tanhβx˙+δsignx,vm=x˙2−x¨⋅x,α=α1I+α2⋅α3vm2+α4vm+α5,β=β1⋅vmβ2,δ=δ1I+δ2,c0=c01I+c02.The number of parameters to be identified increases from 6 (α, β, δ, F0, c0, and k0) to 13 (a1, a2, a3, a4, a5, ß1, ß2, d1, d2, F0, c01, c02, and k0). The parameters were reidentified under every input current and vibration condition, using the shuffled frog-leaping algorithm. The recalculated parameters are shown in Table 4.Table 4
Summary of identification results.
ParameterValueParameterValueα11.4478δ20.4354α20.1233F0−88.7903α3−0.0002c011.4962α40.0608c010.9767α5702.0492k00.0208β12.9106β2−0.8621δ10.5472The hyperbolic tangent model only considering current variable was established through the same process as above:(17)z=tanh0.09276x˙+0.361signx,F=cx˙+0.1897x+αz+F0,α=938.1199I+58.0515,δ=0.6160I+0.4273,c0=2.3285I+0.9799.Figure18 shows velocity-damping force characteristics of different models under 5 mm, 2 Hz sine excitation. The black curve represents the damping force measured by the MR damper mechanical test, the red curve represents the hyperbolic tangent model only considering the input currents, and the blue curve represents the hyperbolic tangent model considering both excitations and currents.Figure 18
Damper force calculation results, hyperbolic tangent model, at 2 Hz excitation with amplitude of 5 mm. (a) Input current = 0.4 A. (b) Input current = 0.6 A. (c) Input current = 0.8 A. (d) Input current = 1 A.
(a)(b)(c)(d)We can notice that the damper force calculated with the hyperbolic tangent model which only considers current variable has a large margin of error whether in low or high velocity area. The reason leading to great error is that the hysteretic model parameters do not include the elements related to excitation characteristics, while the parametersα, β, and c0 are significantly influenced by the excitation frequency and amplitude as well.The damper force calculated by the detailed model including input currents and excitation conditions is more in line with the experimental curve under different working conditions. The output damping force error RMS of the two models compared to the test damping force is 83.7974 N and 65.3787 N, respectively, under 0.4 A current input. The error RMS of the two models is 205.77 N and 137.75 N under 1 A current input. The error RMS of modified model are reduced by 21.98%, 24.77%, 29.33%, and 33.06% compared with the hyperbolic tangent model only considering current. Therefore, it can be found that the improvement effect of the model is better as the input current increases.


## 6.2. Model Validation

A two-degree-of-freedom seat suspension model is established to verify the damping control of the proposed model applied to the semiactive suspension. The model structure is shown in Figure19.Figure 19
Semiactive seat suspension system model.The motion differential equation of the seat suspension model shown in Figure19 is shown as follows:(18)Mpz¨p+kpzp−zs+cpz˙p−z˙s=0,Msz¨s+kszs−z0+Fd−kpzp−zs−cpz˙p−z˙s=0,Msz¨s+kszs−z0+cpz˙s−z˙0−kpzp−zs−cpz˙p−z˙s=0.The parameter meanings in equation (18) are listed in Table 5.Table 5
Parameter meaning of seat suspension.
ParameterRepresentationValuesMpPassenger weight50 kgMsSeat weight17 kgkpCushion stiffness19556 N/mksSuspension stiffness78000 N/mcpCushion damping1865 Nm/scsPassive suspension damping942 Nm/sFdMR damper force—zpPassenger displacement—zsSeat displacement—z0Excitation displacement—The sky-hook control strategy is used for semiactive control simulation. The damping coefficient of sky-hook control was selected as 2000 Nm/s. The control current of MR damper was calculated by the lookup table method. The database of lookup table was obtained from the established forward model. Figure20 shows the simulation results under random excitation. The included frequency range is 1–10 Hz. The acceleration of the seat upper surface is reduced when the semiactive control is applied. When the modified MR damper model is applied to the control system, the acceleration is further reduced, as the black curve in the figure shows.Figure 20
The acceleration contrast of seat upper surface.In order to see the vibration reduction effect more clearly, the vibration transmissibility under different frequencies is shown in Figure21. Semiactive suspension of the MR damper can effectively reduce the transmission of vibration, especially at frequencies sensitive to the human body (4–8 Hz). At the resonance peak frequency, the original tanh model reduces the transmissibility by 0.79, and the modified model further reduces the transmissibility by 0.18.Figure 21
Comparison of vibration transmissibility of different models.Figure22 shows the comparison of the control current of the two models. In order to clearly see the difference, the sampling time of the current is set as 0.1 s. It can be seen that the output current of the two models is different. According to Figures 20–22, it can be determined that the control current calculated by the modified model is more accurate. The accurate model can further improve the vibration reduction effect of the whole semiactive suspension control system.Figure 22
Control current for two kinds of model.


## 7. Conclusions

This paper proposed a detailed, accurate damper model, which was based on the hyperbolic tangent model, to solve the problems of nonlinearity and hysteresis characteristics of an MR damper.By employing the mechanical test of the MR damper, the fitting effect of the different hysteresis model was verified. The hyperbolic tangent model was selected as the research object.The calculation result of sensitivity method was analysed comprehensively to confirm the parameters that have important impacts on the output damper force of the hyperbolic tangent hysteretic model, in order to provide a guidance for the design of the extended model. It is proved that the neural network model is appropriate to analyse parameter sensitivity.After analysing and comparing different hysteretic model, we identified the model parameters using the shuffled frog-leaping algorithm. By establishing the relationship between sensitive parameters and loading characteristics (input current, excitation amplitude, and frequency), the modified model and its parameters were determined. The results show that the modified model is an efficient hysteretic model to characterize the damping property of the MR damper under a wide range of control current and excitation conditions. The precise model can greatly improve the damping effect of the whole semiactive suspension control system.



---

*Source: 6691650-2021-04-08.xml*