Fuzzy Sliding Mode Control of Vehicle Magnetorheological Semi-Active Air Suspension

: In order to reduce vehicle vibration during driving conditions, a fuzzy sliding mode control strategy (FSMC) for semi-active air suspension based on the magnetorheological (MR) damper is proposed. The MR damper used in the semi-active air suspension system was tested and analyzed. Based on the experimental data, the genetic algorithm was used to identify the parameters of the improved hyperbolic tangent model, which was derived for the MR damper. At the same time, an adaptive neuro fuzzy inference system (ANFIS) was used to build the reverse model of the MR damper. The model of a quarter vehicle semi-active air suspension system equipped with a MR damper was established. Aiming at the uncertainty of the air suspension system, fuzzy control was used to adjust the boundary layer of the sliding mode control, which can effectively suppress the inﬂuence of chattering on the control accuracy and ensure system stability. Taking random road excitation and impact road excitation as the input signal, the simulation analysis of passive air suspension, semi-active air suspension based on SMC and FSMC was carried out, respectively. The results show that the semi-active air suspension based on FSMC has better vibration attenuating performance and ride comfort.


Introduction
The suspension system is the transmission device connecting the wheels and the car body, which mainly undertakes the body mass and ensures good contact between the tire and the ground, as well as reducing the uneven excitation from the road surface [1]. Compared with the traditional suspension system, the air spring suspension system has the advantages of being light weight, being low noise and having favorable vibration isolation performance. The air spring has the characteristics of non-linearity and self-adaptation, which makes the deviation frequency of the vehicle body mass in a stable state under the condition of load changes and obtains favorable ride comfort and handling stability [2][3][4][5].
As a new type of intelligent semi-active control device, the MR damper has become a research hotspot and can be well applied to semi-active air suspension, due to its good performance, i.e., simple structure, controllable damping force, low power consumption as well as rapid response time [6][7][8][9]. In industrial application, the accurate establishment of the mathematical model of MR damper is an important step for the design and analysis of the semi-active controller. However, the influences of movement condition of the piston, applied magnetic field, displacement amplitude and frequency on the dynamic performance of MR damper are inevitable, and the MR damper has the strong nonlinear characteristics, such as yield, hysteresis and saturation. It is difficult to establish the accurate dynamic model of the MR damper. So far, different mechanical models of MR damper have been proposed, which are mainly divided into two categories: non-parametric models and parametric models [10]. Parametric modeling of the MR damper is a process Appropriate control algorithms can make the semi-active air suspension almost achieve the same vibration attenuation effect as the active air suspension. In recent years, the sliding mode control (SMC) with strong anti-interference ability and good robustness attracted more and more attention in related fields. Choi et al. designed a sliding mode controller for the vehicle seat and considered the human body mass as the uncertainty factor of the system [15]. Through simulation analysis, it is concluded that the sliding mode controller has favorable robustness to the uncertainty of the system, and the ride comfort is also improved. Pang et al. proposed an improved adaptive sliding mode-based fault-tolerant control design, the results showing that the dynamics performance of half vehicle active suspensions with parametric uncertainties and actuator faults in the context of external road disturbances are improved [16]. Wang et al. presented a practical terminal sliding mode control framework based on an adaptive disturbance observer for the active suspension systems, which requires no exact feedback linearization about the suspension dynamics [17]. For the tracking control problem of vehicle suspension, Fei et al. derived a robust design method of adaptive sliding mode control [18]. Under this approach, the influence of parameter uncertainties and external disturbances on the system performance can be reduced, and the system robustness is also improved. Du et al. proposed a terminal sliding mode control approach for an active suspension system, which has an ability to reach the sliding surface in a finite time to achieve high control accuracy [19]. Combining the principle of fractional order control with sliding mode control, Nguyen et al. proposed a corresponding fractional order sliding mode control strategy for vehicle suspension and verified its effectiveness through simulation and experiment [20]. Zheng et al. combined the advantages of fuzzy control and sliding mode control and confirmed the effectiveness of the fuzzy sliding mode controller [21]. In addition, many other control strategies are applied to the control of suspension systems, which include switching control, adaptive robust control, linear-quadratic-Gaussian (LQG) control, and intelligent control methods [22][23][24][25][26][27].
Sliding mode variable structure control has good robustness, strong anti-interference ability, and can be well adapted to air suspension systems. However, when the actual controlled system switches the system structure at a higher frequency, the chattering phenomenon often occurs, which affects the accuracy of the sliding mode variable structure control and the system stability [28]. Fuzzy control can compensate for the uncertainty of the system and adjust the boundary layer of the sliding mode control by making fuzzy control rules to realize the imprecise control. In order to suppress the influence of chattering on the control precision and ensure the system stability, a FSMC controller was developed on the combination of fuzzy control and sliding mode variable structure control.
In this paper, the vibration test system was built to test and analyze the performance of the MR damper used in the air suspension system. Using the obtained test data, the genetic algorithm was applied to identify the parameters of the improved hyperbolic tangent model, and the inverse model of the MR damper was built, using an adaptive neurofuzzy inference system simultaneously. Subsequently, the air spring and MR damper were combined in vehicle suspension. The air spring has good nonlinear elastic characteristics and has different elastic characteristics under a heavy load and light load, which can improve ride comfort and handling stability. Then, a quarter car MR semi-active air suspension system model was established. In order to reduce the chattering problem in the sliding mode control and ensure system stability, a fuzzy sliding mode control strategy of semi-active air suspension was proposed. The fuzzy control can compensate for the uncertainty of the system and reduce the chattering problem in sliding mode control and ensure the stability and safety of the vehicle suspension. Finally, the numerical simulation was carried out to verify the effectiveness of the proposed control strategy under the conditions of excitations of the random road profile and impact road profile.

Damping Performance Test of MR Damper
In order to test the damping characteristics of MR damper, a vibration test system was established as shown in Figure 1, which is mainly composed of a damping suspension test bench, DC power supply, electro-hydraulic servo controller, computer and MR damper. The DC power supply was used to supply the required DC current to the coil in the MR damper. The hydraulic station provided hydraulic power for the vibration test bench. The computer was used to collect experimental data, and the vibration test bench was controlled by the electro-hydraulic servo controller. gent model, and the inverse model of the MR damper was built, using an adaptive neu fuzzy inference system simultaneously. Subsequently, the air spring and MR damp were combined in vehicle suspension. The air spring has good nonlinear elastic charact istics and has different elastic characteristics under a heavy load and light load, which c improve ride comfort and handling stability. Then, a quarter car MR semi-active air s pension system model was established. In order to reduce the chattering problem in t sliding mode control and ensure system stability, a fuzzy sliding mode control strategy semi-active air suspension was proposed. The fuzzy control can compensate for the u certainty of the system and reduce the chattering problem in sliding mode control a ensure the stability and safety of the vehicle suspension. Finally, the numerical simulati was carried out to verify the effectiveness of the proposed control strategy under the co ditions of excitations of the random road profile and impact road profile.

Damping Performance Test of MR Damper
In order to test the damping characteristics of MR damper, a vibration test syste was established as shown in Figure 1, which is mainly composed of a damping suspensi test bench, DC power supply, electro-hydraulic servo controller, computer and M damper. The DC power supply was used to supply the required DC current to the coil the MR damper. The hydraulic station provided hydraulic power for the vibration t bench. The computer was used to collect experimental data, and the vibration test ben was controlled by the electro-hydraulic servo controller.  In this experiment, a sinusoidal excitation signal with a frequency of 1 Hz and a plitudes of 5 mm, 7.5 mm, 10 mm was adopted. The currents were 0 A, 0.25 A, 0.5 A, 0 A and 1 A, provided by a DC regulated power supply. As shown in Figure 2, the variati of the output damping force with displacement and velocity of the MR damper und different excitation currents is obtained. Figure 2a shows that the output damping fo of the MR damper increases with the increase in the excitation current. This is main because the magnetic flux density in the damping channel increases with the increase the excitation current. The output damping force varies little with the displacement, a it is distributed symmetrically with respect to the relative equilibrium position. As show in Figure 2b, the damping force of the MR damper increases with the increase in the ex tation current and vibration speed, but the influence of the speed is much smaller th that with the current. The proposed MR damper has a wide range of damping force a justments, and thus, it can be used in a vehicle suspension system. In this experiment, a sinusoidal excitation signal with a frequency of 1 Hz and amplitudes of 5 mm, 7.5 mm, 10 mm was adopted. The currents were 0 A, 0.25 A, 0.5 A, 0.75 A and 1 A, provided by a DC regulated power supply. As shown in Figure 2, the variation of the output damping force with displacement and velocity of the MR damper under different excitation currents is obtained. Figure 2a shows that the output damping force of the MR damper increases with the increase in the excitation current. This is mainly because the magnetic flux density in the damping channel increases with the increase in the excitation current. The output damping force varies little with the displacement, and it is distributed symmetrically with respect to the relative equilibrium position. As shown in Figure 2b, the damping force of the MR damper increases with the increase in the excitation current and vibration speed, but the influence of the speed is much smaller than that with the current. The proposed MR damper has a wide range of damping force adjustments, and thus, it can be used in a vehicle suspension system.

Improved Hyperbolic Tangent Model of MR Damper
The improved hyperbolic tangent model can comprehensively describe the damping characteristics of the MR damper. Compared with the hyperbolic tangent model, the physical quantities of this model have a clearer physical significance, and there are fewer unknown parameters, which makes it easy for subsequent parameter identification. The improved hyperbolic tangent model is shown in Figure 3, and its mathematical expression is as follows: where a1 and k are both proportional factors, a1 is related to the current and k affects the width of the hysteresis loop; a2 and a3 are damping coefficients that affect the curve variation trend of the pre-yield zone and post-yield zone, respectively; and f0 denotes the offset damping force.
. Improved hyperbolic tangent model. Figure 4 shows the parameter identification process using the genetic algorithm. Taking the mean square deviation of the difference between the actual damping force and the simulated damping force as the objective function, the smallest mean square deviation is the optimal solution, and the optimal model of the MR damper is obtained through genetic operation selection, crossover and mutation.

Improved Hyperbolic Tangent Model of MR Damper
The improved hyperbolic tangent model can comprehensively describe the damping characteristics of the MR damper. Compared with the hyperbolic tangent model, the physical quantities of this model have a clearer physical significance, and there are fewer unknown parameters, which makes it easy for subsequent parameter identification. The improved hyperbolic tangent model is shown in Figure 3, and its mathematical expression is as follows: where a 1 and k are both proportional factors, a 1 is related to the current and k affects the width of the hysteresis loop; a 2 and a 3 are damping coefficients that affect the curve variation trend of the pre-yield zone and post-yield zone, respectively; and f 0 denotes the offset damping force.

Improved Hyperbolic Tangent Model of MR Damper
The improved hyperbolic tangent model can comprehensive characteristics of the MR damper. Compared with the hyperb physical quantities of this model have a clearer physical significa unknown parameters, which makes it easy for subsequent param improved hyperbolic tangent model is shown in Figure 3, and its m is as follows: where a1 and k are both proportional factors, a1 is related to the c width of the hysteresis loop; a2 and a3 are damping coefficients th tion trend of the pre-yield zone and post-yield zone, respectively; damping force.  Figure 4 shows the parameter identification process using the ing the mean square deviation of the difference between the actua simulated damping force as the objective function, the smallest m the optimal solution, and the optimal model of the MR damper netic operation selection, crossover and mutation.  Figure 4 shows the parameter identification process using the genetic algorithm. Taking the mean square deviation of the difference between the actual damping force and the simulated damping force as the objective function, the smallest mean square deviation is the optimal solution, and the optimal model of the MR damper is obtained through genetic operation selection, crossover and mutation.  The fitness function of the improved hyperbolic tangent model is denoted as follows:

Improved Parameter Identification of Hyperbolic Tangent Model
where Fsimi and Fexpi are the simulation damping forces and experimental damping forces of the ith point, respectively, and m is the number of experimental points. The identification parameters were processed by using the real coding method. The scattered crossover method was used to determine the crossover probability by repeated experiments, and the crossover probability was set to 0.8. To improve the local search ability of the genetic algorithm, the Gaussian function was used for mutation with a scale of 0.5 as well as a shrink of 0.7, and the number of iterations was 100.
The improved hyperbolic tangent model has five parameters to be identified, which is expressed as follows: The test results of the mechanical characteristics of the MR damper under sinusoidal excitation with an amplitude of 10 mm and frequency of 1 Hz were selected as the identification data. Parameter identification was carried out for the data of 5 groups under different control currents so as to obtain parameter identification results under different current, as shown in Table 1. By comparing and analyzing the identification results of different current parameters, the variation rules of each parameter were found, and the fitting values of each parameter in the improved hyperbolic tangent model were also obtained. Observing Table 1, parameters a2, k, and f0 are not changed significantly under the different current. Therefore, a2, k, and f0 can be considered constant values, and their mean values a2 = 3874.38, k = 0.31, and f0 = 36.12 are selected. Parameters a1 and a3 increase linearly with the increase in current, and the relationship between the parameters and current is obtained through fitting, as shown in Figure 5. The fitness function of the improved hyperbolic tangent model is denoted as follows: where F simi and F expi are the simulation damping forces and experimental damping forces of the ith point, respectively, and m is the number of experimental points. The identification parameters were processed by using the real coding method. The scattered crossover method was used to determine the crossover probability by repeated experiments, and the crossover probability was set to 0.8. To improve the local search ability of the genetic algorithm, the Gaussian function was used for mutation with a scale of 0.5 as well as a shrink of 0.7, and the number of iterations was 100.
The improved hyperbolic tangent model has five parameters to be identified, which is expressed as follows: The test results of the mechanical characteristics of the MR damper under sinusoidal excitation with an amplitude of 10 mm and frequency of 1 Hz were selected as the identification data. Parameter identification was carried out for the data of 5 groups under different control currents so as to obtain parameter identification results under different current, as shown in Table 1. By comparing and analyzing the identification results of different current parameters, the variation rules of each parameter were found, and the fitting values of each parameter in the improved hyperbolic tangent model were also obtained. Observing Table 1, parameters a 2 , k, and f 0 are not changed significantly under the different current. Therefore, a 2 , k, and f 0 can be considered constant values, and their mean values a 2 = 3874.38, k = 0.31, and f 0 = 36.12 are selected. Parameters a 1 and a 3 increase linearly with the increase in current, and the relationship between the parameters and current is obtained through fitting, as shown in Figure 5. As shown in Figure 5, parameters a1 and a3 are linearly correlated with the current, which can be obtained by linear fitting as shown below: The expression of the improved hyperbolic tangent model after parameter identification is obtained as follows: According to the improved hyperbolic tangent model expressed in Equation (5) after parameter identification, the mechanical properties of the MR damper was simulated and compared with the test results, as shown in Figure 6. The force-displacement curve of the improved hyperbolic tangent model after parameter identification has a high degree of fit with the test curve. For the force-velocity curve, the nonlinear hysteretic characteristics of the MR dampers in low, and the high-speed regions are well described. Therefore, the hyperbolic tangent model can be applied to the subsequent semi-active control.  As shown in Figure 5, parameters a 1 and a 3 are linearly correlated with the current, which can be obtained by linear fitting as shown below: The expression of the improved hyperbolic tangent model after parameter identification is obtained as follows: According to the improved hyperbolic tangent model expressed in Equation (5) after parameter identification, the mechanical properties of the MR damper was simulated and compared with the test results, as shown in Figure 6. The force-displacement curve of the improved hyperbolic tangent model after parameter identification has a high degree of fit with the test curve. For the force-velocity curve, the nonlinear hysteretic characteristics of the MR dampers in low, and the high-speed regions are well described. Therefore, the hyperbolic tangent model can be applied to the subsequent semi-active control. As shown in Figure 5, parameters a1 and a3 are linearly correlated with the current, which can be obtained by linear fitting as shown below: The expression of the improved hyperbolic tangent model after parameter identification is obtained as follows: According to the improved hyperbolic tangent model expressed in Equation (5) after parameter identification, the mechanical properties of the MR damper was simulated and compared with the test results, as shown in Figure 6. The force-displacement curve of the improved hyperbolic tangent model after parameter identification has a high degree of fit with the test curve. For the force-velocity curve, the nonlinear hysteretic characteristics of the MR dampers in low, and the high-speed regions are well described. Therefore, the hyperbolic tangent model can be applied to the subsequent semi-active control.

ANFIS Inverse Model of MR Damper
The ANFIS is a multilayer feedforward network that uses neural network learning algorithms and fuzzy reasoning to map an input space to an output space [29]. The ANFIS algorithm has a strong adaptive ability and self-learning ability, and it has a strong approximation ability to the nonlinear system. Therefore, ANFIS is suitable for the modeling of complex systems, such as MR dampers. Figure 7 shows the structure of ANFIS. There are five layers in total, and the function types of the nodes in each layer are the same.

ANFIS Inverse Model of MR Damper
The ANFIS is a multilayer feedforward network that uses neural netwo algorithms and fuzzy reasoning to map an input space to an output space [29]. algorithm has a strong adaptive ability and self-learning ability, and it has a proximation ability to the nonlinear system. Therefore, ANFIS is suitable for th of complex systems, such as MR dampers. Figure 7 shows the structure of AN are five layers in total, and the function types of the nodes in each layer are the Layer 1 is used to fuzzify the input signal, and its node output is defined where o1,i is the output of the ith node of layer 1, and x represents the input o respectively. Ai (or Bi−2) represents fuzzy sets (negative big, negative small, p positive small, etc.). The membership function of the fuzzy set is as follows: c is the premise parameter used to adjust the shape of the fun Layer 2 is represented by π, whose function is to calculate the applicabi rule wi. The specific operation is to take the product of the membership deg input signal as the applicability of this rule, as shown in the equation below.   Layer 1 is used to fuzzify the input signal, and its node output is defined as follows: where o 1,i is the output of the ith node of layer 1, and x represents the input of the node, respectively. A i (or B i−2 ) represents fuzzy sets (negative big, negative small, positive big, positive small, etc.). The membership function of the fuzzy set is as follows: The above equation is bell-shaped functions (the maximum is 1, the minimum is 0), where {a i , b i , c i } is the premise parameter used to adjust the shape of the function. Layer 2 is represented by π, whose function is to calculate the applicability of each rule w i . The specific operation is to take the product of the membership degree of each input signal as the applicability of this rule, as shown in the equation below.
Layer 3 is represented by circle N as shown in the Figure 7. Layer 3 is used to normalize the applicability of each rule, as shown in Equation (9).
The function of Layer 4 is to calculate the output of each rule. The node function is expressed as follows: where parameters {p i , q i , r i } represent an adjustable parameter set. Layer 5 sums up all the inputs of the node and outputs the results to calculate the total output of the system.
ANFIS can obtain the required fuzzy rules and membership function after repeated training based on a large amount of data. Usually, the BP algorithm or hybrid algorithm is used to train the data.
The adaptive neuro fuzzy inference system (ANFIS) was used to establish the inverse model of the MR damper in this paper [30]. The training schematic diagram of the ANFIS inverse model of the MR damper is shown in Figure 8. In order to establish the inverse model of ANFIS, the velocity, displacement and expected damping force at the current moment as well as the velocity and desired damping force provided by the forward model at the previous moment were selected as the input, while the control current at the current moment was taken as the output.
Appl. Sci. 2021, 11, x FOR PEER REVIEW 8 of 22 Layer 5 sums up all the inputs of the node and outputs the results to calculate the total output of the system.
ANFIS can obtain the required fuzzy rules and membership function after repeated training based on a large amount of data. Usually, the BP algorithm or hybrid algorithm is used to train the data.
The adaptive neuro fuzzy inference system (ANFIS) was used to establish the inverse model of the MR damper in this paper [30]. The training schematic diagram of the ANFIS inverse model of the MR damper is shown in Figure 8. In order to establish the inverse model of ANFIS, the velocity, displacement and expected damping force at the current moment as well as the velocity and desired damping force provided by the forward model at the previous moment were selected as the input, while the control current at the current moment was taken as the output.  Two sets of data were needed to model a MR damper, using an adaptive neuro-fuzzy reasoning system. One set was used to train the inverse model of ANFIS, and the other set was used to test the accuracy of the model. A sinusoidal function with an amplitude of 20 mm was used for the displacement signal, and a sinusoidal function with an amplitude of 0~2 A was used for the control current signals. There were 1500 data points of which 1000 points were used for training and 500 points were used for validation. Figures 9 and 10 show the comparison between the experimental and predicted values as well as the error of the ANFIS inverse dynamics model of the MR damper, respectively. The results illustrate that the predictive control current of the ANFIS inverse model can track the target control current well, and verify the effectiveness of the model. Two sets of data were needed to model a MR damper, using an adaptive neuro-fuzzy reasoning system. One set was used to train the inverse model of ANFIS, and the other set was used to test the accuracy of the model. A sinusoidal function with an amplitude of 20 mm was used for the displacement signal, and a sinusoidal function with an amplitude of 0~2 A was used for the control current signals. There were 1500 data points of which 1000 points were used for training and 500 points were used for validation. Figures 9 and 10 show the comparison between the experimental and predicted values as well as the error of the ANFIS inverse dynamics model of the MR damper, respectively. The results illustrate that the predictive control current of the ANFIS inverse model can track the target control current well, and verify the effectiveness of the model.

Elastic Model of Air Spring
The principle of the air suspension vehicle is that the air spring uses the air compressibility to play a supporting role. The vertical supporting force Fv produced by the air spring is the following: where P is the air pressure in the airbag, and A is the effective bearing area. The air spring moves vertically in the process of working, and the air spring bag is stretched or compressed, due to the change in volume. This process exchanges heat with the outside world. Assuming that the time required for the movement of the air spring is short enough, it can be regarded as an adiabatic process [31]. The state equation of the gas at this time can be regarded as follows: where V0 is the volume of the air spring in the initial state, V is the working volume of the air spring, P0 is the gas pressure in the initial state of the air spring, and Pa is the standard atmospheric pressure. According to Equation (13), the pressure in the air bag can be obtained as follows: Substituting Equations (13) and (14) to Equation (12), it can be derived as follows: Figure 10. Error of ANFIS inverse dynamics model.

Elastic Model of Air Spring
The principle of the air suspension vehicle is that the air spring uses the air compressibility to play a supporting role. The vertical supporting force F v produced by the air spring is the following: where P is the air pressure in the airbag, and A is the effective bearing area. The air spring moves vertically in the process of working, and the air spring bag is stretched or compressed, due to the change in volume. This process exchanges heat with the outside world. Assuming that the time required for the movement of the air spring is short enough, it can be regarded as an adiabatic process [31]. The state equation of the gas at this time can be regarded as follows: where V 0 is the volume of the air spring in the initial state, V is the working volume of the air spring, P 0 is the gas pressure in the initial state of the air spring, and P a is the standard atmospheric pressure.
According to Equation (13), the pressure in the air bag can be obtained as follows: Substituting Equations (13) and (14) to Equation (12), it can be derived as follows: When the air spring works, its volume and bearing area are related to the displacement, and thus, taking the derivative of Equation (15) with respect to the displacement, the stiffness k a of the air spring is yielded, and it can be regarded as follows: When the air spring is in an equilibrium state, there are the following relations: Therefore, the air spring stiffness k 0 under an equilibrium state is obtained as follows: Appl. Sci. 2021, 11, 10925 10 of 21 The natural frequency of the air spring is the following:

Modeling of Quarter Vehicle MR Semi-Active Air Suspension System
Passive suspension is mainly composed of a passive damper and spring with immutable stiffness, and its structure is relatively simple, as shown in Figure 11. The damping force and stiffness of passive suspension cannot be adjusted in real time with different road surface and vehicle parameters. Under a certain suspension relative speed, only a fixed suspension output damping force can be provided. Therefore, the vibration damping capacity of passive suspension decreases significantly, and its road adaptability is poor under a certain range of road surface conditions, especially on a bad road surface. Therefore, the air spring stiffness k0 under an equilibrium state is obtained as follows: The natural frequency of the air spring is the following:

Modeling of Quarter Vehicle MR Semi-Active Air Suspension System
Passive suspension is mainly composed of a passive damper and spring with immutable stiffness, and its structure is relatively simple, as shown in Figure 11. The damping force and stiffness of passive suspension cannot be adjusted in real time with different road surface and vehicle parameters. Under a certain suspension relative speed, only a fixed suspension output damping force can be provided. Therefore, the vibration damping capacity of passive suspension decreases significantly, and its road adaptability is poor under a certain range of road surface conditions, especially on a bad road surface.  Figure 12 shows the model of a quarter vehicle MR semi-active air suspension system, which is composed of a car body (sprung mass), wheel (unsprung mass), tire stiffness, damping and suspension stiffness.  Figure 12 shows the model of a quarter vehicle MR semi-active air suspension system, which is composed of a car body (sprung mass), wheel (unsprung mass), tire stiffness, damping and suspension stiffness. According to Newton's second law, the motion equation of the system is as follows: where mb and mw are the mass of the sprung body and the unsprung wheel, respectively, ka is the stiffness of the air spring, kt is the stiffness of the tire, xb is the displacement of the mass of the body, xw is the displacement of the wheel, xt is the displacement of the road surface excitation, and fd depicts the output damping force of the MR damper.
The system state variable can be denoted as follows: According to Newton's second law, the motion equation of the system is as follows: ..
where m b and m w are the mass of the sprung body and the unsprung wheel, respectively, k a is the stiffness of the air spring, k t is the stiffness of the tire, x b is the displacement of the mass of the body, x w is the displacement of the wheel, x t is the displacement of the road surface excitation, and f d depicts the output damping force of the MR damper.
The system state variable can be denoted as follows: The output variable can be depicted as follows: where ..
x b is the acceleration of the car body, x b − x w is the dynamic deflection of the suspension, and k t (x w − x t ) is the dynamic load of the tire.
The state space expression of the suspension system is as follows: where Figure 13 is the reference model of sliding model controller based on the knowledge of optimal control theory proposed by Karnopp [32], which has excellent robustness and excellent control performance. According to Newton's second law, the dynamic differential equation is as follows:

Reference Model of Sliding Mode Controller
where m b and m w represent sprung body mass and unsprung wheel mass, respectively; k and k t denote the spring stiffness and tire stiffness, respectively; x b1 , x w1 , and x t1 are the body displacement, wheel displacement and displacement of the road surface excitation, respectively; and c b1 and c w1 are the damping coefficients.  Figure 13 is the reference model of sliding model controller based on the knowledge of optimal control theory proposed by Karnopp [32], which has excellent robustness and excellent control performance. According to Newton's second law, the dynamic differential equation is as follows:

Reference Model of Sliding Mode Controller
where mb and mw represent sprung body mass and unsprung wheel mass, respectively; k and kt denote the spring stiffness and tire stiffness, respectively; xb1, xw1, and xt1 are the body displacement, wheel displacement and displacement of the road surface excitation, respectively; and cb1 and cw1 are the damping coefficients. The state variables are selected as follows: Substituting Equation (24) into Equation (23), it can be obtained the following: The state variables are selected as follows: Substituting Equation (24) into Equation (23), it can be obtained the following: where

Design of Sliding Mode Controller
The design principle of the sliding mode controller (SMC) is to make the actual controlled system track the motion of the reference system and produce the sliding mode in the error dynamics system between the actual controlled system and the reference system. The error vector e of the air suspension system is defined as the velocity error, displacement error and integral of displacement error, and its matrix form can be expressed as follows: . e = [ ..
Therefore, the error dynamics equation is defined as follows: In the design of switching the surface of the sliding mode variable structure, the pole assignment method is generally used to set the function of the switching surface. The output value of the controller is the controllable damping force, and thus its switching surface can be expressed as follows: Parameter n is 3, so the sliding mode switching surface function is as follows: The state Equation (29) is simplified as follows: All control points of sliding mode motion must have good dynamic quality and be asymptotically stable when entering the sliding mode switching surface. Equation (31) is the motion differential equation of the sliding mode, and its characteristic equation is as follows: The values of c 1 and c 2 can be obtained by making the characteristic root equal to the given pole.
The transfer function in standard form for a second-order system is as follows: The two closed-loop poles of the system are s 1,2 = −ξω n ± jω n 1 − ξ 2 . Only when the pole is on the left half of the complex plane does the control system have partial oscillations. The desired number of poles n is set as 3. Among them, the dynamic performance of the control system is determined by the dominant pole, and the far pole has little influence on it. It can be obtained that the characteristic root of the system is −3.9 ± 5.2j, under the condition of σ ≤ 15%, t p ≤ 0.7, ω n = 6.5, and ξ = 0.6. Then, c = [42. 25, 7.8, 1] can be obtained.
Therefore, the switching surface function of MR semi-active air suspension is as follows: According to Equation (30), it can be obtained as follows: .
From the above, it can be developed as follows: To minimize the error between the actual controlled model and the reference model, let . s be 0, and Equation (36) can be deduced as follows: The equivalent output damping force of the sliding mode controller can be obtained as follows: To improve sliding mode motion, the reaching law is defined as follows: where sgn(s) is a sign function, and ε is the rate of control point approaching the switching surface, which affects the system stability. The smaller ε is, the smaller the approaching rate is, while the larger ε is, the greater the approaching rate is, and the more severe the system chattering is. In the sliding mode control, ε is 0.1 according to the experience and debugging results. The system stability was verified according to the Lyapunov stability criterion. If V(x) is positive definite and .
V(x) is negative definite, the system is stable. Select the Lyapunov function as follows: Then it can be derived as follows: .
According to the Lyapunov stability criterion, it can be inferred that the system is stable. Therefore, the final sliding mode control law of the system is as follows: where m s is the coefficient related to s and .
s. Finally, the variable damping force by the sliding mode controller is as follows:

Design of Fuzzy Sliding Mode Controller
The structure diagram of FSMC is shown in Figure 14. In order to make the controlled system track the motion of the reference system, the sliding mode switching surface function s is designed according to the error between the controlled system and the reference system. The equivalent output damping force U eq of the sliding mode controller can be obtained. At the same time, ε and the switching control force U sw are obtained by the fuzzy controller. Finally, the force required for the damper (predictive damping force) can be obtained. According to the inverse model of the damper, the control current is obtained, and the actual control force is generated by controlling the MR damper.

Design of Fuzzy Sliding Mode Controller
The structure diagram of FSMC is shown in Figure 14. In order to make the controlled system track the motion of the reference system, the sliding mode switching surface function s is designed according to the error between the controlled system and the reference system. The equivalent output damping force Ueq of the sliding mode controller can be obtained. At the same time, ε and the switching control force Usw are obtained by the fuzzy controller. Finally, the force required for the damper (predictive damping force) can be obtained. According to the inverse model of the damper, the control current is obtained, and the actual control force is generated by controlling the MR damper.  In the fuzzy system, s and s  are defined as inputs of the fuzzy controller, and the output value of the fuzzy controller is defined as ε of the sliding mode controller. As shown in Figure 15  If the controller satisfies the conditions of < 0, the fuzzy control rules shown in Table 2 are established. The control target can be expressed as follows:

Design of Fuzzy Sliding Mode Controller
The structure diagram of FSMC is shown in Figure 14. In order to make the con system track the motion of the reference system, the sliding mode switching surfac tion s is designed according to the error between the controlled system and the ref system. The equivalent output damping force Ueq of the sliding mode controller obtained. At the same time, ε and the switching control force Usw are obtained by the controller. Finally, the force required for the damper (predictive damping force) obtained. According to the inverse model of the damper, the control current is ob and the actual control force is generated by controlling the MR damper.   If the controller satisfies the conditions of < 0, the fuzzy control rules sho Table 2 are established. The control target can be expressed as follows: If the controller satisfies the conditions of s . s < 0, the fuzzy control rules shown in Table 2 are established. The control target can be expressed as follows: where ε f uzzy is the output of the fuzzy sliding mode controller.

Random Road Excitation Response
In order to verify the effectiveness of the FSMC method, the random road excitation was selected as the road input for simulation analysis. The white noise signal after the forming filter was usually used as the time domain signal of the random road surface. The real road surface of vehicle driving can be regarded as follows: .
where f 1 is the lower cutoff frequency of the road filter, q(t) donates the road surface excitation, G q (n 0 ) is the road roughness coefficient, w(t) is the white noise with a covariance of 1, and v is the vehicle speed.
In the simulation analysis, MATLAB software was used for numerical simulation of the system, in which the sampling period is 0.001 s. Figure 16 shows the random input waveform of the road surface, using the filtering white noise generation method, where f 1 is 0.01 Hz, the C-level road roughness coefficient G q (n 0 ) is 256 × 10 −6 m 3 , and the vehicle velocity v is 20 m/s. Table 3 shows the parameters of semi-active air suspension system required for simulation.

Random Road Excitation Response
In order to verify the effectiveness of the FSMC method, the random road was selected as the road input for simulation analysis. The white noise signa forming filter was usually used as the time domain signal of the random road su real road surface of vehicle driving can be regarded as follows: where f1 is the lower cutoff frequency of the road filter, q(t) donates the road sur tation, Gq(n0) is the road roughness coefficient, w(t) is the white noise with a cov 1, and v is the vehicle speed.
In the simulation analysis, MATLAB software was used for numerical sim the system, in which the sampling period is 0.001 s. Figure 16 shows the rand waveform of the road surface, using the filtering white noise generation method is 0.01 Hz, the C-level road roughness coefficient Gq(n0) is 256 × 10 −6 m 3 , and th velocity v is 20 m/s. Table 3 shows the parameters of semi-active air suspensio required for simulation.    Figure 17 shows the time domain of the 2-DOF quarter vehicle semi-active air suspension controlled by the FSMC method. The random road excitations of passive suspension, SMC control and FSMC control suspension are the same, shown in Figure 16. It can be seen that the semi-active air suspension controlled by FSMC switching has significantly improved the damping performance, compared with the passive suspension and the semi active air suspension controlled by SMC.  Figure 17 shows the time domain of the 2-DOF quarter vehicle semi-active air suspension controlled by the FSMC method. The random road excitations of passive suspension, SMC control and FSMC control suspension are the same, shown in Figure 16. It can be seen that the semi-active air suspension controlled by FSMC switching has significantly improved the damping performance, compared with the passive suspension and the semi active air suspension controlled by SMC.      Figure 18 shows the power spectral density (PSD) under the input of random road surface. In the high frequency range, the power spectral density of FSMC control and sliding mode control is similar to that of the passive system in terms of body acceleration, suspension dynamic travel and tire dynamic load, due to frequency conversion interference. However, in the low frequency domain where the body is sensitive to the resonance of the car body, the 2-DOF quarter car semi-active air suspension controlled by FSMC is obviously superior to the controlled suspension and the sliding mode control suspension, the occupant's comfort is enhanced, and the resonance of the car body is effectively prevented.   Figure 18 shows the power spectral density (PSD) under the input of random road surface. In the high frequency range, the power spectral density of FSMC control and sliding mode control is similar to that of the passive system in terms of body acceleration, suspension dynamic travel and tire dynamic load, due to frequency conversion interference. However, in the low frequency domain where the body is sensitive to the resonance of the car body, the 2-DOF quarter car semi-active air suspension controlled by FSMC is obviously superior to the controlled suspension and the sliding mode control suspension, the occupant's comfort is enhanced, and the resonance of the car body is effectively prevented.

Impacted Pavement Excitation Response
Vehicle driving is not only the existence of random continuous vibration; sometimes, there will be low-lying, protrusion sections, which bring strong impact to the vehicle. Therefore, in order to fully reflect the vibration reduction ability of the vehicle suspension system, the impact road surface model was established to replace the impact vibration of the vehicle on the road surface. Figure 19 shows the simulation curve of impact road surface input, and its mathematical model is as follows: where x 0 (t) represents the input displacement of the impact road, H is the height of the impact road with a value of 0.07 m, W 0 represents the width of the impact road with a value of 0.08 m, u donates the vehicle speed, and the value is 3.08 km/h.

Impacted Pavement Excitation Response
Vehicle driving is not only the existence of random continuous vibra there will be low-lying, protrusion sections, which bring strong impac Therefore, in order to fully reflect the vibration reduction ability of the ve system, the impact road surface model was established to replace the im the vehicle on the road surface. Figure 19 shows the simulation curve of i face input, and its mathematical model is as follows: where x0(t) represents the input displacement of the impact road, H is t impact road with a value of 0.07 m, W0 represents the width of the imp value of 0.08 m, u donates the vehicle speed, and the value is 3.08 km/h. The simulation curve of the 2-DOF quarter vehicle semi-active air trolled by FSMC impacting the road surface is shown in Figure 20. It refle value of the FSMC-controlled semi-active air suspension is significantly pared with that of the passive and SMC-controlled semi-active air suspe same impact road excitation. The simulation curve of the 2-DOF quarter vehicle semi-active air suspension controlled by FSMC impacting the road surface is shown in Figure 20. It reflects that the peak value of the FSMC-controlled semi-active air suspension is significantly reduced, compared with that of the passive and SMC-controlled semi-active air suspension under the same impact road excitation.
Compared with passive suspension, the peak reduction percentages of body acceleration, suspension dynamic travel and wheel dynamic load of semi-active suspension under SMC control and FSMC control are shown in Figure 21, in which the body acceleration under SMC control and FSMC control are 13.9% and 27.5%, respectively, lower than that of passive suspension. The suspension dynamic travel is reduced by 12.8% and 26.5%, respectively, while the wheel dynamic load is reduced by 5% and 9.6%, respectively.
The simulation curve of the 2-DOF quarter vehicle semi-active air suspension con-trolled by FSMC impacting the road surface is shown in Figure 20. It reflects that the peak value of the FSMC-controlled semi-active air suspension is significantly reduced, compared with that of the passive and SMC-controlled semi-active air suspension under the same impact road excitation. Compared with passive suspension, the peak reduction percentages of body acceleration, suspension dynamic travel and wheel dynamic load of semi-active suspension under SMC control and FSMC control are shown in Figure 21, in which the body acceleration under SMC control and FSMC control are 13.9% and 27.5%, respectively, lower than that of passive suspension. The suspension dynamic travel is reduced by 12.8% and 26.5%, respectively, while the wheel dynamic load is reduced by 5% and 9.6%, respectively. Compared with passive suspension, the peak reduction percentages of body acceleration, suspension dynamic travel and wheel dynamic load of semi-active suspension under SMC control and FSMC control are shown in Figure 21, in which the body acceleration under SMC control and FSMC control are 13.9% and 27.5%, respectively, lower than that of passive suspension. The suspension dynamic travel is reduced by 12.8% and 26.5%, respectively, while the wheel dynamic load is reduced by 5% and 9.6%, respectively.

Conclusions
Based on the experimental data, the parameters of the improved hyperbolic tangent model were identified by the genetic algorithm, and the magnetic parameter identification results were consistent with the experimental data. The identified parameters of the improved hyperbolic tangent model could be used in semi-active control. The predictive control current of ANFIS inverse model can track the target control current well, which

Conclusions
Based on the experimental data, the parameters of the improved hyperbolic tangent model were identified by the genetic algorithm, and the magnetic parameter identification results were consistent with the experimental data. The identified parameters of the improved hyperbolic tangent model could be used in semi-active control. The predictive control current of ANFIS inverse model can track the target control current well, which verifies the effectiveness of the model.
The semi-active air suspension system model of quarter vehicle equipped with MR damper was established, and the feasibility of the FSMC strategy was verified by taking the random road excitation and impact road excitation as input. Compared with the passive suspension, the vehicle acceleration under SMC control and FSMC control is reduced by 15.9% and 26%, respectively; the suspension dynamic travel is reduced by 16.9% and 29.4%, respectively; and the wheel dynamic load is reduced by 7.4% and 8.4%, respectively, under random road excitation. With the impact of road excitation, the peak acceleration of the vehicle based on SMC control and FSMC control decreases by 13.9% and 27.5%, respectively; the peak dynamic travel of suspension decreases by 12.8% and 26.5%, respectively; and the peak dynamic load of tire decreases by 5% and 9.6%, respectively, compared with the passive suspension. It can be concluded that the control strategy proposed in this article has a better damping effect and can improve the driving comfort and stability of the vehicle.