Nonlinear Damping Curve Control of Semi-Active Suspension Based on Improved Particle Swarm Optimization

This paper studies the semi-active suspension damping control algorithm. The shortcomings of traditional damping control algorithms such as the sky-hook algorithm and the acceleration-driven damper algorithms are analyzed in this article. For the shortcomings of the traditional damping control algorithm, a semi-active suspension damping control strategy based on improved particle swarm optimization is proposed. The proposed algorithm uses the dynamic nonlinear inertia coefficient instead of the fixed inertia coefficient of the traditional particle swarm algorithm, which improves solution efficiency and accuracy. Nonlinear damping curves for different forms of expression are optimized using the algorithm, and the optimal nonlinear damping curve for the white noise pavement is obtained. A simulation model is established to verify the effect of the proposed algorithm. The simulation results show that the nonlinear semi-active control strategy can achieve the target of high vibration isolation coefficient at high frequency and strong inhibitory resonance ability at low frequency. Moreover, the jerk of the spring mass does not deteriorate by a large margin. The results show that the nonlinear semi-active control strategy improves the driving comfort of the vehicle.


I. INTRODUCTION
With the progress of technology, people have higher requirements for comfort during vehicle driving. The suspension system, as an important segment in the vibration reduction system, has become one of the research hotspots. According to the control mode, the suspension system can be divided into passive suspension, active suspension, and semi-active suspension. The stiffness and damping characteristics of the passive suspension cannot be adjusted and fixed in real time. It is difficult to adapt to the comfort requirements of different driving conditions. Although the active suspensions can attain a better effect, achieving large-scale promotion and application in the middle and low-end markets is difficult to due to its The associate editor coordinating the review of this manuscript and approving it for publication was Meng Huang . complex structure, large energy consumption, and high costs. Semi-active suspensions first appeared in 1973 [1], [2]. As an alternative to the high-cost active suspension system, the vibration isolation performance of the semi-active suspension is close to that of the active suspension and has the advantages of relatively simple structure, low energy consumption, low cost, and high reliability. It has become the current research and industrial application hotspot.
The damping control algorithm is the key to restrict the comfort of the semi-active suspension and has become the research focus of the semi-active suspension. Common classical control algorithms of the semi-active suspension are skyhook (SH) control [3], ground-hook (GH) control [4], [5], and acceleration-driven damper (ADD) control [6], [7]. This kind of control is based on an idealized model and adapts to single working conditions. For example, the SH algorithm has a good optimization performance at low frequencies; the ADD algorithm has a good optimization performance at high frequencies. On this basis, many semi-active control algorithms are based on SH control ideas, such as fractional SH, adaptive-SH [8], [9], [10], [11], and gain-scheduled SH [12]. The SH control algorithm is simple to realize and has great practical application value [13]. The hybrid algorithm of SH-GH, which combines the operational stability of the GH algorithm with the ride comfort of the SH algorithm, was proposed and experimentally verified by Ahmadian et al. [14], [15]. The modified ADD algorithm, also known as the SH-ADD algorithm [16], was proposed by Guo et al. The algorithm uses a frequency divider to control high and low frequencies separately: The low-frequency band uses the SH algorithm, and the high-frequency band uses the ADD algorithm. It combines the advantages of the two algorithms and achieves good optimization results. However, such algorithms only adjust at the limit position of the damping and do not make full use of the performance of the stepless adjustable dampers. Its frequent switching between high and low damping will lead to a very large jerk (the differential of acceleration) and the appearance of a ''tremor'' phenomenon. It will bring drivers a poor driving experience [17]. In recent years, semi-active suspension control methods have gradually emerged, such as model predictive control based on modern control method [18], [19], [20], [21], fuzzy control algorithm [22], [23], [24], [25], [26], and control algorithm based on the genetic algorithm and neural network optimization [27], [28]. Such methods have good control effect, but the algorithm structure is complex, which requires more computing resources and high cost.
The nonlinear vibration isolation system can improve the vibration isolation performance of the system itself by using its own stiffness nonlinearity and damping nonlinearity [29], [30], [31], [32]. In the nonlinear vibration isolation system, the optimal parameters of stiffness and damping can be designed according to its own characteristics and input excitation to achieve the target of high vibration isolation coefficient at high frequency and strong inhibitory resonance ability at low frequency. For semi-active suspension, the damped adjustable shock absorber can adjust a variety of nonlinear damping curves according to its own relative speed in the damping limit, which is suitable for multiple pavement conditions. Generally, in order to suppress the resonance of system, the nonlinear damping curve has a large damping coefficient when the relative speed is low, and it has a small damping coefficient to improve the vibration isolation coefficient when the relative speed is high [33]. The damping curve can be expressed as damping force versus relative velocity (F − v d ) or damping coefficient versus relative velocity (C − v d ). The C − v d curve expression is used as an example, as shown in Fig. 1.
The advantages and disadvantages of the classical control algorithms, SH and ADD algorithms, are analyzed. For the shortcomings of the SH and ADD algorithms, this paper proposes a nonlinear damping curve control strategy based on improved particle swarm optimization (PSO). Through the improved PSO, the damping curve under many typical roads is optimized to achieve the target of high vibration isolation coefficient at high frequency and strong inhibitory resonance ability at low frequency. Moreover, this algorithm eliminates the ''tremor'' phenomenon of the switching algorithm. In addition, the algorithm only needs the car's high sensor in the engineering application. It is great for reducing the cost problem of the semi-active suspension system application.

II. MODELING OF QUARTER-SUSPENSION SYSTEM A. KINETIC MODEL OF QUARTER SUSPENSION SYSTEM
The quarter-suspension model is shown in Fig. 2 (a). This model is a two degree of freedom simplified model. Its structure is simple and can qualitatively reflect the movement characteristics of the vehicle in the vertical direction. Therefore, it is commonly used in the preliminary simulation of the semi-active suspension control algorithm.
A kinetic analysis of this model is performed. The kinetic model is obtained as follows: For easy analysis, the force generated by the damper is separated because the damping can be changed.
The formula can be rewritten as follows: where m c is the sprung mass; m f is the under-sprung mass; c s is the damping factor of the passive suspension; k s is the stiffness coefficient of the damper spring; k t is the stiffness coefficient of the tire; f c is the damped force of the damper output, which replaces the damping force of the semi-active suspension in the model; z r , z f , and z c are the displacement of the road, under-sprung mass, and sprung mass, respectively. The specific parameters are listed in the following Table 1:

B. MODELING OF ROAD SURFACE
According to the file ISO/TC108/SC2N67 and GB7031-86, this white noise method is used to build the modeling of road surface. Its equation is as follows: where n 00 is the frequency of the lower cutoff space, u is the speed of the vehicle, z r (t) is the displacement of the random pavement height, n 0 is the frequency of the reference space, G q (n 0 ) is the coefficient of road unevenness, and W (t) is the Gaussian white noise, whose mean is zero. This model agrees with the standard pavement spectrum and clear physical importance, which can be used as an input incentive for vehicle ride comfort analysis [34].

III. CLASSICAL ALGORITHM ANALYSIS
The input signal is a separate frequency, and the under-spring mass velocity is as follows: The velocity and acceleration of the sprung mass are as follows:ż where A is the amplitude of the velocity of the under-sprung mass; |H (ω)| and ϕ(ω) are the amplitude-frequency characteristic function and phase-frequency characteristic function from under-sprung mass to sprung mass. Laplace Transform is used to transform (1) to obtain the following formula:ż ċ z f = c s s + k s m c s 2 + c s s + k s (8) The frequency characteristic is as follows: The bode plots of this system are shown as Fig. 3.

A. SKY-HOOK ALGORITHM
The SH algorithm is one of the classical algorithms. It has a good optimization effect in the low-frequency band but deteriorates in the high-frequency band. The causes of this phenomenon are analyzed. The formula of the SH algorithm is as follows: whereż def =ż c −ż f . Formulas (5), (6), and (10) are integrated to obtain the following: Upon further simplification, the following formula is obtained: Thus,ż def ·ż c > 0 can be equally replaced with the following formula: It can be discussed based on the frequency range of the input signal. The frequency range is divided into three levels.
When the frequency of the input signal is much smaller than the resonance frequency, the amplitude frequency characteristic is close to 1, and the values of the phase lag to 0. The formula is as follows: Formulas (12) and (14) are integrated to obtain the following:ż The state of damping is at the softest state. The frequency of the input signal cannot be zero, so the state of damping switches between the softest and the hardest. The effect on comfort is minimal at this time.
When the input frequency and resonance frequency of the signal are close, the amplitude frequency characteristic is greater than 1, and the values of the phase lag to −45 • . The formula is as follows: Formula (16) is incorporated into Formula (13) to obtain the following: When |H (ω)| = 1, the large state of damping is taken for 62.5% of a period. The proportion of large damped states increases as |H (ω)| improves. This situation corresponds with the expectation that the damper state takes the large damped state in the resonance frequency. The SH algorithm has a good optimization effect in the resonance region.
When the frequency of the input signal is much larger than the resonance frequency, the amplitude frequency characteristic is less than 1, and the values of the phase lag to −90 • . The formula is as follows: Formula (18) is incorporated into Formula (13) to obtain the following: When |H (ω)| = 1, the large state of damping is taken for 75% of a period. This situation is not consistent with the expectation of the damper state to take a small damping in the high-frequency range. Therefore, the SH algorithm has a poor optimization effect in the high-frequency range.
In conclusion, the SH algorithm has a good optimization effect in the low-frequency range but a poor optimization effect in the high-frequency range.

B. ACCELERATION DAMPING ALGORITHM
In contrast to the SH algorithm, the ADD algorithm has good optimization results in the high-frequency range. However, in the low-frequency range, it is not as effective as the SH algorithm. The formula of ADD algorithm is as follows: Similarly, Formulas (5), (7), and (20) are integrated to obtain the following: Thus,ż def ·z c > 0 can be equally replaced with the following formula: Again, it can be discussed based on the frequency range of the input signal.
When the frequency of the input signal is much smaller than the resonance frequency, the amplitude frequency characteristic and the values of the phase lag are as Formula (14). VOLUME 10, 2022 Formula (14) is integrated into Formula (22) to obtain the following: The ADD algorithm is similar to the SH algorithm in the low-frequency band range.
When the input frequency and resonance frequency of the signal are close, the amplitude frequency characteristic and the values of the phase lag are as Formula (16). Formula (16) is integrated into Formula (22) to obtain the following: When |H (ω)| = 1, the large state of damping is taken for 12.5% of a period. Although this proportion decreases when the values of |H (ω)| increase, the large state of damping is taken for 10% of a period when |H (ω)| = 10. This situation is not consistent with the expectation of the damper state to take a large damping in the low-frequency range. Therefore, the ADD algorithm has a poor optimization effect in the lowfrequency range.
When the frequency of the input signal is much larger than the resonance frequency, the amplitude frequency characteristic and the values of the phase lag are as Formula (18). Formula (18) is integrated into Formula (22) to obtain the following: When |H (ω)| = 1, the large state of damping is taken for 25% of a period. This proportion decreases as the value of |H (ω)| decreases. This situation corresponds with the expectation that the damper state takes the small damped state in the high-frequency range.
In conclusion, the ADD algorithm has a good optimization effect in the high-frequency state and a poor optimization effect in the low-frequency state.

IV. NONLINEAR DAMPING CURVE CONTROL
The optimization effects of the SH and ADD algorithms for different frequencies are analyzed from a statistical perspective above. Other algorithms, like SH and ADD algorithms, switch at the maximum and minimum damping quickly and produce a serious ''tremor'' phenomenon. The nonlinear damping curve control is presented in this paper. The damping curve of the nonlinear damping algorithm is continuous, so it can effectively eliminate the ''tremor'' phenomenon. This optimal damping curve can be obtained through PSO.

A. IMPROVED PARTICLE SWARM OPTIMIZATION
The idea of the PSO algorithm stems from the predation behavior of birds. The algorithm sets multiple particles. Each particle needs to optimize the parameters. Each particle shares location information, namely, using their own current parameters of the target to determine which particles are closest to the ''food,'' that is, the closest to the optimal solution.
Other particles are in the position of the particles closest to the optimal solution. In the process of approach, positions closer to the optimal solution than before may be found. This situation results in all particles replacing the target and moving toward positions with a better target. This movement is also called social learning and individual learning. The core of this algorithm is to iterate through continuous learning until the true ''optimal solution'' is reached. The conditions for the algorithm ends can be defined as follows: • The specified maximum number of iterations is retained. • The results remain unchanged after multiple iterations. • Data optimized for multiple times are within the specified error range. The advantages of this algorithm are its fast search speed, high optimization efficiency, and simple algorithm. It is suitable for finding the optimal damping curve under specific working conditions. The velocity and position update formula of the particle group algorithm are as follows: where n is the number of particles, v is the velocity of the particles, xis the position of the particles, w is the weight coefficient, c 1 , c 2 is the learning factor, r 1 , r 2 is a random number from 0 to 1, pbest is the optimal location of the individual, and gbest is the global optimal location. The velocity update formula shows the following: The value of the weight coefficient increases the effect of the current particle velocity on the next moment, so the weight coefficient is also known as the inertia factor. When the parameters are too large, the search space of the particles expands because the movement direction is difficult to change. The overall convergence rate of the algorithm is greatly reduced. Conversely, when the parameters decrease, the particle speed can easily change. The search space is smaller, and the overall convergence rate of the algorithm is relatively higher. However, when the range of particles does not include the optimal solution, target confusion or even absence of a target can occur, and the convergence rate of the algorithm is affected. Generally, the weight coefficients range from 0 to 1.
The above discussion reveals that the weight coefficient affects the optimization effect and convergence speed of the algorithm. Therefore, the weight coefficient of the next optimization adopts the adaptive value. At the beginning of the algorithm, the goal is to increase the particle search range and avoid the particles falling into the local optimal solution, so the weight coefficient should take a larger value. When the algorithm runs later, the algorithm should pursue a faster convergence rate, so the weight coefficient should take a smaller value. That conventional adaptive algorithms allow the weight coefficient to decrease with the number of iterations linearly is unreasonable. The weight coefficient should be slow in the early stage to improve the global search  ability. The late decline speed should be faster to accelerate the overall convergence speed of the algorithm. The adaptive algorithm for the weight coefficients is optimized as the firstorder inertia link of the inversion, with the specific formula as follows: where t is the number of current iterations; w t is the current weighting coefficient; w ini is the initial weight coefficient; w des is the target weight coefficient; T is the weight coefficient convergence time, 1/4 times the iterations of the weight coefficient change. When the algorithm is initialized, t = 0, and w t = 0.982w ini + 0.018w des ≈ w ini . When the number of iterations reaches the desired number, t = 4T , and w t = w des . If expected, after the 100 iterations, the weight coefficient is reduced from w ini to w des , and parameter T should be set to 25. The change of the weight coefficient is shown as Fig. 4.
The optimization of the improved PSO algorithm is as follows: 1. The particle population is initialized, and each particle obtains the initial position and velocity. 2. Fitness and global optimal fitness values are calculated for each particle.   (27). 5. The position and velocity of the current particle are updated using the particle state and the coefficient of inertia calculated in step 4. 6. Whether the algorithm end condition is met is determined. 7. If the algorithm meets the conditions, the results are output, and the algorithm ends. Otherwise, return to step 2. The flow chart of the algorithm is shown in Fig. 5.

B. SELECTION OF OPTIMIZED TARGETS AND OBJECTS
The target raises the riding comfort of the sprung mass. The ISO2631-1-1997 (E) presents the calculation criteria for ride comfort [35]. This standard document takes into account the impact of vibration at different frequencies, different directions and different positions on ride comfort. The optimization object is a damping curve. The damping curves have multiple forms of expression. The damping curves have three common expression forms to optimization.

1) FORM OF EMPIRICAL FORMULA
The empirical formula of damping curve is as follows: where v D is the relative velocity of the vibration damper, c D is the damping coefficient of the vibration damper, and β i (i = 1, 2, 3, 4) are the parameters to optimize. The optimization parameters are divided into compression and recovery stroke parameters due to the damping curve asymmetry corresponding to the compression and recovery stroke. The optimization parameters have upper and lower limits. On this basis, the damping coefficient of the semi-active suspension is limited, so the damping coefficient calculated by the empirical formula needs to be limited by the maximum minimum of the  Table 2.
The shape of the damping curve depends on the parameters β i . The magnitude of β 3 determines when β 3 v D cones into play. These parameters will work when the relative speed is large popularly. When the speed is not large, the previous one in the empirical formula plays a major role. When the speed is small, β 0 v D accounts for the main contribution, and as the speed gradually increases, the proportion of β 1 + β 2 v D gradually increases.

2) FORM OF A MULTIPOINT INTERPOLATION
The form of a multipoint interpolation is the shape of the entire damping curve using the damping coefficient at a specific relative velocity. The optimal curve is obtained by optimizing the damping coefficient at a specific relative velocity. The relative speeds of the vibration absorber used in engineering are ±0.052, ±0.131, ±0.262, and ±0.393. Therefore, this parameter is selected according to this standard as Table 3: c v i represents the damping coefficient at the relative velocity of the shock absorber is v i .(v i = ±0.052, ±0.131, ±0.262, and ±0.393)

3) FORM OF POLYNOMIAL
The form of polynomial is taken using polynomials to formulate the shape of the entire polynomial curve. The optimal curve is obtained by optimizing the polynomial equation. To ensure the same number of parameters as the above two methods, the highest time of the polynomial is 3. It produces three poles and can realize the expression of multiclass damping curves. The polynomial formula is expressed as follows: where a i (i = 0, 1, 2, 3) are the coefficients of the polynomials and the parameters to optimize next. The compression and stretching states are also subdivided into stretching coefficient and compression coefficient because they require different curves.

C. OPTIMIZATION PROCESS AND RESULTS
To verify the effectiveness of the optimization algorithm, and to compare the optimization effect of the nonlinear damping curves for the three different forms of expression. The white noise pavement conditions were used for simulation optimization. In the real vehicle ride comfort test, according to the experience of the company, 50 km/h is the most representative working condition, and the real vehicle data is the closest to the simulation results of class D white noise pavement. In order to better guide the real vehicle test through simulation, the input signal selects the Class D white noise pavement at a vehicle speed of 50 km/h (G q (n 0 ) = 1024 × 10 −6 m 3 ) and is optimized using an improved particle group algorithm.   There are some key parameters in the improved PSO algorithm need to be discussed. n is the number of particles, too much number of particles will lead to excessive computation in one iteration and reduce the iteration speed, and too small the number of particles will lead to a small search range and easy to fall into the local optimal solution. c 1 , c 2 are the individual and global learning factor, the range of this parameter is 0∼4 generally. w ini , w des are the initial and expected weight coefficients, the range of this parameter is 0∼1 generally. The improved PSO algorithm parameters are set in Table 4.
Three damping curves are iterated 500 times with the same parameters of the improved PSO algorithm. The results are shown in Fig. 6. VOLUME 10, 2022 When the optimization object is a multipoint interpolation method curve, the optimization effect and the sum convergence rate are good. It is stabilized by about 60 iterations. Most of the parameters in the polynomial curve cannot directly react to the adjustment of the damping curve. The various parameters in the damping curve of the multipoint interpolation can directly correspond to the shape of the damping curve, which facilitates the convergence of the PSO algorithm.
In conclusion, the expression form of the damping curve selects the multipoint interpolation form. The damping curve is optimized when the input signals are different working conditions. The optimization results are shown in Fig. 7.

D. APPLICATION TO THE WHOLE VEHICLE MODEL
Compared with the quarter suspension model, the whole vehicle model increases the movement of pitch and heel direction. This is closer to the real car driving situation. Take the whole vehicle comfort in the standard ISO2631 as the target function. The nonlinear damping curve control strategy is extended to the whole vehicle model. And the effect is also verified through simulation.
The 12 degrees of freedom dynamics model of the whole commercial vehicle is shown in Fig. 8. In this model, the chassis suspension damping is passive, and the cabin suspension damping is variable.
The same method is applied to separate the damping force of the cabin suspension. The kinetic equations are obtained as follows: Vertical vibration differential equations of the six chassis suspension: (30) where m wi is the mass of chassis suspension, z wi is the displacement of the chassis suspension, k wi is the stiffnesses coefficient of the tire, q i is the displacement of road, k si is the stiffnesses coefficient of the chassis suspension, z fi is the displacement of the sprung mass of the chassis suspension, c si is the damping coefficient of the chassis suspension. Three degrees of freedom vibration differential equation of the chassis: where m f is the mass of chassis, z f is the displacement of the chassis, z bi is the displacement of the under-sprung mass of  the cabin suspension, z ci is the displacement of the sprung mass of the cabin suspension, f i is the damping force of the cabin suspension, I fxx and I fyy are the moment of inertia of the chassis, k ci is the stiffnesses coefficient of the cabin suspension, α f and γ f are the Angular displacement of the chassis heel and pitch, respectively. (P fxi , P fyi ) and (P bxi , P byi ) are the coordinate of the chassis suspension and the cabin suspension in the coordinate system with the chassis centroid as the origin, respectively. Three degrees of freedom vibration differential equation of the cabin: where m f is the mass of cabin, z c is the displacement of the cabin, I cxx and I cyy are the moment of inertia of the cabin, (P cxi , P cyi ) is the coordinate of the cabin suspension in the coordinate system with the cabin centroid as the origin, α c and γ c are the Angular displacement of the cabin heel and pitch, respectively. The kinematic relationship of the three directions of the chassis and cabin to each suspension is as follows: 90966 VOLUME 10, 2022 The numerical values of model are shown in Table 5 and  Table 6.

V. SIMULATION ANALYSIS A. THE QUARTER SUSPENSION SIMULATION
The semi-active suspension model is constructed in Simulink according to the parameters in Table 1. The grade D white noise pavement at 50 km/h is used as input. The simulation  time is 50 s. The acceleration and urgency of the spring mass are analyzed. The effects of passive suspension, SH algorithm control, ADD algorithm control, and nonlinear damping curve control are compared. The passive suspension takes the maximum and minimum damping states as a reference. The effect comparison of spring mass acceleration for different strategies is shown in Fig. 9.
The jerk values of spring mass under different strategies are shown in Fig. 10.
The spectrum map of level D random pavement is messy. To compare the effect of using a nonlinear damping curve control strategy easily, it is optimized again using the sweep as the input pavement. The simulation results are also compared with other algorithms, as shown in Fig. 11.
Because there is only vertical motion in the simulation of quarter suspension model, only the effects of different frequencies are considered. Table 7 shows the comparison of acceleration weighted RMS and jerk peak values for different control strategies.
Compared with the SH algorithm and the ADD algorithm, the nonlinear semi-active control strategy can achieve the target of high vibration isolation coefficient at high-frequency state and strong inhibitory resonance ability at low-frequency state under the condition that the jerk of the spring mass does not deteriorate by a large margin.

B. THE WHOLE VEHICLE SIMULATION
The whole vehicle model is constructed in Simulink according to the parameters in Table 5 and Table 6. The grade D white noise pavement at 50 km/h is used as input. The simulation time is 50 s. The effects of passive suspension, SH algorithm control, ADD algorithm control, and nonlinear damping curve control are compared, the results are shown in Table 8. The results show that the nonlinear damping curve control can still achieve a good effect in the whole vehicle simulation.

VI. CONCLUSION
To improve vehicle ride comfort, the paper first analyzes the principles, advantages, and disadvantages of the SH algorithm and the ADD algorithm. An algorithm for nonlinear damping control based on the improved PSO is proposed. For the PSO used, the dynamic weight coefficient is introduced, improving the early search ability of the particle group algorithm and accelerating the late convergence rate of the algorithm. The damping curves are obtained and optimized by the improved PSO algorithm.
Ignoring the influence of high-order dynamics, the quarter suspension model and the 12 degrees of freedom dynamics model of the whole commercial vehicle is established, and the grade D white noise random pavement and sweep pavement are adopted for the input pavement. The simulation results from the passive suspension, SH algorithm, and ADD algorithm reveal that the nonlinear semi-active control strategy can achieve the target of high vibration isolation coefficient at high-frequency state and strong inhibitory resonance ability at low-frequency state under the condition that jerk of the spring mass does not deteriorate by a large margin.
The simulation results show that the nonlinear damping curve control based on the improved PSO has good results in improving vehicle ride comfort. The proposed algorithm requires relatively minimal hardware compared with other algorithms. It only needs the car's high-level sensor to obtain the relative speed. When the algorithm is applied to practice, it can effectively reduce the cost of hardware and expand the scope of application. His research interests include semi-active suspension, mechanical kinetics, and its application.