# Estimating Vehicle Stability Region Based on Energy Function

Authors: Yu-guang Yan; Hong-guo Xu; Hong-fei Liu
Journal: Discrete Dynamics in Nature and Society
 (2015)

---

## Abstract

In order to improve the deficiency of vehicle stability region, according to vehicle nonlinear dynamic model, method of estimating vehicle spatial stability region was proposed. With Pacejka magic formula tire model, nonlinear 3DOF vehicle model was deduced and verified though vehicle test. Detailed detecting system and data processing were introduced. In addition, stability of the vehicle system was discussed using Hurwitz criterion. By establishing energy function for vehicle system, the vehicle’s stability region in 20 m/s was estimated based on Lyapunov theorem and vehicle system characteristics. Vehicle test in the same condition shows that the calculated stability region defined by Lyapunov and system stability theorem has good effect on characterized vehicle stability and it could be a valuable reference for vehicle stability evaluation.


---

## Body

## 1. Introduction

It is clear that stability region for vehicle stability evaluation has been extensively studied in our country and abroad now. Up to date, some researchers have analyzed stability region applying simple vehicle model [1, 2]. Typically, in the 1990s, Inagaki, Samsundar, Yasui, and Stotsky proposed analytical expression for vehicle plane stability region based on simplified cubic tire model and deduced plane oval stability region by simulation [3–6]. In 2002, Ko and Lee estimated plane stability region of vehicle by side slip angle and yaw rate according to topology theory and 2DOF vehicle model [7]. Similarly, in 2007, Shen et al. estimated singular points of vehicle and studied stability characteristics of vehicle [8].Many reference classical conclusions play a positive role for vehicle stability region; however there are still some limitations: first, most of the traditional studies were based on simplified 2DOF vehicle models, so vehicle state could not be comprehensively considered. In addition, it is crystal clear that traditional vehicle dynamic model with simplified tire formula, such as linear tire model, was conveniently calculated but it still had great deviation compared with fact, and as a result, the accuracy of vehicle region was unsatisfied. Furthermore, traditional researches for vehicle region boundary conditions had insufficient basis, which affect the accuracy as well.Therefore, an improved nonlinear vehicle dynamic model with Pacejka magic formula is introduced, and the stability region of the vehicle is estimated based on energy method. It is valuable for evaluating the stability of vehicle and improving traffic safety.


## 2. Vehicle Nonlinear Dynamic Model

### 2.1. Vehicle Differential Equations

Figures1(a) and 1(b) illustrate the 3DOF vehicle model or double-axis vehicle and coordinate systems. Suppose the vehicle is steering at a speed without braking and driving forces, pitch movement, and longitudinal tire force and align toque could be ignored in this condition. So the established vehicle 3DOF model consists of lateral motion, yaw, and roll motions.Figure 1
Schematic diagram of body movement.
(a)
Plane movement
(b)
Roll movementCorresponding equations of the vehicle model are described as follows:(1)m(v˙y+vxω)-mshsψ¨s=Fyf+Fyr,Izω˙=aFyf-bFyr,Ixsψ¨s+mshsv˙y+vxω-hsψ¨s=msghsψs-Mψ,(2)Mψ=Mψf+Mψr,Mψf=Kfψs+Dfψ˙s,Mψr=Krψs+Drψ˙s.


### 2.2. Nonlinear Tire Model

The lateral tire force depending on the side slip angle is represented as the Pacejka magic formula [9]:(3)Fy=DsinCarctanBα-EBα-arctanBα,where α is the slip angle of tire and B, C, D, and E are, respectively, stiffness factor, shape factor, peak factor, and curvature factor. These factors are obtained by fitting. According to the structure of vehicle tire and vertical load, the relevant factor values of tire on normal road (μ=0.7) are listed (Table 1) [10].Table 1
Curve factors’ values of nonlinear tire.
Tire position (number)
B
C
D
E
Front wheel (2)
10.30
1.30
−11821.63
−1.05
Rear wheel (4)
10.51
−13437.49


### 2.3. Vehicle Dynamic Model Validation

Vehicle test on normal road is to verify the vehicle dynamic model. The details of tests environment and related equipment are showed (Figure2).Figure 2
Details of vehicle tests and related equipment.As shown in Figure2, the integrated test system mainly involves HCZ-1 vehicle handling force and angle measuring instrument, inertial sensors IMU02, vehicle motion parameters measuring instrument by GPS called RLVB20SL, data collection device RLVB-ADC03, and computer. Function of the system is to measure the steering wheel angle input and vehicle motion parameters, such as longitudinal, lateral, and vertical velocity, acceleration, angle, and angular velocity.The test data is processed by VBOX Tools and Matlab. The position deviations between the actual locations of the test equipment and the vehicle’s centroid should be taken into account, so processing the sampling data is necessary and a key step [11].(1) Inertial Sensors IMU02 Direction Correction. Coordinate system of IMU02 (Figure 3) is corrected by right-hand rule and detail (Table 2).Table 2
Coordinate system of IMU02 correction.
Parameter
Longitudinal
Lateral
Vertical
Yaw
Roll
Pitch
Direction correction
+
+
−
−
+
−Figure 3
IMU02 coordinate.(2) Vehicle Centroid State Parameters Correction. According to the deviation between the actual locations of the test equipment and the vehicle’s centroid, sampling data of RLVB20SL and IMU02 should be processed as follows.(1) Vehicle Longitudinal Velocity and Acceleration. Consider(4)vx=vxGPS-sz1γ˙INS,ax=axINS-sz2γ¨INS,where vxGPS is longitudinal velocity tested by RLVB20SL, axINS andγINS are longitudinal acceleration and pitch angle tested by IMU02, sz1 is vertical distance between the main antenna of RLVB20SL and vehicle centroid, and sz2 is vertical distance between the location of IMU02 and vehicle centroid.(2) Vehicle Lateral Velocity and Acceleration. Consider(5)vy=vyGPS-sx1ωINS+sz1ψ˙INS,ay=ayINS-sx2ω˙INS+sz2ψ¨INS,where vyGPS is lateral velocity tested by RLVB20SL, ayINS, ωINS, and ΨINS are lateral acceleration, yaw rate, and roll angle tested by IMU02, sx1 is longitudinal distance between the main antenna of RLVB20SL and vehicle centroid, and sx2 is longitudinal distance between the location of IMU02 and vehicle centroid.As shown in Figure2, the tested vehicle is driven along the double change race on normal road [12]. The velocity and steering wheel angle operated by driver are showed in Figure 4(a). Then, under the same condition, vehicle lateral acceleration by vehicle tests and Matlab simulation applying the established model is obtained, as shown in Figure 4(b).Figure 4
Comparing simulation with vehicle test.
(a)
Vehicle’s velocity and steering wheel angle input
(b)
Vehicle’s lateral acceleration outputClearly, vehicle’s lateral acceleration in Figure4 shows the good consistency between simulation and vehicle test, and the error is less than 10%. So the reliability of the established nonlinear dynamic model is definitely concluded and could be used for the vehicle stability analysis described in the next section.


## 2.1. Vehicle Differential Equations

Figures1(a) and 1(b) illustrate the 3DOF vehicle model or double-axis vehicle and coordinate systems. Suppose the vehicle is steering at a speed without braking and driving forces, pitch movement, and longitudinal tire force and align toque could be ignored in this condition. So the established vehicle 3DOF model consists of lateral motion, yaw, and roll motions.Figure 1
Schematic diagram of body movement.
(a)
Plane movement
(b)
Roll movementCorresponding equations of the vehicle model are described as follows:(1)m(v˙y+vxω)-mshsψ¨s=Fyf+Fyr,Izω˙=aFyf-bFyr,Ixsψ¨s+mshsv˙y+vxω-hsψ¨s=msghsψs-Mψ,(2)Mψ=Mψf+Mψr,Mψf=Kfψs+Dfψ˙s,Mψr=Krψs+Drψ˙s.


## 2.2. Nonlinear Tire Model

The lateral tire force depending on the side slip angle is represented as the Pacejka magic formula [9]:(3)Fy=DsinCarctanBα-EBα-arctanBα,where α is the slip angle of tire and B, C, D, and E are, respectively, stiffness factor, shape factor, peak factor, and curvature factor. These factors are obtained by fitting. According to the structure of vehicle tire and vertical load, the relevant factor values of tire on normal road (μ=0.7) are listed (Table 1) [10].Table 1
Curve factors’ values of nonlinear tire.
Tire position (number)
B
C
D
E
Front wheel (2)
10.30
1.30
−11821.63
−1.05
Rear wheel (4)
10.51
−13437.49


## 2.3. Vehicle Dynamic Model Validation

Vehicle test on normal road is to verify the vehicle dynamic model. The details of tests environment and related equipment are showed (Figure2).Figure 2
Details of vehicle tests and related equipment.As shown in Figure2, the integrated test system mainly involves HCZ-1 vehicle handling force and angle measuring instrument, inertial sensors IMU02, vehicle motion parameters measuring instrument by GPS called RLVB20SL, data collection device RLVB-ADC03, and computer. Function of the system is to measure the steering wheel angle input and vehicle motion parameters, such as longitudinal, lateral, and vertical velocity, acceleration, angle, and angular velocity.The test data is processed by VBOX Tools and Matlab. The position deviations between the actual locations of the test equipment and the vehicle’s centroid should be taken into account, so processing the sampling data is necessary and a key step [11].(1) Inertial Sensors IMU02 Direction Correction. Coordinate system of IMU02 (Figure 3) is corrected by right-hand rule and detail (Table 2).Table 2
Coordinate system of IMU02 correction.
Parameter
Longitudinal
Lateral
Vertical
Yaw
Roll
Pitch
Direction correction
+
+
−
−
+
−Figure 3
IMU02 coordinate.(2) Vehicle Centroid State Parameters Correction. According to the deviation between the actual locations of the test equipment and the vehicle’s centroid, sampling data of RLVB20SL and IMU02 should be processed as follows.(1) Vehicle Longitudinal Velocity and Acceleration. Consider(4)vx=vxGPS-sz1γ˙INS,ax=axINS-sz2γ¨INS,where vxGPS is longitudinal velocity tested by RLVB20SL, axINS andγINS are longitudinal acceleration and pitch angle tested by IMU02, sz1 is vertical distance between the main antenna of RLVB20SL and vehicle centroid, and sz2 is vertical distance between the location of IMU02 and vehicle centroid.(2) Vehicle Lateral Velocity and Acceleration. Consider(5)vy=vyGPS-sx1ωINS+sz1ψ˙INS,ay=ayINS-sx2ω˙INS+sz2ψ¨INS,where vyGPS is lateral velocity tested by RLVB20SL, ayINS, ωINS, and ΨINS are lateral acceleration, yaw rate, and roll angle tested by IMU02, sx1 is longitudinal distance between the main antenna of RLVB20SL and vehicle centroid, and sx2 is longitudinal distance between the location of IMU02 and vehicle centroid.As shown in Figure2, the tested vehicle is driven along the double change race on normal road [12]. The velocity and steering wheel angle operated by driver are showed in Figure 4(a). Then, under the same condition, vehicle lateral acceleration by vehicle tests and Matlab simulation applying the established model is obtained, as shown in Figure 4(b).Figure 4
Comparing simulation with vehicle test.
(a)
Vehicle’s velocity and steering wheel angle input
(b)
Vehicle’s lateral acceleration outputClearly, vehicle’s lateral acceleration in Figure4 shows the good consistency between simulation and vehicle test, and the error is less than 10%. So the reliability of the established nonlinear dynamic model is definitely concluded and could be used for the vehicle stability analysis described in the next section.


## 3. Vehicle Steering Stability Analysis

Suppose(Xe,δfe) and Xe=(vye,ωe,ψ˙se,ψse) are one of the balance points of vehicle system. Then the nonlinear formula x˙=f(x) could be described as (6)x˙=Ax+Fx,δf,where A is Jacobian matrix located on the point of (Xe,δfe).Ifδf≠0, the Jacobian matrix A is expressed as(7)A=DxfX,δf(vye,ωe,ψ˙se,ψse)=A11A12A13A14A21A22A23A24A31A32A33A34A41A42A43A44.Characteristic equation of Jacobian matrixA is as follows:(8)λ4+c1λ3+c2λ2+c3λ+c4=0.Hurwitz principle shows the characteristic of the equation [13, 14](9)Δ1=c1>0,Δ2=c1c2-c3>0,Δ3=c3(c1c2-c3)>0,Δ4=c3c4(c1c2-c3)>0.According to Hurwitz principle, if thec1~c4 satisfy (9), the characteristic roots are all negative. If the c1~c4 do not satisfy (9) and c4≠0, positive characteristic root will exist, which predict the unbalance of vehicle system. When the vehicle speed is 20 m/s, the front wheel steering angle is in [0,0.042] rad and interval 0.005 rad, Jacobian matrix characteristic roots with the changes of front wheel steering angle shown in Figure 5.Figure 5
Characteristic roots of Jacobian matrix.Note that the arrow indicates the change direction of the Jacobian matrix characteristic roots in Figure5. It is shown that the characteristic roots are all negative when the front wheel steering angle is small, but they gradually approach 0 as the angle increases and then one of the characteristic roots is 0 if the angle reaches a certain value called critical steering angle. The results given by (a) and (b) in Figure 4 represent characteristic root which is a good way to judge the stability and obtain the critical steering angle and state parameters of vehicle system.


## 4. Stability Regions Estimation

### 4.1. Energy Function Establishing

Assume a quadratic energy functionV(x) as follows [15]:(10)Vx=xTPx,V˙x=x˙TPx+xTPx˙,where xT=[vy,ω,ψ˙s,ψs] and “P” is a positive definite symmetric matrix: (11)P=p11p12p13p14p12p22p23p24p13p23p33p34p14p24p34p44.The matrix “P” is obtained from the next equation:(12)ATP+PA=-I,where “I” is unit matrix and “A” is the system’s Jacobian matrix in origin: A=∂f(x)/∂xx=0.Fact is obvious that the Jacobian matrix is a function of vehicle velocity, so, based on the reference equations, the matrices “A” and “P” in different driving velocities and the corresponding energy functions V(x) and V˙x could be obtained. The matrices “A” and “P” when the vehicle’s velocity is 20 m/s are listed:(13)A=3.725-2.007e1-7.956-1.384e2-6.640e-2-7.503-9.950e-2-1.7321.031e14.140e-1-1.048e1-1.823e20010,P=2.466e-1-2.999e-1-1.395e-1-2.0534-2.999e-18.899e-13.815e-11.638-1.395e-13.815e-12.976e-11.5465-2.05341.6381.54653.563e1.


### 4.2. Boundary Conditions of Stability Region

Generally speaking, traditional principle of the boundary conditions of driving stability region depends on energy attenuation called Lyapunov’s second theorem. Actually, energy attenuation is not enough to calculate the stability region accurately because vehicle is a regular system when it is stable and the space estimated by the energy attenuation characteristic is unreasonably expanded. So the critical vehicle state based on the vehicle system stability analysis should be introduced to determine the vehicle stability region.Boundary conditions mainly lie in two aspects: first, according to Lyapunov’s second theorem,V˙x≤0 is the sufficient condition for system asymptotically stable [16]. Second, steady critical steering angle of front wheel and the corresponding related vehicle state parameters could be obtained by Hurwitz principle, as listed in Table 3.Table 3
Critical values ofV(x).
u(m/s)
δ
c (rad)
Stable equilibrium points of state variables
V
C
β (rad)
ω (rad/s)
ψ
˙
s (rad/s)
ψ
s (rad)
20
0.042
−0.067
0.240
0
0.066
0.15Consequently, the boundary conditions of the vehicle driving stability regions could be described as follows:(14)V˙(x)=0,Vx=Vc.


### 4.3. Vehicle Driving Stability Regions Estimation

According to the proposed boundary conditions of the vehicle stability regions, Figure7 vividly describes the critical surfaces in terms of the boundary conditions.As shown in Figure6, for V˙x=0 and any value of vector x, V˙x is a complexity function of vector x due to the Pacejka tire magic formula, so irregular symmetrical ridge surfaces are formed. Besides, the V(x)=Vc, as quadratic functions of vector x, are described by symmetrical ellipsoids.Figure 6
Critical surfaces in terms of the boundary conditions.Figure 7
Stability region defined byV(x)=0.15 in 20 m/s.Obviously, the space region determined byV˙x=0 driving stability regions is obtained by the intersection of V˙x=0 and V(x)=Vc. So the stability regions are, respectively, showed in Figure 7.


## 4.1. Energy Function Establishing

Assume a quadratic energy functionV(x) as follows [15]:(10)Vx=xTPx,V˙x=x˙TPx+xTPx˙,where xT=[vy,ω,ψ˙s,ψs] and “P” is a positive definite symmetric matrix: (11)P=p11p12p13p14p12p22p23p24p13p23p33p34p14p24p34p44.The matrix “P” is obtained from the next equation:(12)ATP+PA=-I,where “I” is unit matrix and “A” is the system’s Jacobian matrix in origin: A=∂f(x)/∂xx=0.Fact is obvious that the Jacobian matrix is a function of vehicle velocity, so, based on the reference equations, the matrices “A” and “P” in different driving velocities and the corresponding energy functions V(x) and V˙x could be obtained. The matrices “A” and “P” when the vehicle’s velocity is 20 m/s are listed:(13)A=3.725-2.007e1-7.956-1.384e2-6.640e-2-7.503-9.950e-2-1.7321.031e14.140e-1-1.048e1-1.823e20010,P=2.466e-1-2.999e-1-1.395e-1-2.0534-2.999e-18.899e-13.815e-11.638-1.395e-13.815e-12.976e-11.5465-2.05341.6381.54653.563e1.


## 4.2. Boundary Conditions of Stability Region

Generally speaking, traditional principle of the boundary conditions of driving stability region depends on energy attenuation called Lyapunov’s second theorem. Actually, energy attenuation is not enough to calculate the stability region accurately because vehicle is a regular system when it is stable and the space estimated by the energy attenuation characteristic is unreasonably expanded. So the critical vehicle state based on the vehicle system stability analysis should be introduced to determine the vehicle stability region.Boundary conditions mainly lie in two aspects: first, according to Lyapunov’s second theorem,V˙x≤0 is the sufficient condition for system asymptotically stable [16]. Second, steady critical steering angle of front wheel and the corresponding related vehicle state parameters could be obtained by Hurwitz principle, as listed in Table 3.Table 3
Critical values ofV(x).
u(m/s)
δ
c (rad)
Stable equilibrium points of state variables
V
C
β (rad)
ω (rad/s)
ψ
˙
s (rad/s)
ψ
s (rad)
20
0.042
−0.067
0.240
0
0.066
0.15Consequently, the boundary conditions of the vehicle driving stability regions could be described as follows:(14)V˙(x)=0,Vx=Vc.


## 4.3. Vehicle Driving Stability Regions Estimation

According to the proposed boundary conditions of the vehicle stability regions, Figure7 vividly describes the critical surfaces in terms of the boundary conditions.As shown in Figure6, for V˙x=0 and any value of vector x, V˙x is a complexity function of vector x due to the Pacejka tire magic formula, so irregular symmetrical ridge surfaces are formed. Besides, the V(x)=Vc, as quadratic functions of vector x, are described by symmetrical ellipsoids.Figure 6
Critical surfaces in terms of the boundary conditions.Figure 7
Stability region defined byV(x)=0.15 in 20 m/s.Obviously, the space region determined byV˙x=0 driving stability regions is obtained by the intersection of V˙x=0 and V(x)=Vc. So the stability regions are, respectively, showed in Figure 7.


## 5. Stability Region Validation

Appling the integrated test system, vehicle snake shape dynamic behavior is tested. Details of the vehicle’s velocity and steering wheel angles with time are showed in Figure8.Figure 8
Test velocity and steering angle.Figure9 shows the corresponding driving stability region and the obtained distribution of slip angle, yaw rate, and roll angle. Besides, for sake of well comparing the test results with the calculated driving stability region, details of yaw rate and roll angle phase plane are illustrated in Figure 11.Figure 9
Comparing the stability region to vehicle tests results.As shown in Figure9, it is clear that when the vehicle is stable, the distributions of the relevant state variables of vehicle are both in the stability region. However, as the steering wheel angle increases and vehicle system has approached its critical state, trend from stability to instability is obvious. When the steering wheel angle reaches the critical value, the distribution of the state variables while skidding is outside the stability region. Details are showed in Figure 10.Figure 10
Details of vehicle’s dynamic parameters’ distribution.Figure 11
Vehicle’s lateral acceleration with time.Figure11 represents the lateral acceleration of the vehicle and the maximum value of the acceleration is up to 0.5 g which has beyond the critical state. Actually, there is slip phenomenon during the test.The contrast results indicate that the stability region of the vehicle’s stability has good consistency with the test results. Therefore, the proposed driving stability region is reliable for judging and evaluating the stability of vehicle.


## 6. Conclusion

(1)
3DOF vehicle dynamic model with nonlinear tire formulas is deduced and simulated by Matlab. Applying advanced GPS and inertial measurement technology, detecting equipment and data process for vehicle handing stability is integrated and developed, and vehicle in double lane change conditions is tested and the established vehicle model is verified. The reliability of the established nonlinear dynamic model is clearly concluded and could be used for the vehicle stability analysis described in the next section.
(2)
Law of vehicle system’s Jacobian matrix characteristic roots changing with increased front wheel steering angles are fundamentally revealed based on Hurwitz principle. Consequences of the characteristics introduce a good method for defining the boundary of stability region.
(3)
The driving stability region of vehicle verified by vehicle snake shape test provides an effective reference for stability evaluation. The consequences of the vehicle space stability region introduce the new insight into the vehicle performance of stability.



---

*Source: 805063-2015-03-24.xml*