Comparative simulation study on two estimation methods and two control strategies for semi-active suspension

With regard to the structural characteristics of the McPherson suspension system, when a vehicle is being driven on a rough road surface, the force direction of the suspension varies. This poses challenges to the vehicle’s driving safety and handling stability. Based on Lagrangian equations, this paper proposes a new nonlinear semi-vehicle suspension model and presents comparative studies, conducted through simulation, on the estimated accuracy and computational overhead of the small-computational-overhead extended Kalman filter (EKF) and unscented Kalman estimation (UKF) methods, and on the effectiveness of the skyhook sliding mode control (SHSMC) and nonlinear skyhook-sliding mode control (NSHSMC) semi-active suspension control methods. The response of the vehicle to the state estimation algorithm was evaluated through computer simulations using the Carsim vehicle dynamic software. The simulation results reveal that the vehicle dynamic states were satisfactorily estimated when the vehicle was driven on a rough road surface. Compared with the small-computational-overhead EKF algorithm, the estimated results of these variables based on the UKF algorithm have higher accuracy. However, the UKF algorithm requires longer computation time compared with the EKF algorithm. The SHSMC control algorithm achieved greater improvement for the vehicle’s drive handling stability in the 6–10-Hz vibration region compared with the NSHSMC control algorithm. In a high-frequency region over 10Hz, the semi-active suspension controlled by the SHSMC method had a more adverse effect on the driving comfort.


Problem statement
Vehicle suspension is a device that controls the forces and moments induced by the road roughness on the wheels and acting on the vehicle body. Different suspension types have different mechanical mechanisms. McPherson's independent suspension system is widely used as the front suspension of a vehicle. However, owing to this system's structural characteristics, when the wheel moves up and down, the lateral swing of the wheel occurs, and the suspension system cannot resist the lateral impacts. The traditional vector linear model only considers the upward and downward movement of its sprung and un-sprung mass, and the lateral motion characteristics of the McPherson suspension cannot be effectively characterized.
When a vehicle is being driven on a rough road surface, the vehicle suspension system plays an important role in the vehicle's ride comfort and handling stability. Notably, the ride comfort and handling stability are a pair of conflicting indices, that is, when the handling stability is improved, the ride comfort often deteriorates. Compared with passive suspension, semi-active suspension can be used to improve the handling stability without compromising the ride comfort.
Owing to the mechanical complexity of McPherson's independent suspension system, it is necessary to accurately obtain the system's mechanical dynamic states to realize the semi-active suspension control algorithm. With the least number of sensors, the states estimation algorithm can be used to estimate the mechanical signal, which is difficult to measure.

State of the art
Various studies have focused on vehicle dynamic models that include suspension mechanisms to investigate the dynamic performance of vehicles. Owing to the nonlinear geometrical effect of these mechanisms, the quarter-vehicle McPherson suspension nonlinear dynamic model has been investigated. [1][2][3][4][5] The abovementioned studies are based on the dynamic model of a quarter-vehicle suspension system, assumed that the vehicle sprung mass only performs vertical motion, and ignored the roll motion of the vehicle sprung mass.
Based on the linear vector vehicle suspension model, the ''skyhook damper'' scheme, whereby a virtual damper is set between the vehicle sprung mass and the inertial reference, has been proposed. 6 In recent years, based on the ''skyhook damper'' scheme, novel semiactive suspension control strategies and structures have been proposed. [7][8][9][10] Most modern vehicle control systems employ a feedback control structure. Therefore, the real-time estimation of the vehicle's dynamic states is required to enhance the performance of the vehicle's control system. Moreover, vehicle state estimation is becoming ever more relevant as chassis control systems are becoming increasingly sophisticated. 11 However, most existing studies on vehicle state estimation are based on the quarter traditional linear vector model, which is different to actual nonlinear suspension motion, and therefore results in estimation error. Furthermore, by comparing actual vehicles to the quarter-vehicle model used in the above-mentioned studies, it can be found that the system noise and measurement noise, which form owing to the wheel-road interaction on both sides of the vehicle, have an additive effect. Hence, the system noise and measurement noise variance must be re-determined.
In this study, by considering the mechanical structure characteristics of the McPherson suspension system, a semi-vehicle suspension dynamics model that considers the roughness of the road surface was established. Considering the superposition effect of the left and right road surface roughness excitation, the smallcomputational-overhead extended Kalman filter (EKF) and unscented Kalman estimation (UKF) state estimation algorithms of the semi-vehicle suspension dynamics model were developed. Additionally, the Semi-active suspension control algorithms, skyhook sliding mode control (SHSMC) and nonlinear skyhook sliding mode control (NSHSMC) were established. A road level classification algorithm based on the difference of dynamic states in the time and frequency domains has been proposed. 12,13 The road process noise variance matrices Q and measurement noise variance matrices R of the different levels of road roughness were provided to the above-mentioned Kalman filter algorithm. This study was structured as shown in Figure 1.

Contribution
(1) For the feasibility of implementation in onboard hardware, this paper discusses the estimated accuracy and computational overhead of the small-computational-overhead EKF and UKF state estimation methods for the vehicle suspension system when the vehicle is being driven on a rough road surface. (2) To determine the advantages, disadvantages, and application field of the NSHSMC and SHSMC control algorithms, the vehicle's drive handling stability and ride comfort with the semi-active suspension controlled by the NSHSMC and SHSMC control algorithms are discussed with regard to different frequency regions.
The rest of this paper is organized as follows. In Section 2, the wheel-road excitation model and a nonlinear semi-vehicle McPherson suspension model are established. In Section 3, the small-computationaloverhead EKF and UKF algorithms are constructed, and the system noise and measurement noise variance are re-determined. In Section 4, the new nonlinear control algorithms of semi-active suspension, namely, NSHSMC and SHSMC, are proposed and their stability is verified. Finally, the conclusions drawn from this study are presented in Section 5.

Road excitation
Road roughness function built using harmonic superposition method. As a source of system excitation, the road irregularities are analyzed using a typical ergodic stationary stochastic process. 14,15 The statistical characteristics of road irregularities are described by the power spectral density G q (n). According to GB7031, the fitting expression of the road power spectral density is expressed as follows: In equation (1), n 0 is the reference space frequency, and n 0 = 0.1 m 21 ; G q (n 0 ) is the road power spectral density under the reference space-frequency n 0 ; W is the frequency index, and W = 2. The range of the spatial frequency n is (n 1 , n 2 ).The road roughness variance is expressed as follows: The interval (n 1 , n 2 ) is divided into m small intervals. The power spectral density G q (n mid-i ) at the center frequency n mid-i (i = 1, 2, ÁÁÁ, m) of each interval replaces the spectral density values between the intervals. Equation (2) is approximately discretized as follows: For each interval, the sinusoidal wave function with frequency n mid-i (i = 1,2, ÁÁÁ, m) and standard deviation ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi G q (n midÀi )Dn i p is expressed as follows: The sinusoidal wave function corresponding to each interval is superimposed to obtain the road surface roughness function model, as follows: In equation (5), I is the length of the road (m); Dn is the frequency interval (m 21 ); u i is a random variable with [0, 2p] uniform distribution and independence. Then, equation (5), is converted to a time-domain model, as follows: In equation (6), t = I/u is time (s); f mid-i = n mid-i u is the time-frequency corresponding to n mid-i (Hz), where u is the velocity of the vehicle (m/s).
Coherence of bilateral track excitation. The left and right road tracks have the same road roughness statistical characteristics and self-power spectral density on the same road. However, the random process of the road roughness of the two tracks has a cross-spectrum. The two tracks are coherent and the coherence strength is weak. In equation (6), u i is the key parameter of the coherence of the left and right track; therefore, u i_right can be obtained by fitting. The coherence exponentially decreases as follows 16,17 : Where r is an adjustment parameter and r = 1. Then, u i_right can be calculated as follows:

New nonlinear semi-vehicle suspension model
In this study, a new nonlinear semi-vehicle suspension model was established for the Macpherson suspension system, as shown in Figure 2. The model includes a rigid chassis (1), control arm (6), wheel (5), damper (4), and coil spring (3) assemblies, and steering knuckle (2). The control arm is connected to the chassis with a rotational joint at point O and O#. The spherical joint connects the control arm at spindle B. The wheel is installed onto the steering knuckle. A cylindrical joint connects the steering knuckle to the damper, and the damper is connected to the chassis with a spherical joint at point A, as shown in Figure 2.
Because the mass of the control arm is much smaller than the mass of the spring and un-spring mass, the mass of the control arm can be ignored. The entire system has four degrees of freedom. In this case, the generalized coordinates are the upward vertical displacement z s , roll angle f of the vehicle sprung mass, rotation angle of the right control arm u r counter-clockwise, and rotation angle of the left control arm u l clockwise.
The semi-vehicle suspension dynamic model is constructed based on the following assumptions: (1) The un-sprung mass is linked to the chassis in two ways: through the damper and through the lower arm; u r and u l denote the angular displacement of the control arm.
(2) The sprung mass distribution coefficient e of the front and rear axles of the vehicles is equal to one; that is, the front and rear half of the vehicle move independently in the vertical direction. The longitudinal effect of the road profile and the effect of the vehicle pitching motion are ignored. (3) The bending and torsional deformation of the wheels and the effect of the wheel width are ignored.
Let z ul and z ur denote the vertical displacement of the center of mass of the left wheel C and right wheel C#, respectively, as follows: where z s is the vertical displacement of the vehicle body; O m is the point of the vehicle roll center; O and O# are the point of the mounting position of the rotational joint between the control arm and the vehicle chassis, respectively; r 0 is the initial angular displacement between O m O and the horizontal direction; l 0 is the distance between O m and O; l c is the distance between O and C.
Here, u 0 is the initial angular displacement of the control arm at the equilibrium state. Let a# = a 0 + u 0 . e`0 be the initial angular displacement between the control arm and line OA or O#A#. Then, the following relationships can be obtained from triangle OAB and O#A#B#: where l r and l l are the initial distance of A-B and A#-B# at the equilibrium state, and l r # and l l # are the changes in the distance of A-B and A#-B# with the rotation of the control arm by u r and u l ; l a is the distance between Then, the deflection of the spring can be expressed as follows: The relative velocity of the damper is expressed as follows: Equations of motion. The motion equations of the new nonlinear semi-vehicle suspension model are derived using the second Lagrange equation. Among them, T, V, and D are the kinematic energy, potential energy, and damping function of the system, respectively.
The state variables are introduced as q 1 = z s , q 2 = f, q 3 = u r , andq 4 = u l . Then, the dynamical model can be written in the state equation as follows: Then, equation (12) can be converted as follows: It is assumed that the vehicle system state vector is selected as follows: Then, the motion equations of the new nonlinear semivehicle suspension model can be written as follows.
The nonlinear system equation (15) can be written in vector form, as follows: Vehicle dynamic state estimator

Small-computational-overhead EKF method
Regarding the nonlinear discrete-time system expressed in equation (16), for the EKF, in each sampling time, the state transition matrix of the system is replaced by the Jacobian matrix of f at the k-step prediction value. Compared with the linear Kalman estimation method, the computational overhead is much larger. In this study, the system state variables fluctuated around the equilibrium state (x 1e , x 2e , x 3e , x 4e , x 5e , x 6e , x 7e , x 8e ) = (0, 0, 0, 0, 0, 0, 0, 0). The small-computationaloverhead EKF method is proposed to reduce the computational overhead. The highly nonlinear equations of the semi-vehicle suspension model are linearized at the equilibrium state, as follows: where C, E 1 , E 2 , D 1 , D 2 , and H are as described in the appendix.
Based on the above analysis, by combining the discrete-time state equations (21) and (22) for each state, the Kalman filter algorithm is reduced to the following expressions: 1006188, ð25Þ Equation (24) describes the prediction process from the estimated state at time k to the state at time k + 1. Equation (25) quantitatively describes the prediction process, where P (k/k) is the estimation error variance matrix, K (k + 1) is the Kalman gain, m 0 is the system initial state, and P 0 is the initial estimation error variance matrix.

Unscented Kalman estimation method
In Julier and Uhlmann, 19 the unscented Kalman estimation method (UKF) was proposed for nonlinear systems. The nonlinear system equation (16) is converted to a nonlinear discrete-time state equation, as follows: The measurement equation is expressed as follows: The UKF is expressed as follows: Initialize with: Calculate the sigma points: where l = a 2 L + k ð ÞÀL is a scaling parameter; a is set to 1e-3, k is set to zero, and b is set to 2.The relevant details are provided in Wan and Merwe. 20 Time update: Measurement update: K k + 1 ð Þ= P x k y k P À1 y k y k , ð43Þ Estimation of road process noise variance matrices Q and measurement noise variance matrices R In equations (25), (26), (39), and (41), the tuning operation of Q and R is a critical step for the accuracy of the estimation and control algorithm. These two parameters should be defined using a priori knowledge about the perturbationw k acting on the nonlinear system equation, andṽ k + 1 acting on the measurement equation; Q and R are the corresponding process noise variance and measurement noise variance matrices, respectively. Notably, Q and R are bounded positive definite matrices (that is Q . 0, R . 0). To obtain an optimal solution with the Kalman filter, it is assumed that the roughness of the road profile is Gaussian white noise. The noisew k is Gaussian white noise. Using the road excitation model proposed in Sections 2.1, the system process noise variance Q under different road excitation levels can easily be obtained, and the upper and lower spatial frequencies n 2 and n 1 are set as 2.83 m 21 and 0.011 m 21 .
Based on the above analysis, the Q q = E(w k w k T ) of the state estimator under different road levels(A-H) of the quarter vehicle are listed in Table 1.
The semi-vehicle model established in this study was compared to the quarter-vehicle model. While the vehicle is being driven, the wheels on both sides of the vehicle interact with the road at the same time, and there are multiple process noise w k and measurement noise v k + 1 sources of input to the vehicle during the driving process, owing to the superposition effect.
Let us assume that the distribution of the road surface roughness in all directions is isotropic Gaussian distribution. Then, the correlation function between two values of road surface roughness is a function only of the distance between them.
Let us assume that R q B ð Þ is constant under a certain level of road surface roughness.
The estimation algorithm processes noise variance Q, as follows: Q=diag q 1 , q 2 , q 3 , q 4 , q 5 , q 6 , q 7 , q 8 ½ : ð50Þ Many previous studies assumed that the value of R is constant and can be determined through empirical assessments for Kalman filter applications. However, the test equipment operates in a vibration environment with a rough road surface and the above-mentioned superimposed effects. The statistical characteristics of the measurement noise v k are affected and change, which increases the deviation between the estimated value and the actual value, and results in the divergence of the estimation algorithm. The method used to determine R in this study is described below. 18 In the measurement equation, the actual values of v k + 1 are unknown. Let us assume that Eṽ ð Þ =r and Eṽṽ T À Á = R. The measurement innovation is defined as follows:r Let us assume thatr k f g is a zero-mean white noise sequence, whose variance is expressed as follows: Where R k + 1 is the variance matrix for the measured noise, and P k + 1 is the estimation error variance matrix in the estimation method. Let us assume thatṽ k (k = 1, 2, ÁÁÁN) in the measurement equation are not correlated, andr and R are constant. Then, the sampling value of the measured noisẽ r k (k = 1, 2, ÁÁÁN)is an element of V R = {ṽ k }. Based on the noise samplingr k obtained in V R , the unbiased estimation ofr can be carried out as follows: The unbiased estimation of C r can be carried out as follows:Ĉ The expectation ofĈ r is expressed as follows: The estimation algorithm processes noise variance R k + 1 of the test equipment as follows: In Qin et al., 12 a road classification algorithm is proposed to classify the road roughness levels based on the time and frequency domain signals for different road roughness levels. In this study, this method was used to classify the uneven levels of the road online. When the vehicle is driven on different road levels, it provides real-time Q and R values to the estimator. The parameters considered in the investigation of the dynamic state estimation algorithm and control algorithm of the semi-active suspension system are listed in Table 2.

State estimation results
The effectiveness of the extended and unscented Kalman estimation methods proposed in Section 3.1 for a vehicle being driven on different road levels was verified. The system definition in Table 2 was adopted, and the vehicle velocity was set to 80 km/h (22.2 m/s). Based on the vehicle model constructed in Section 2.2, the vehicle was driven on ISO Level-B, -D, and -F roads. As presented in Table 3, the 3 3 1 measured signal vector [_ z s, k + 1 ;u lk ;u rk ] of the vehicle to the state estimation algorithm was evaluated using the vehicle dynamics software Carsim. The estimated results of the vehicle dynamic states, which were obtained while the vehicle was being driven on constructed roads, are presented in Figure 3.
Based on the small-computational-overhead EKF and UKF methods, as shown in Figure 3, when the   Table 3. Measured signal and estimated state of estimation algorithm. vehicle was driven on ISO Level-B, -D, and -F rough road surfaces, the actual dynamic state variables of the vehicle and the estimated results of these variables exhibited the same change pattern. The abovementioned vehicle dynamic states (a)z s , (b)u, (c) _ u, (d) _ u r , and (e) _ u l were satisfactorily estimated when the vehicle was driven on a rough road surface. Compared with the small-computational-overhead EKF algorithm, the estimation results for these variables based on the UKF algorithm were more accurate. However, the UKF algorithm required longer computation time for each step compared with the small-computationaloverhead EKF algorithm.

Measured signal Estimated state
For the proposed estimation method, the sampling time was set to 0.001 s, and the computation times of the small-computational-overhead EKF and UKF are shown in Figure 4.

Nonlinear reference mode sliding mode control algorithm
Karnopp proposed the ''skyhook damper'' scheme of passive suspension as a reference model to achieve the minimum , where x br and x w are the vertical displacements of the reference model sprung mass and un-sprung mass, respectively; E ( ) denotes the ''expected value.'' However, it is not possible to connect a damper from the isolated mass to an inertial reference. In this study, the control objective was to force the semi-active suspension system of the semi-vehicle to behave as the ''skyhook damper'' reference model by introducing an appropriate control input. The difference between the response provided by the reference model and the actual system was calculated using the sliding mode control (SMC) strategy.

Skyhook sliding mode control algorithm (SHSMC)
For a nonlinear semi-vehicle model, a semi-active suspension control algorithm is proposed by combining the ''skyhook'' reference mode and SMC strategy algorithms, based on the above Kalman state estimation method. The dynamical equation of the quarter vehicle sprung mass is expressed as follows: Equation (57) can be expressed in the form of a state equation, as follows: where m b is the sprung mass of the quarter vehicle model, x b and x w are the vertical displacement of the sprung mass and un-sprung mass, k s and c p are the suspension stiffness and damping coefficient, u is the semi-active suspension control force acting on the vehicle. In actual applications, it is impossible to realize the synchronous excitation of the reference model and actual model on the actual road surface. During the simulation, the reference Skyhook model is expressed as follows: Equation (58) can be expressed in the form of a state equation, as follows: where c sky-hook is the damping coefficient between the reference model sprung mass and the inertial reference.
The system tracking error is defined as follows: The system tracking error equation (59) can be rewritten in a state-space form, as follows: where: The system switching surface is defined as follows: where c 1 and c 0 satisfy Hurwitz's conditions. Equation (61) is differentiated once, as follows: To force the system dynamics to reach the sliding surface, the following system reaching law is used: To weaken the system chatter, when s is sufficiently small, the following relationships are used: At this time, e 1 is the movement speed gain when the states reach the sliding mode surface s. When the states pass through the sliding mode surface, a small e 1 can ensure that the approach speed is relatively small, and thereby small system chatter is ensured. However, because the system reaching law _ s consists of e 1 and e 2 , if e 1 is simply reduced, the overall approach speed outside the sliding mode surface will also be reduced. Therefore, when e 1 is reduced, e 2 should be appropriately increased to ensure the system's overall approach speed. In this study, e 1 = 1 and e 2 = 15 were selected. Then, the control force of the semi-active suspension controlled by the SHSMC can be expressed as follows: The control force of the semi-active suspension controlled by the SHSMC control algorithm was used as input to the vehicle dynamics software Carsim. According to the vehicle's structural geometry presented in Table 2, the SHSMC control algorithm's signal input x l b and x wl on the left side, and x r b and x wr on the right side, can be calculated as follows: x wl = z s À l oû + l c u l , ð67Þ x r b = z s + l aû , ð68Þ where x b in equation (57) corresponds to x l b in equation (66) and x r b in equation (68) of the semi-active suspension systems;x w corresponds to x wl in equation (67) and x wr in equation (69) of the semi-active suspension systems;û is the roll angular of the sprung mass. The Kalman state estimation method is used to estimate the dynamic state parameters of the vehicle for the control algorithm.

New nonlinear Skyhook-SMC control algorithm
This section proposes the new NSHSMC control algorithm, whose control structure is shown in Figure 5. According to the mechanical structural characteristics of the Macpherson suspension system, new quarter reference models and new quarter decomposition models of the nonlinear Macpherson suspension system are established on the left and right sides of the vehicle. As shown in Figure 5, there are four sub-models in total. In this study, the root mean square of the sprungmass acceleration rms € z s ð Þ and the ratio a w between the root mean square of the tire dynamic load and static tire load were adopted to further verify the ride comfort and handing stability of the suspension system, respectively, as follows: where F d is the tire dynamic load, and F s is the static tire load. In Table 4, for an ISO-E road surface, the ride comfort and handling stability indices of a vehicle with a passive suspension and semi-active suspension system  controlled by the SHSMC and NSHSMC control algorithm were obtained. Compared with passive suspension, the a w index of the semi-active suspension system controlled by the SHSMC control algorithm decreased, while the handing stability of the semi-vehicle improved. However, the value of rms € z s ð Þ increased, and the vehicle ride comfort decreased. The semi-active suspension system controlled by the NSHSMC control algorithm achieved better ride comfort compared with the SHSMC control algorithm and also improved the handling stability of vehicles.

Conclusion
In this paper, a nonlinear semi-vehicle model with a McPherson semi-active suspension system was developed using Lagrange equations. The small-computational-overhead EKF and UKF state estimation algorithms were developed, and the semi-active suspension control algorithms SHSMC and NSHSMC were established. The following conclusions were drawn from this study: (1) Based on the above-mentioned state estimation methods, the vehicle's dynamic states were satisfactorily estimated while the vehicle was being driven on a rough road surface.
Compared with the small-computationaloverhead EKF algorithm, the estimation of these variables based on the UKF algorithm was more accurate. However, the UKF algorithm requires longer computation time for each step compared with the EKF algorithm. (2) Compared with the NSHSMC control algorithm, in the 6-10-Hz vibration region, the SHSMC control algorithm achieved greater improvement of the vehicle's drive handling stability. In a high-frequency region greater than 10 Hz, the semi-active suspension controlled by the SHSMC control algorithm exerted a more adverse effect on the ride comfort.