Optimized Fuzzy Skyhook Control for Semi-Active Vehicle Suspension with New Inverse Model of Magnetorheological Fluid Damper

: To improve the performance of vehicle suspension, this paper proposes a semi-active vehicle suspension with a magnetorheological ﬂuid (MRF) damper. We designed an optimized fuzzy skyhook controller with grey wolf optimizer (GWO) algorithm base on a new neuro-inverse model of the MRF damper. Because the inverse model of the MRF damper is difﬁcult to establish directly, the Elman neural network was applied. The novelty of this study is the application of the new inverse model for semi-active vibration control and optimization of the semi-active suspension control method. The calculation results showed that the new inverse model can accurately calculate the required control current. The fuzzy skyhook control method optimized by the grey wolf optimizer (GWO) algorithm was established based on the inverse model to control the suspension vibration. The simulation results showed that the optimized fuzzy skyhook control method can simultaneously reduce the amplitude of vertical acceleration, suspension deﬂection, and tire dynamic load.


Introduction
Magnetorheological fluid (MRF) is an intelligent material that can transform between a Newtonian fluid state and a semi-solid state according to the strength of the magnetic field. The ferromagnetic particles of MRF are randomly distributed without an external magnetic field, so MRF performs as a Newtonian fluid and can flow under only small shear. Furthermore, the ferromagnetic particles in MRF rearrange along the direction of the magnetic field. The higher the strength of the magnetic field, the denser the arrangement of the ferromagnetic particles. Finally, MRF exhibits semi-solid properties, and much greater shear force is required to make it flow. When the applied magnetic field disappears suddenly, the MRF material can return to its original state in milliseconds [1,2].
Due to the fact of its variable and controllable properties, an MRF damper is considered to have great potential in vibration mitigation. The damping force of an MRF damper can vary with the applied current or voltage due to the change in the magnetic field strength [3,4].
For vehicle vibration reduction, to improve the performance of vehicle suspension, the main schemes include active and semi-active suspensions. Although active suspension better reduces vibrations, it also has some disadvantages, such as complex structure, high cost, and high energy consumption, which limit its application in mass production vehicles [5]. With the help of a control algorithm, semi-active suspension represented by an MRF damper can achieve the desirable effect of active control without high-cost hardware.
When the MRF damper is applied to vehicle suspension, a precise mechanical model of the damping force is essential for vibration control. The mechanical model of the damper. Therefore, the fuzzy skyhook control method was selected to control the MRF suspension system.
In the fuzzy control method, the two control inputs must be multiplied by quantitative factors, K a and K v , then calculated using fuzzy rules. The calculation result of the fuzzy controller must be multiplied by a scale factor K u and then output as the control result.
To perform well, fuzzy controllers must have a reasonable choice of quantitative factors K a , K v , and scale factor K u . These parameters have a certain regulating effect on the control performance. Increasing the quantization factor, K a , can speed up the rising rate of the system, but can also lead to an increase in overshoot. Increasing the quantization factor, K v , can improve the stability of the system. However, if K v is too large, the system rise rate process is too slow. Increasing the scale factor K u will make the system respond faster, but too large a value of K u can cause the system to overshoot and undergo oscillation dispersion. However, these three parameters are generally selected based on experience, so it is difficult to obtain the best control effect. The grey wolf algorithm is a new intelligent optimization algorithm with the characteristics of a simple principle, few adjustment parameters, and low computational requirement, and it has been widely used in the field of parameter optimization [22]. Therefore, the grey wolf optimizer algorithm was used in this study to optimize the parameters, or quantification factors K a and K v , and the scale factor K u in the fuzzy control to achieve the best control effect.
The research route of this paper can be summarized as follows: Because the semiactive control algorithm is used for the MRF damper, the MRF damper control force output is determined by the input current. Therefore, it is necessary to establish a reliable method to control the current through the output damping force of the MRF damper. Because of the complexity of the dynamic model of the MRF damper, the neural network method was used. At present, the BP neural network is usually used to establish the inverse model, and the accuracy of the model needs to be improved. Therefore, in this study, an Elman neural network was used to establish the inverse model of the MRF damper, which can obtain a control current according to the control force that is calculated by the control algorithm.
Regarding semi-active control, the fuzzy skyhook control method is a suitable control method for the MRF damper. However, the selection of key parameters K a , K v , and K u in fuzzy control is not accurate and mainly depends on experience. Therefore, to address this problem, this study optimized fuzzy control to obtain a better suspension damping effect. The schematic of the research workflow is shown in Figure 1. force can only be the same as the direction of the damping force produced by the MRF damper. The output law of the fuzzy skyhook control method coincides with the characteristics of the damper. Therefore, the fuzzy skyhook control method was selected to control the MRF suspension system. In the fuzzy control method, the two control inputs must be multiplied by quantitative factors, Ka and Kv, then calculated using fuzzy rules. The calculation result of the fuzzy controller must be multiplied by a scale factor Ku and then output as the control result.
To perform well, fuzzy controllers must have a reasonable choice of quantitative factors Ka, Kv, and scale factor Ku. These parameters have a certain regulating effect on the control performance. Increasing the quantization factor, Ka, can speed up the rising rate of the system, but can also lead to an increase in overshoot. Increasing the quantization factor, Kv, can improve the stability of the system. However, if Kv is too large, the system rise rate process is too slow. Increasing the scale factor Ku will make the system respond faster, but too large a value of Ku can cause the system to overshoot and undergo oscillation dispersion. However, these three parameters are generally selected based on experience, so it is difficult to obtain the best control effect. The grey wolf algorithm is a new intelligent optimization algorithm with the characteristics of a simple principle, few adjustment parameters, and low computational requirement, and it has been widely used in the field of parameter optimization [22]. Therefore, the grey wolf optimizer algorithm was used in this study to optimize the parameters, or quantification factors Ka and Kv, and the scale factor Ku in the fuzzy control to achieve the best control effect.
The research route of this paper can be summarized as follows: Because the semiactive control algorithm is used for the MRF damper, the MRF damper control force output is determined by the input current. Therefore, it is necessary to establish a reliable method to control the current through the output damping force of the MRF damper. Because of the complexity of the dynamic model of the MRF damper, the neural network method was used. At present, the BP neural network is usually used to establish the inverse model, and the accuracy of the model needs to be improved. Therefore, in this study, an Elman neural network was used to establish the inverse model of the MRF damper, which can obtain a control current according to the control force that is calculated by the control algorithm.
Regarding semi-active control, the fuzzy skyhook control method is a suitable control method for the MRF damper. However, the selection of key parameters Ka, Kv, and Ku in fuzzy control is not accurate and mainly depends on experience. Therefore, to address this problem, this study optimized fuzzy control to obtain a better suspension damping effect. The schematic of the research workflow is shown in Figure 1.   The main work of this paper is shown as follows: the mechanical properties of the MRF damper were tested, and the dynamic model of the MRF damper was established based on the Bouc-Wen model. To obtain accurate results of the control current, the inverse model of the MRF damper was modeled and trained using the Elman neural network. Finally, the neural network inverse model was applied to the semi-active control of the quarter vehicle suspension system. The fuzzy skyhook method was used to control the suspension vibration. Aiming to overcome the shortcomings of the fuzzy control approach, the GWO algorithm was used to optimize the quantization factors and the scale factor in the fuzzy controller to obtain a satisfactory control effect. The analysis results indicate that the dynamic performance of the vehicle suspension improved, and the inverse model of the MRF damper proposed in this paper was validated.

Mechanical Test of the MRF Damper
The structure of the MRF damper is shown in Figure 2. It is composed of a cylinder block, piston, floating piston, magnetorheological fluid, coil, sealing ring, gas chamber, and coil lead wire. There is an orifice on the piston and a coil wound around the piston to generate a strong magnetic field. The main work of this paper is shown as follows: the mechanical properties of the MRF damper were tested, and the dynamic model of the MRF damper was established based on the Bouc-Wen model. To obtain accurate results of the control current, the inverse model of the MRF damper was modeled and trained using the Elman neural network. Finally, the neural network inverse model was applied to the semi-active control of the quarter vehicle suspension system. The fuzzy skyhook method was used to control the suspension vibration. Aiming to overcome the shortcomings of the fuzzy control approach, the GWO algorithm was used to optimize the quantization factors and the scale factor in the fuzzy controller to obtain a satisfactory control effect. The analysis results indicate that the dynamic performance of the vehicle suspension improved, and the inverse model of the MRF damper proposed in this paper was validated.

Mechanical Test of the MRF Damper
The structure of the MRF damper is shown in Figure 2. It is composed of a cylinder block, piston, floating piston, magnetorheological fluid, coil, sealing ring, gas chamber, and coil lead wire. There is an orifice on the piston and a coil wound around the piston to generate a strong magnetic field. To obtain the damping characteristics of the MRF damper, the magnetorheological fluid damper was installed on the Instron E10000 electronic dynamic testing machine (manufactured by Instron company in Boston, US) as shown in Figure 3. In the test, the excitation frequency and amplitude could set to make the MRF damper piston move up and down under the action of the dynamic testing machine, and the input current of the MRF damper was controlled by a DC power supply. In this experiment, the set current was 0, 0.25, 0.5, 0.75, and 1 A; the displacement input excitation was 5 and 10 mm; the frequency was 0.5, 1, and 2 Hz, for a total of 30 working conditions. Due to the large number of working conditions, we used the 5 mm, 1 Hz working condition as an example for analysis as shown in Figure 4. To obtain the damping characteristics of the MRF damper, the magnetorheological fluid damper was installed on the Instron E10000 electronic dynamic testing machine (manufactured by Instron company in Boston, US) as shown in Figure 3. In the test, the excitation frequency and amplitude could set to make the MRF damper piston move up and down under the action of the dynamic testing machine, and the input current of the MRF damper was controlled by a DC power supply. In this experiment, the set current was 0, 0.25, 0.5, 0.75, and 1 A; the displacement input excitation was 5 and 10 mm; the frequency was 0.5, 1, and 2 Hz, for a total of 30 working conditions. Due to the large number of working conditions, we used the 5 mm, 1 Hz working condition as an example for analysis as shown in Figure 4. Figure 4 shows that the shape of the damping force and displacement curve of the MRF damper was rectangular and full, which indicates that the damper has good vibration attenuation ability. Under the same amplitude and frequency excitation, the higher the input current, the greater the output damping force of the MRF damper. The forcevelocity curve of the MRF damper was approximately a symmetric hyperbola with obvious hysteresis characteristics. When the relative velocity of the MRF damper piston was small, the MRF damper presented the hysteresis phenomenon, and the area of the hysteresis loop increased with the increase in the input current; when the relative velocity of the damper piston was large, the curve shows shear phenomenon, revealing that the MRF damper has strong nonlinear dynamic characteristics.    Figure 4 shows that the shape of the damping force and displacement curve of the MRF damper was rectangular and full, which indicates that the damper has good vibration attenuation ability. Under the same amplitude and frequency excitation, the higher the input current, the greater the output damping force of the MRF damper. The forcevelocity curve of the MRF damper was approximately a symmetric hyperbola with obvious hysteresis characteristics. When the relative velocity of the MRF damper piston was small, the MRF damper presented the hysteresis phenomenon, and the area of the hysteresis loop increased with the increase in the input current; when the relative velocity of the damper piston was large, the curve shows shear phenomenon, revealing that the MRF damper has strong nonlinear dynamic characteristics.

Bouc-Wen Model
To fully utilize the semi-active characteristic of the MRF damper, we need to build an accurate dynamic model to achieve the real-time performance of semi-active control. Among various models, the Bouc-Wen model is the most commonly used [6,7], and describes the hyperbolic hysteresis characteristics of the MRF damper. It is made of a linear spring, a linear damper, and a hysteretic element in parallel. As shown in Figure 5, the      Figure 4 shows that the shape of the damping force and displacement curve of the MRF damper was rectangular and full, which indicates that the damper has good vibration attenuation ability. Under the same amplitude and frequency excitation, the higher the input current, the greater the output damping force of the MRF damper. The forcevelocity curve of the MRF damper was approximately a symmetric hyperbola with obvious hysteresis characteristics. When the relative velocity of the MRF damper piston was small, the MRF damper presented the hysteresis phenomenon, and the area of the hysteresis loop increased with the increase in the input current; when the relative velocity of the damper piston was large, the curve shows shear phenomenon, revealing that the MRF damper has strong nonlinear dynamic characteristics.

Bouc-Wen Model
To fully utilize the semi-active characteristic of the MRF damper, we need to build an accurate dynamic model to achieve the real-time performance of semi-active control. Among various models, the Bouc-Wen model is the most commonly used [6,7], and describes the hyperbolic hysteresis characteristics of the MRF damper. It is made of a linear spring, a linear damper, and a hysteretic element in parallel. As shown in Figure 5, the

Bouc-Wen Model
To fully utilize the semi-active characteristic of the MRF damper, we need to build an accurate dynamic model to achieve the real-time performance of semi-active control. Among various models, the Bouc-Wen model is the most commonly used [6,7], and describes the hyperbolic hysteresis characteristics of the MRF damper. It is made of a linear spring, a linear damper, and a hysteretic element in parallel. As shown in Figure 5, the hysteresis curve function was added to describe the hysteretic characteristics of the MRF damper. The damping force of the MRF damper is shown as follows: where F is the damping force, c 0 is the viscosity coefficient, k 0 is the stiffness coefficient, α is the proportion adjustment parameter of the hysteretic force in the damping force, x 0 is the relative equilibrium displacement, z is the hysteretic variable, x and .
x are the relative displacement and relative velocity at both ends of the damping, respectively, n is the smoothing coefficient, and γ, β, and A are parameters used to adjust the smoothness of the relative equilibrium displacement, z is the hysteretic variable, x and x  are the relative displacement and relative velocity at both ends of the damping, respectively, n is the smoothing coefficient, and γ, β, and A are parameters used to adjust the smoothness of the model and the linear characteristics from the yield area to the post-yield area of the curve.
Because this model can better describe the relationship between the output damping force and the piston speed, it can also better reflect the hysteresis characteristics of the magnetorheological fluid damper, so we selected this method to fit the MRF damper dynamic model.

The Dynamic Model of MRF Damper
The eight unknown parameters of the Bouc-Wen model were difficult to identify due to the relationship between the Bouc-Wen model and the current. Using the Simulink Design Optimization toolbox [23], the measured displacement and velocity were taken as inputs to calculate the damping force using the landmark programming of interconnection by adjusting the vector of unknown parameters. Then, the calculated result was compared with the test results.
The parameter identification steps were as follows: (1) Build the Bouc-Wen model in Simulink. The model is shown in Figure 6.
(2) Import the test data. The test data of sinusoidal excitation with 5 mm, 1 Hz shown in Figure 3 under different current conditions were imported into the Simulink model.
The input values of the model were displacement and velocity, and the output value was damping force. (3) Select the parameters to identify. The Bouc-Wen model has eight unknown parameters; first, set an initial value first, and then start to identify. The initial values are listed in Table 1.  Because this model can better describe the relationship between the output damping force and the piston speed, it can also better reflect the hysteresis characteristics of the magnetorheological fluid damper, so we selected this method to fit the MRF damper dynamic model.

The Dynamic Model of MRF Damper
The eight unknown parameters of the Bouc-Wen model were difficult to identify due to the relationship between the Bouc-Wen model and the current. Using the Simulink Design Optimization toolbox [23], the measured displacement and velocity were taken as inputs to calculate the damping force using the landmark programming of interconnection by adjusting the vector of unknown parameters. Then, the calculated result was compared with the test results.
The parameter identification steps were as follows: (1) Build the Bouc-Wen model in Simulink. The model is shown in Figure 6.
(2) Import the test data. The test data of sinusoidal excitation with 5 mm, 1 Hz shown in Figure 3 under different current conditions were imported into the Simulink model.
The input values of the model were displacement and velocity, and the output value was damping force. (3) Select the parameters to identify. The Bouc-Wen model has eight unknown parameters; first, set an initial value first, and then start to identify. The initial values are listed in Table 1. of the transition curve from the elastic zone to the plastic zone [8]. From Figure 4, it can be seen that the magnitude of the current has a great influence on the damping force amplitude. Through the preliminary identification of the model parameters, we found that the values of parameters A, α, and c0 vary greatly with the current, so it was assumed that a relationship can be developed between A, α, c0 and the magnitude of current. The other parameters were fixed as constants. The parameters identified are listed in Table 2.  After the parameter recognition of the Bouc-Wen model was completed, the output damping force can be calculated by inputting the displacement, velocity, and current of the MRF damper, and then the displacement-force curve and velocity-force curve can be simulated.
To verify the model identification results, the simulation results were compared with the test data under different currents as shown in Figure 7. In the Figure 7, the experimental data are shown as solid lines and the simulation results are dotted lines. In the Figure 7, we can see that the calculated displacement-force curve and the velocity-force (4) In the process of parameter identification, the parameters have different effects on the Bouc-Wen model. It was found that γ and β only affect the shape of the hysteresis loop but not the magnitude of the damping force, and n only affects the smoothness of the transition curve from the elastic zone to the plastic zone [8]. From Figure 4, it can be seen that the magnitude of the current has a great influence on the damping force amplitude. Through the preliminary identification of the model parameters, we found that the values of parameters A, α, and c 0 vary greatly with the current, so it was assumed that a relationship can be developed between A, α, c 0 and the magnitude of current. The other parameters were fixed as constants. The parameters identified are listed in Table 2. After the parameter recognition of the Bouc-Wen model was completed, the output damping force can be calculated by inputting the displacement, velocity, and current of the MRF damper, and then the displacement-force curve and velocity-force curve can be simulated.
To verify the model identification results, the simulation results were compared with the test data under different currents as shown in Figure 7. In the Figure   The parameters of the Bouc-Wen model were obtained under a specific case of 1 Hz sinusoidal excitation of 5 mm amplitude, so the generality of the model (i.e., the dynamic response under other excitations) was also tested. The working condition of 0.5 Hz, 10 mm with 0.25 A current was selected for analysis. The sinusoidal excitation signal of this working condition was input into the Bouc-Wen model identified in Table 2 to calculate the damping force. The displacement damping force current and velocity damping force current were determined, then compared with the test results as shown in Figure 8. From Figure 8, we can see that the proposed model can accurately reflect the dynamic response of the MRF damper under sinusoidal excitation in other cases, it has good universality.

The Inverse Dynamic Model of MRF Damper
In the semi-active control, the control force can be calculated by the control algorithm, which needs to be output by the MRF damper. At this time, it is necessary to regulate the current of the MRF damper to change the output damping force. Therefore, it is necessary The parameters of the Bouc-Wen model were obtained under a specific case of 1 Hz sinusoidal excitation of 5 mm amplitude, so the generality of the model (i.e., the dynamic response under other excitations) was also tested. The working condition of 0.5 Hz, 10 mm with 0.25 A current was selected for analysis. The sinusoidal excitation signal of this working condition was input into the Bouc-Wen model identified in Table 2 to calculate the damping force. The displacement damping force current and velocity damping force current were determined, then compared with the test results as shown in Figure 8. The parameters of the Bouc-Wen model were obtained under a specific case of 1 Hz sinusoidal excitation of 5 mm amplitude, so the generality of the model (i.e., the dynamic response under other excitations) was also tested. The working condition of 0.5 Hz, 10 mm with 0.25 A current was selected for analysis. The sinusoidal excitation signal of this working condition was input into the Bouc-Wen model identified in Table 2 to calculate the damping force. The displacement damping force current and velocity damping force current were determined, then compared with the test results as shown in Figure 8. From Figure 8, we can see that the proposed model can accurately reflect the dynamic response of the MRF damper under sinusoidal excitation in other cases, it has good universality.

The Inverse Dynamic Model of MRF Damper
In the semi-active control, the control force can be calculated by the control algorithm, which needs to be output by the MRF damper. At this time, it is necessary to regulate the current of the MRF damper to change the output damping force. Therefore, it is necessary From Figure 8, we can see that the proposed model can accurately reflect the dynamic response of the MRF damper under sinusoidal excitation in other cases, it has good universality.

The Inverse Dynamic Model of MRF Damper
In the semi-active control, the control force can be calculated by the control algorithm, which needs to be output by the MRF damper. At this time, it is necessary to regulate the current of the MRF damper to change the output damping force. Therefore, it is necessary to build an inverse dynamic model of the MRF damper to calculate the control current from the expected damping force.
Because the relationship between the output force and the control current of the MRF damper is complex, the mathematical expression of the relationship between them Energies 2021, 14, 1674 9 of 21 cannot be calculated directly. Thus, the neural network method was used to establish the relationship between them to predict the required control current according to the expected damping force. In this study, the Elman neural network was used to establish the inverse model of the MRF damper.

Data Generation
To ensure the training samples can cover the working condition of the MRF damper, the following principals were used: (1) the current range of training samples should include the application range from 0 to 1 A; (2) the samples of displacement and speed range should cover a wide range of applications; and (3) the number of samples should be large enough.
Based on the above principles, because this MRF damper was mainly used on car suspensions, to make the training data relevant to practical applications, a quarter vehicle suspension model was established. The suspension displacement calculated from the road spectrum and its differential result suspension velocity were used as the training data of the neural network. The current data were set as a random white noise signal with a magnitude from 0 to 1 A. The output damping force was obtained using the forward model of the MRF damper established above.

Elman Neural Network
To realize the real-time control of the MRF damper, the structure of the neural network should be as simple as possible. The three inputs were the piston displacement x, the piston velocity v, and the output damping force F. The output variable was the control current I.
The Elman neural network is a multilayer feedforward neural network and is based on the BP neural network and negative feedback system. BP neural networks generally have an input layer, several hidden layers, and an output layer. In addition, the Elman neural network adds an undertaking layer. Compared with BP neural networks, the Elman neural network has better stability and is better at classifying the original data. The model structure is shown in Figure 9 [24].
to build an inverse dynamic model of the MRF damper to calculate the control current from the expected damping force.
Because the relationship between the output force and the control current of the MRF damper is complex, the mathematical expression of the relationship between them cannot be calculated directly. Thus, the neural network method was used to establish the relationship between them to predict the required control current according to the expected damping force. In this study, the Elman neural network was used to establish the inverse model of the MRF damper.

Data Generation
To ensure the training samples can cover the working condition of the MRF damper, the following principals were used: (1) the current range of training samples should include the application range from 0 to 1 A; (2) the samples of displacement and speed range should cover a wide range of applications; and (3) the number of samples should be large enough.
Based on the above principles, because this MRF damper was mainly used on car suspensions, to make the training data relevant to practical applications, a quarter vehicle suspension model was established. The suspension displacement calculated from the road spectrum and its differential result suspension velocity were used as the training data of the neural network. The current data were set as a random white noise signal with a magnitude from 0 to 1 A. The output damping force was obtained using the forward model of the MRF damper established above.

Elman Neural Network
To realize the real-time control of the MRF damper, the structure of the neural network should be as simple as possible. The three inputs were the piston displacement x, the piston velocity v, and the output damping force F. The output variable was the control current I.
The Elman neural network is a multilayer feedforward neural network and is based on the BP neural network and negative feedback system. BP neural networks generally have an input layer, several hidden layers, and an output layer. In addition, the Elman neural network adds an undertaking layer. Compared with BP neural networks, the Elman neural network has better stability and is better at classifying the original data. The model structure is shown in Figure 9 [24]. The basic mathematical model of an Elman neural network is shown as follow [25]: Equations (2) and (3) are the calculation expressions of the output values of the hidden layer and input layer of an Elman network, respectively. w (1) the is weight matrix connecting the hidden layer and the input layer, and w (2) is the weight matrix connecting the undertaking layer and the hidden layer. b (h) is the hidden neuron bias vector, and t is the sampling time. x t is the output of hidden layer. x (c) t is feedback data. The excitation function of the Elman neural network is a nonlinear sigmoid function, which is shown in the following formula: The output layer neurons of the network are linear transfer functions. The final network output is as follows: In this equation, w (3) is the weight matrix connecting the output layer and the hidden layer, b (0) is the output neuron bias vector, and y t is the final network output.
The training process of the Elman neural network is shown in Figure 10.
Equations (2) and (3) are the calculation expressions of the output values of the hidden layer and input layer of an Elman network, respectively. w (1) the is weight matrix connecting the hidden layer and the input layer, and w (2) is the weight matrix connecting the undertaking layer and the hidden layer. b (h) is the hidden neuron bias vector, and t is the sampling time. xt is the output of hidden layer. x (c) t is feedback data. The excitation function of the Elman neural network is a nonlinear sigmoid function, which is shown in the following formula: The output layer neurons of the network are linear transfer functions. The final network output is as follows: In this equation, w (3) is the weight matrix connecting the output layer and the hidden layer, b (0) is the output neuron bias vector, and yt is the final network output.
The training process of the Elman neural network is shown in Figure 10.  Figure 10. Elman neural network training process.

The Inverse Model of the MRF Damper
In the inverse model of the MRF damper, the inputs were piston displacement, piston velocity, and damping force, and the output was current. Therefore, the neuron number of input layer m was 3, and the neuron number of the output layer n was 1. The Elman neural network for the neurons number of the hidden layer is generally calculated using the formula m n a = + + a∈ [1,10], where b is a positive integer. Through multiple trials, the number of neurons in the hidden layer b was 10.
In the network training, the parameters of the Elman neural network, such as learning efficiency, training error, and the number of neurons in each layer, should be initialized first. Then, samples are input to calculate the output and training error of each layer and judge the training error. If the requirements are met or the maximum training times are reached, the training ends. If the requirements are not met, the network parameters

The Inverse Model of the MRF Damper
In the inverse model of the MRF damper, the inputs were piston displacement, piston velocity, and damping force, and the output was current. Therefore, the neuron number of input layer m was 3, and the neuron number of the output layer n was 1. The Elman neural network for the neurons number of the hidden layer is generally calculated using the formula b = (m + n) + a a∈ [1,10], where b is a positive integer. Through multiple trials, the number of neurons in the hidden layer b was 10.
In the network training, the parameters of the Elman neural network, such as learning efficiency, training error, and the number of neurons in each layer, should be initialized first. Then, samples are input to calculate the output and training error of each layer and judge the training error. If the requirements are met or the maximum training times are reached, the training ends. If the requirements are not met, the network parameters are updated and the training continues until the requirements are met or reach the maximum training times. In this training, the target of convergence error was set to 0.001, and the training times was set to 200.
To verify the effectiveness of the Elman neural network, the BP neural network was used to train and predict the MRF damper inverse model, and the prediction results were compared with those of the Elman neural network. Figure 11 shows the original current and the predicted control current using the Elman and BP neural networks, and Figure 12 shows the error comparison of the control current using the Elman and BP neural networks. From the prediction results and error comparison given in Figures 11 and 12, the fitting effect of the control current predicted by the Elman neural network to the original current sample was higher than that of the BP network, and the root mean square (RMS) value of the Elman neural network prediction error was 0.0236, which was less than that of the BP neural network. The prediction error of BP network was 0.0301. In addition, the fitting error using Elman neural network was reduced by 21.59%, which indicates that the Elman neural network inverse model has better performance.
To verify the effectiveness of the Elman neural network, the BP neural network was used to train and predict the MRF damper inverse model, and the prediction results were compared with those of the Elman neural network. Figure 11 shows the original current and the predicted control current using the Elman and BP neural networks, and Figure 12 shows the error comparison of the control current using the Elman and BP neural networks. From the prediction results and error comparison given in Figures 11 and 12, the fitting effect of the control current predicted by the Elman neural network to the original current sample was higher than that of the BP network, and the root mean square (RMS) value of the Elman neural network prediction error was 0.0236, which was less than that of the BP neural network. The prediction error of BP network was 0.0301. In addition, the fitting error using Elman neural network was reduced by 21.59%, which indicates that the Elman neural network inverse model has better performance.

The Dynamic Model of the Vehicle Suspension
The reverse model of the MRF damper was applied to the semi-active control system of a vehicle suspension model. The system can realize on-line control of vehicle suspension vibration. In the process of vehicle driving, the occupant of the vehicle is excited by  To verify the effectiveness of the Elman neural network, the BP neural network was used to train and predict the MRF damper inverse model, and the prediction results were compared with those of the Elman neural network. Figure 11 shows the original current and the predicted control current using the Elman and BP neural networks, and Figure 12 shows the error comparison of the control current using the Elman and BP neural networks. From the prediction results and error comparison given in Figures 11 and 12, the fitting effect of the control current predicted by the Elman neural network to the original current sample was higher than that of the BP network, and the root mean square (RMS) value of the Elman neural network prediction error was 0.0236, which was less than that of the BP neural network. The prediction error of BP network was 0.0301. In addition, the fitting error using Elman neural network was reduced by 21.59%, which indicates that the Elman neural network inverse model has better performance.

The Dynamic Model of the Vehicle Suspension
The reverse model of the MRF damper was applied to the semi-active control system of a vehicle suspension model. The system can realize on-line control of vehicle suspension vibration. In the process of vehicle driving, the occupant of the vehicle is excited by

The Dynamic Model of the Vehicle Suspension
The reverse model of the MRF damper was applied to the semi-active control system of a vehicle suspension model. The system can realize on-line control of vehicle suspension vibration. In the process of vehicle driving, the occupant of the vehicle is excited by vibration in all directions but, generally, the amplitude of vertical vibration is the largest. Therefore, the vertical vibration of the vehicle has great significance for the vehicle ride comfort. Thus, the two degrees of freedom vehicle suspension model, quarter vehicle suspension model, can be used to test the effect of the vehicle semi-active suspension control method for vertical vibration control. It can also be used to verify the control effect of the MRF damper inverse model to achieve continuous adjustable damping. Therefore, this study established a quarter vehicle semi-active suspension model for analysis. The differential equations of vehicle suspension model are shown as follows: where m b is the sprung mass, m w is the unsprung mass, k s is the suspension spring stiffness, k t is the tire stiffness, F d is the adjustable damping force, in the semi-active control, this force is the sum of the basis damping force (input current is zero) and the semi-active control force, x w is the displacement of the unsprung mass, x b is the displacement of the sprung mass, and x r is the road excitation. The parameters in the system were set as follows: m b = 400 kg, m w = 50 kg, k s = 35,000 N/m, and k t = 200,000 N/m.

Semi-Active Control Algorithm
As a classical control strategy of the semi-active suspension, the skyhook control method has many advantages such as simple implementation, fewer control variables and easy measurement. The skyhook control strategy is: where c sky is the skyhook damping controlled by the current, V rel is the relative velocity of vehicle suspension (sprung mass velocity minus unsprung mass velocity), and V s is the sprung mass velocity of vehicle suspension. For traditional skyhook control, the damping force is generated by a fixed current. That is, when the motion of suspension meets the requirements of the output damping force shown in Equation (7), turn on the current switch, input a fixed current, and output the damping force. Otherwise, turn off the current switch and only output the minimum damping force. However, the skyhook control also has disadvantages: (1) it produces impact load on the suspension system during the damping force switching; (2) the MRF damper can only change between the fixed current and no current under the skyhook control and is thus unable to fully utilize the ability of the MRF damper to adjust damping force continuously and consumes more electric energy. Therefore, a fuzzy algorithm-based skyhook control was examined to maximize the advantage of the MRF damper's continuously adjustable damping [26].
The input variables of the fuzzy controller were the sprung mass velocity vs. and the relative velocity V rel , and the output was the damping force factor. This factor was multiplied by the skyhook damping coefficient to obtain the fuzzy control damping coefficient, and then the coefficient was multiplied by V rel at this moment to obtain the output control force. Based on multiple calculations on the suspension model, the skyhook damping coefficient was selected as 2000 Ns/m.
The input and output of the fuzzy controller must be represented by fuzzy language variables, so the numerical variables should be clear. According to the characteristics of the vehicle suspension system and MRF damper damping force output law, the inputs of V rel and V s . were divided into five levels: negative large, negative small, zero, positive small, positive large. These five levels are defined as NB, NS, ZE, PS, and PB, respectively. In addition, the output of fuzzy control was divided into four levels: zero, small, medium, and large. These four levels are defined as ZE, S, M and B, respectively.
The fuzzy logic control rules are based on the skyhook strategy, as shown in Table 3, and the structure of the fuzzy controller is shown in Figure 13.  Figure 13. Structure of the fuzzy controller As shown in Figure 13, the inputs of the fuzzy control were Vrel and Vs, and the discourse domains of the input variables were set as [−1, 1]. The output was set as [0, 1], and then the output was multiplied by the scale factor, Ku, to obtain the fuzzy control output. The control rules and reasoning mechanism of a designed fuzzy controller cannot be adjusted, but the quantization factors Ka and Kv and the scale factor, Ku, can be adjusted, and have significant influence on the controller performance.
At present, many methods exist to select the quantification factors Ka and Kv and the scale factor Ku; however, there is no unified standard. Selection mainly depends on experience. For the semi-active suspension system, it needs to be adjusted repeatedly to achieve a better control effect, which has certain limitations. Therefore, it is necessary to find a feasible and effective method to select quantification factors and scale factors. Because the most important performance index of vehicle suspension is to control the acceleration of the vehicle body, the optimization design method can be used to search for the best quantification factors and scale factors according to the optimization objective and thus obtain the best control effect.
This study used the GWO method to optimize the quantification factors Ka, Kv, and the scale factor Ku, and obtain the best control effect.
The GWO algorithm is a population intelligent optimization method first developed by Mirjalili [27]. The GWO algorithm has the characteristics of simple principle, less calculation requirements. It has been applied in face recognition [27], traffic flow prediction, power system control, and other fields.
The grey wolf population can be divided into four levels: α wolf, β wolf, δ wolf, and ω wolf. The GWO algorithm defines the optimal solution, suboptimal solution, and the third optimal solution as α wolf, β wolf, and δ wolf, which guide hunting; other candidate solutions are defined as ω wolf, which follows α wolf, β wolf, and δ wolf [28]. The predation behavior of grey wolf is divided into encircling, hunting, and attacking. The behavior of wolves surrounding their prey is expressed as [29]: Figure 13. Structure of the fuzzy controller.
As shown in Figure 13, the inputs of the fuzzy control were V rel and V s , and the discourse domains of the input variables were set as [−1, 1]. The output was set as [0, 1], and then the output was multiplied by the scale factor, K u , to obtain the fuzzy control output. The control rules and reasoning mechanism of a designed fuzzy controller cannot be adjusted, but the quantization factors K a and K v and the scale factor, K u , can be adjusted, and have significant influence on the controller performance.
At present, many methods exist to select the quantification factors K a and K v and the scale factor K u ; however, there is no unified standard. Selection mainly depends on experience. For the semi-active suspension system, it needs to be adjusted repeatedly to achieve a better control effect, which has certain limitations. Therefore, it is necessary to find a feasible and effective method to select quantification factors and scale factors. Because the most important performance index of vehicle suspension is to control the acceleration of the vehicle body, the optimization design method can be used to search for the best quantification factors and scale factors according to the optimization objective and thus obtain the best control effect.
This study used the GWO method to optimize the quantification factors K a , K v , and the scale factor K u , and obtain the best control effect.
The GWO algorithm is a population intelligent optimization method first developed by Mirjalili [27]. The GWO algorithm has the characteristics of simple principle, less calculation requirements. It has been applied in face recognition [27], traffic flow prediction, power system control, and other fields.
The grey wolf population can be divided into four levels: α wolf, β wolf, δ wolf, and ω wolf. The GWO algorithm defines the optimal solution, suboptimal solution, and the third optimal solution as α wolf, β wolf, and δ wolf, which guide hunting; other candidate solutions are defined as ω wolf, which follows α wolf, β wolf, and δ wolf [28]. The predation behavior of grey wolf is divided into encircling, hunting, and attacking. The behavior of wolves surrounding their prey is expressed as [29]: where D is the distance between the wolf and its prey; X(t) and X p (t) are the current positions of the wolf and prey respectively; t is the current number of iterations; A and C are the coefficient vectors. The specific definitions of A and C are as follows: X(t + 1) = (X 1 + X 2 + X 3 )/3 (11) where D α , D β , and D δ are the distances of the α wolf, β wolf, and δ wolf from the other grey wolves, respectively; X α , X β , and X δ are the positions of the α wolf, β wolf, and δ wolf, respectively; A 1 , A 2 , and A 3 are coefficient vectors; C 1 , C 2 , and C 3 are random vectors; rand1 and rand2 are random vectors in the interval [0, 1], which are generated randomly in the optimization process; X is the position of the current grey wolf individual. Finally, the wolves attack the prey. The value of the convergence factor a in Formula (8) decreases linearly, and decreases from 2 to 0 with the value of t. When |A| ≤ 1, the wolf pack attacks the prey, highlighting the local search; when A > 1, the wolf pack gives up hunting the current prey and spreads out, highlighting the global search.
The performance of the suspension system has a direct impact on vehicle performance, including ride comfort, handling stability, and ride comfort. Generally, the performance of the suspension system is evaluated based on vertical acceleration, suspension dynamic deflection, roll angle and roll acceleration, pitch angle and pitch acceleration, and tire dynamic load. The vertical displacement and acceleration are usually larger than those in the roll and pitch directions. This study mainly examined the vertical motion caused by the impact when the vehicle is running, and evaluated the suspension performance in terms of the following three aspects: body acceleration, suspension dynamic deflection, and tire dynamic load.
Among the three parameters, the vehicle body acceleration directly represents the vibration intensity transmitted to the vehicle body, which is the most important index to evaluate-the ride comfort of the vehicle. The suspension dynamic deflection affects the probability of hitting the suspension stop; if this parameter is too large, the driving stability of the vehicle will decrease. The tire dynamic load characterizes the ability of the vehicle to resist the road interference; if this parameter is too high, the handling stability of the vehicle will be reduced. In this study, the root mean square (RMS) value of vertical acceleration was chosen as the optimization target. The fitness function is shown as follows: where a b is the vertical acceleration of vehicle body. After establishing the fitness function, the GWO algorithm was used to optimize the quantification factors K a , K v , and the scale factor K u of fuzzy control. The optimization process was as follows: (1) Initialization of parameters, which set the optimized parameters of quantification factors K a , K v , and the scale factor K u , had a range of values of [1,10], the maximum number of iterations was 200, and the position of the initialized wolves was X = (X 1 , X 2 , . . . , X N ).
(2) Calculate the adaptive function values of each wolf f i (i = 1, 2, . . . , N), filter the top three grey wolf individuals with adaptive values (α wolf, β wolf, and δ wolf), and determine their locations X α , X β , and X δ .
(3) The distance between the remaining grey wolf and X α , X β , and X δ was calculated by Formula (9), and the location of grey wolf α, β, and δ was updated in a timely manner based on Formulas (10) and (11).
(4) Update parameters α, A, and C according to Equation (9). (5) If the algorithm exceeds the maximum number of iterations, stop the loop and output the optimal solution X α ; otherwise, return to Step (2) to continue the loop.
After the control force is obtained by the optimized fuzzy skyhook method, the semiactive control force is summed with the basic damping force of the MRF damper to obtain the ideal output damping force. Then the ideal output damping force, the relative velocity V rel , and relative displacement d rel are simultaneously updated to the inverse model of the MRF damper to obtain its control current to realize the semi-active control. The control flow is shown in Figure 14.
output the optimal solution X α ; otherwise, return to Step (2) to continue the loop.
After the control force is obtained by the optimized fuzzy skyhook method, the semiactive control force is summed with the basic damping force of the MRF damper to obtain the ideal output damping force. Then the ideal output damping force, the relative velocity Vrel, and relative displacement drel are simultaneously updated to the inverse model of the MRF damper to obtain its control current to realize the semi-active control. The control flow is shown in Figure 14.  Figure 14. Control process of the semi-active suspension.

Semi-Active Control Simulation
The C class road excitation according to ISO 8608: 2016 [30] was used in this section as shown in Figure 15. The passive vibration reduction method, the traditional skyhook control method, and the optimized fuzzy skyhook method were simulated and compared after optimization using the GWO method. The optimized fuzzy control parameters were Ka = 9.82, Kv = 3.11, and Ku = 1.10. The passive damping method was used such that the current output of the MRF damper was always zero, and the suspension damping could not be changed in the passive vibration reduction state. The traditional skyhook control method used the approach described in Equation (7) to control the vibration. The time histories and RMS values of the vertical acceleration, suspension dynamic deflection, the tire dynamic load under the passive state, the traditional skyhook control, and the GWO fuzzy skyhook control are shown in Table 4 and Figures 16-18.

Semi-Active Control Simulation
The C class road excitation according to ISO 8608: 2016 [30] was used in this section as shown in Figure 15. The passive vibration reduction method, the traditional skyhook control method, and the optimized fuzzy skyhook method were simulated and compared after optimization using the GWO method. The optimized fuzzy control parameters were K a = 9.82, K v = 3.11, and K u = 1.10. The passive damping method was used such that the current output of the MRF damper was always zero, and the suspension damping could not be changed in the passive vibration reduction state. The traditional skyhook control method used the approach described in Equation (7) to control the vibration.

Semi-Active Control Simulation
The C class road excitation according to ISO 8608: 2016 [30] was as shown in Figure 15. The passive vibration reduction method, the t control method, and the optimized fuzzy skyhook method were simul after optimization using the GWO method. The optimized fuzzy contr Ka = 9.82, Kv = 3.11, and Ku = 1.10. The passive damping method was current output of the MRF damper was always zero, and the suspens not be changed in the passive vibration reduction state. The tradition method used the approach described in Equation (7) to control the vib  The time histories and RMS values of the vertical acceleration, suspension dynamic deflection, the tire dynamic load under the passive state, the traditional skyhook control, and the GWO fuzzy skyhook control are shown in Table 4 and Figures 16-18.      According to the time domain diagrams and RMS values of vertical acceleration, suspension dynamic deflection and tire dynamic load under the three control methods (as shown in Table 3), using traditional skyhook control, the RMS value of vertical accelera- According to the time domain diagrams and RMS values of vertical acceleration, suspension dynamic deflection and tire dynamic load under the three control methods (as shown in Table 3), using traditional skyhook control, the RMS value of vertical acceleration was reduced from 1.3816 to 1.3241 m/s 2 , the dynamic deflection acceleration was reduced from 0.0111 to 0.0091 m, and tire dynamic load was reduced from 751.1950 to 713.3522 N. This shows that skyhook control can effectively improve suspension performance. After using the fuzzy skyhook control strategy and optimizing the parameters by GWO algorithm, the RMS value of vertical acceleration was further reduced to 1.1864 m/s 2 , the dynamic deflection was further reduced to 0.0090 m, and the tire dynamic load was further reduced to 705.3631 N. Compared with passive control, these values were reduced by 14.13%, 18.92%, and 6.10%, respectively. This result proves that the optimal control method proposed in this paper was effective.
In addition to the normal road surface, driving on a potholed road often results in low-frequency large amplitude shock vibration, which has an adverse impact on vehicle ride comfort. Thus, we also carried out a simulation analysis of the suspension vibration reduction of this condition. This excitation signal is shown in Figure 19. According to the time domain diagrams and RMS values of vertical acc pension dynamic deflection and tire dynamic load under the three contro shown in Table 3), using traditional skyhook control, the RMS value of ver tion was reduced from 1.3816 to 1.3241 m/s 2 , the dynamic deflection accele duced from 0.0111 to 0.0091 m, and tire dynamic load was reduced from 713.3522 N. This shows that skyhook control can effectively improve suspe mance. After using the fuzzy skyhook control strategy and optimizing the p GWO algorithm, the RMS value of vertical acceleration was further redu m/s 2 , the dynamic deflection was further reduced to 0.0090 m, and the tire was further reduced to 705.3631 N. Compared with passive control, these reduced by 14.13%, 18.92%, and 6.10%, respectively. This result proves tha control method proposed in this paper was effective.
In addition to the normal road surface, driving on a potholed road o low-frequency large amplitude shock vibration, which has an adverse imp ride comfort. Thus, we also carried out a simulation analysis of the suspen reduction of this condition. This excitation signal is shown in Figure 19. The time histories and RMS values of the body vertical acceleration, su namic deflection, and tire dynamic using three control methods under poth citation are shown in Table 5   The time histories and RMS values of the body vertical acceleration, suspension dynamic deflection, and tire dynamic using three control methods under potholed road excitation are shown in Table 5 and Figures 20-22.       From Table 5 and Figures 20-22, it can be seen that under the potholed road condition, using the optimized fuzzy skyhook control strategy, the RMS value of vertical acceleration was reduced from 1.3388 to 1.0356 m/s 2 , dynamic deflection acceleration was reduced  From Table 5 and Figures 20-22, it can be seen that under the potholed road condition, using the optimized fuzzy skyhook control strategy, the RMS value of vertical acceleration was reduced from 1.3388 to 1.0356 m/s 2 , dynamic deflection acceleration was reduced from 0.0110 to 0.0093 m, and tire dynamic load was reduced from 737.9775 to 680.7965 N. Compared with passive control, these values were reduced by 22.65%, 15.46%, and 7.75%, respectively. It can be seen from Figures 20 and 21, that after receiving large amplitude vibration excitation, the optimized fuzzy skyhook control strategy can more quickly lower the acceleration of the vehicle body and decay the dynamic deflection of the suspension. These results prove that the optimal control method proposed in this paper was effective. This shows that the proposed method can also obtain a better damping effect under lowfrequency and large amplitude excitation and can improve the ride comfort of vehicles.

Conclusions
In this study, the vehicle magnetorheological semi-active suspension with continuously adjustable MRF damping was examined. The inverse model of the MRF damper and the fuzzy skyhook control optimized by the GWO algorithm were investigated through experiment with numerical simulation, aiming to improve the ride comfort and handling stability of the vehicle. The main works are shown as follows: 1. The inverse model of the MRF damper developed using the Elman neural network was established to calculate the control current from the expected output damping force. The inverse model can effectively convert the expected damping force into the control current, and apply the control current to the actual MRF damper to obtain the actual damping force. Thus, the damping force can be continuously adjusted to reduce the suspension vibration. Compared with the BP neural network model, the Elman neural network model reduced the local detail error of the inverse model and improved the prediction accuracy of the control current; 2. A quarter vehicle suspension model was established in which the MRF damper was used for the suspension model, and semi-active control analysis was carried out. The results show that the semi-active vibration control of the MRF damper can effectively improve the ride comfort of the suspension. The vertical vehicle acceleration, suspension dynamic reflection, and tire dynamic load can be reduced simultaneously using the skyhook semi-active control method; 3. A fuzzy skyhook control method based on GWO was proposed, aiming to overcome the shortcomings that the quantification factor and scale factor in the fuzzy control are determined through expert experience. The fuzzy skyhook control method optimized by the GWO algorithm was proposed to optimize the quantification factor and scale factor. The normal road and pothole road conditions were selected for simulation analysis. The analysis results show that this optimized skyhook control method can reduce the vertical acceleration, suspension dynamic reflection, and tire dynamic load. For vertical acceleration, the most important component of the ride comfort index, the control effect of the optimized fuzzy skyhook control method proposed in this paper was clearly better than that of the traditional skyhook control method.
Future research will focus on further improving the semi-active control algorithm to obtain a better damping effect. In addition, the research will be extended to the full-vehicle semi-active suspension system, so that the MRF damper can control the vertical, roll, and pitch vibration of the vehicle simultaneously.