# Research on Curved Path Tracking Control for Four-Wheel Steering Vehicle considering Road Adhesion Coefficient

Authors: Runqiao Liu; Minxiang Wei; Nan Sang; Jianwei Wei
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

Curved path tracking control is one of the most important functions of autonomous vehicles. First, small turning radius circular bends considering bend quadrant and travel direction restrictions are planned by polar coordinate equations. Second, an estimator of a vehicle state parameter and road adhesion coefficient based on an extended Kalman filter is designed. To improve the convenience and accuracy of the estimator, the combined slip theory, trigonometric function group fitting, and cubic spline interpolation are used to estimate the longitudinal and lateral forces of the tire model (215/55 R17). Third, to minimize the lateral displacement and yaw angle tracking errors of a four-wheel steering (4WS) vehicle, the front-wheel steering angle of the 4WS vehicle is corrected by a model predictive control (MPC) feed-back controller. Finally, CarSim® simulation results show that the 4WS autonomous vehicle based on the MPC feed-back controller can not only significantly improve the curved path tracking performance but also effectively reduce the probability of drifting or rushing out of the runway at high speeds and on low-adhesion roads.


---

## Body

## 1. Introduction

Path tracking control enables autonomous vehicles to travel accurately, smoothly, and safely along a predetermined path and to remain stable under all driving conditions [1]. As early as 1980, MacAdam [2] proposed an optimal preview path tracking control method, which could minimize the tracking error in channel-changing simulation. Currently, this control method has been applied to the path tracking model in CarSim® software [3]. In addition, the model predictive control (MPC) method can easily handle the nonholonomic constraints of a vehicle system and solve the trade-off problem between path tracking and driving safety [4]. Due to these advantages, many researchers have studied the MPC path tracking method. Borrelli et al. [5] presented a nonlinear MPC approach for active front-wheel steering control. Katriniok and Abel [6] used a predictive controller based on a linear time-varying model (LTV-MPC) to minimize tracking error and response lag. Kim et al. [7] used an MPC path tracking control method based on quadratic programming (QP) optimization to improve the path tracking performance of autonomous vehicles. Schwarting et al. [8] used a receding horizon planner based on nonlinear model predictive control (NMPC) to ensure safety collision avoidance. However, the above control methods ignore the influence of high speeds and low-adhesion roads on curved path tracking.Road adhesion coefficient plays an important role in vehicle electronic control systems. Accurate road adhesion coefficient estimation is critical to improve the path tracking performance of autonomous vehicles. In recent years, to estimate accurate road adhesion coefficient, scholars have proposed various estimation methods. Hahn et al. [9] reported a road friction coefficient estimation algorithm based on a differential global positioning system (DGPS). Wenzel et al. [10] proposed an estimation method based on a double extended Kalman filter (DEKF). Enisz et al. [11] used an extended Kalman filter (EKF) method to estimate the instantaneous and maximum values of road friction coefficient. Zhao et al. [12] designed an EKF method based on a braking dynamics model considering the load transfer of front and rear axles.For steering on a wet and slippery road, when the front-wheel loses most of the lateral force (grip), the vehicle (front-drive) will push ahead or rush out of the runway (understeer); when the rear wheel loses most of the lateral force (grip), the vehicle (rear-drive) will drift out or have a tail flicking (oversteer) [13]. To reduce the understeer and oversteer in path tracking, Hang and Chen [14] used a path tracking integrated control based on four-wheel steering (4WS) and direct yaw-moment control to improve the path tracking performance of vehicles on low-adhesion roads. It is generally known that the handling stability and driving safety of 4WS vehicles are much higher than those of front-wheel steering (FWS) vehicles [15, 16]. As early as 1986, Sano et al. [17] proposed a speed function and a steering angle function for 4WS vehicles. In a small steering wheel angle range or at high speed, the rear wheels of the 4WS vehicle were steered in the same direction as the front ones. In a large steering wheel angle range or at low speed, the rear wheels of the 4WS vehicle were steered in the opposite direction to the front ones. This greatly improves the handling stability and driving safety of 4WS vehicles. In [18, 19], by setting the zero-sideslip angle of 4WS vehicle, the stability at high speed and the sensitivity at low speed were significantly improved. In a recent study by Liu et al. [20], a feed-forward controller that can adjust the steering sensitivity coefficient of the 4WS vehicle based on zero-sideslip angle according to the steering angular frequency and vehicle speed is proposed, but this controller is not suitable for low-adhesion curved roads.The rest of this paper is organized as follows. First of all, small turning radius circular bends (single circle center) considering different restrictions on bend quadrant and travel direction are planned by polar coordinate equations in Section2. Then, in Section 3, an estimator based on the extended Kalman filter for the vehicle state parameter and road adhesion coefficient is designed. To improve the convenience and accuracy of the estimator, we use the combined slip theory, trigonometric function group fitting, and cubic spline interpolation to estimate the longitudinal and lateral forces of the tire model (215/55 R17). Next, in Section 4, a 4WS vehicle based on a MPC feed-back controller is used to improve the curved path tracking performance at high speeds and on low-adhesion roads, and the lateral displacement and yaw angle tracking errors of the 4WS vehicle in path tracking are minimized. Furthermore, CarSim® simulation effects for the optimal preview feed-back and MPC feed-back controllers under different curved paths, vehicle speeds, and road adhesion coefficients are compared in Section 5. Finally, Section 6 presents the conclusions of the work.


## 2. Circular Bend and Vehicle Model

### 2.1. Circular Bend

It is necessary for autonomous vehicles to plan a predetermined curved path in curved path tracking [21–25]. The predetermined path (target path) is usually obtained by high-precision vehicle positioning system, combining the global positioning system (GPS) with the inertial navigation system (INS) [26]. In the absence of electronic map (E-map) data, to facilitate calculations and simulations on small turning radius circular bends (single circle center), the circular bends used as target paths can be planned by polar coordinate equations.A polar coordinate equation can express a circular bend simply and accurately, and the arc lengthL0 of the circular bend is defined as follows:(1)L0=θ0R,where θ0 is the maximum radian of the circular bend and R is the turning radius (m).To plan all forms of circular bends based on the ground inertial coordinate system, we introduce constraints on bend quadrant and travel direction. The equations of the radian with respect tot for the circular bends under different restrictions are described as follows.(1)
Travel clockwise:(2)θi=π2−uctR,i=1,θi=−uctR,i=2,θi=−π2−uctR,i=3,θi=π−uctR,i=4.(2)
Travel anticlockwise:(3)θi=uctR,i=1,θi=−π2+uctR,i=2,θi=π+uctR,i=3,θi=π2+uctR,i=4,where t is the driving time; θi is the radian value corresponding to time t; uc represents the longitudinal speed based on vehicle coordinate system; i=1 is the first quadrant; i=2 is the second quadrant; i=3 is the third quadrant; and i=4 is the fourth quadrant. Then, the coordinate equations for the circular bends are expressed as follows:(4)Xd=X0+Rcosθi,Yd=Y0+Rsinθi,where X0,Y0 is the circle center coordinate; Yd is the desired lateral displacement based on the ground inertial coordinate system; and Xd is the desired longitudinal displacement based on the ground inertial coordinate system.After the derivative of the function ofYd with respect to Xd, the function of the desired yaw angle Ωd with respect to Xd is expressed as follows:(5)ΩdXd=arctanYd′Xd,where Ωd represents the desired yaw angle based on the ground inertial coordinate system.After the derivative of the function ofΩd with respect to t, the function of the desired yaw rate rd with respect to t is formulated as follows:(6)rdt=Ωd′t.Taking the restriction of traveling clockwise as an example, Figure1 shows the right-angle circular bends under four different quadrants, where α is the corner angle and 0≤α≤180 (deg).Figure 1
Right-angle circular bends under four different quadrants: (a) fourth quadrant; (b) first quadrant; (c) third quadrant; (d) second quadrant.Taking the restrictions of traveling anticlockwise and the second quadrant (i=2) as an example, Figure 2 shows the four circular bend diagrams with acute, right, obtuse, and U-shaped angles that meet the restrictions.Figure 2
Four circular bend diagrams under the restrictions of travel anticlockwise and second quadrant (i=2): (a) acute angle; (b) right angle; (c) obtuse angle; (d) U-shaped angle.
(a)
(b)
(c)
(d)Finally, the target paths can be converted into the desired front and rear wheel steering angles and used as initial inputs of autonomous vehicles. For a 4WS vehicle, the yaw rate response is generally set to a first-order lag response and the sideslip angle response is generally set to 0 [27–31]. Its initial inputs can be formulated as follows:(7)δfdsδsws=Izs+mbuc+laCaf/uclCaf⋅Grd1+Trs,δrdsδsws=−Izs−mauc+lbCar/uclCar⋅Grd1+Trs,where s is a Laplace operator; δsw=rd/Grd; Grd is the desired steering sensitivity coefficient and Grd=uc/i0l+i0lKuc2; K is the stability factor and K=mbCar−aCaf/l2CafCar; δfd and δrd are initial inputs of front and rear wheel steering angles for 4WS vehicles, respectively; and Tr=1/l2CafCar1+Kuc2/uc2mIz.


### 2.2. Vehicle Model

A vehicle model (D-class Sedan) of CarSim® is used in the following simulations. The vehicle model is front-drive (Generic 200 KW powertrain with an automatic transmission). A tire model (215/55 R17) used in the vehicle model is a standard lookup table tire model. The main parameters for the vehicle model are shown in Table1.Table 1
Parameters of D-class sedan.
Model parameter
Notation
Unit
Value
Yaw moment of inertia
I
z
kg
⋅
m
2
4607.47
Vehicle mass
m
kg
1530
Distances from vehicle C.G. to front axle
a
m
1.11
Distances from vehicle C.G. to rear axle
b
m
1.66622
Front-wheel track width
T
wf
m
1.55
Rear wheel track width
T
wr
m
1.55
Vehicle wheelbase
l
m
2.77622
Vehicle height
h
0
m
1.471
Vehicle length
l
0
m
4.52
Vehicle width
w
0
m
1.8
Cornering stiffness of front axle
C
af
N
/
rad
97937∗ 2
Cornering stiffness of rear axle
C
ar
N
/
rad
70287∗ 2
Longitudinal stiffness of front axle
C
lf
N
/
rad
85740∗ 2
Longitudinal stiffness of rear axle
C
lr
N
/
rad
57110∗ 2
Initial transmission ratio ofδsw to δf
i
0
23
Front-wheel optimal slip ratio
s
f
0.2
Rear wheel optimal slip ratio
s
r
0.2


## 2.1. Circular Bend

It is necessary for autonomous vehicles to plan a predetermined curved path in curved path tracking [21–25]. The predetermined path (target path) is usually obtained by high-precision vehicle positioning system, combining the global positioning system (GPS) with the inertial navigation system (INS) [26]. In the absence of electronic map (E-map) data, to facilitate calculations and simulations on small turning radius circular bends (single circle center), the circular bends used as target paths can be planned by polar coordinate equations.A polar coordinate equation can express a circular bend simply and accurately, and the arc lengthL0 of the circular bend is defined as follows:(1)L0=θ0R,where θ0 is the maximum radian of the circular bend and R is the turning radius (m).To plan all forms of circular bends based on the ground inertial coordinate system, we introduce constraints on bend quadrant and travel direction. The equations of the radian with respect tot for the circular bends under different restrictions are described as follows.(1)
Travel clockwise:(2)θi=π2−uctR,i=1,θi=−uctR,i=2,θi=−π2−uctR,i=3,θi=π−uctR,i=4.(2)
Travel anticlockwise:(3)θi=uctR,i=1,θi=−π2+uctR,i=2,θi=π+uctR,i=3,θi=π2+uctR,i=4,where t is the driving time; θi is the radian value corresponding to time t; uc represents the longitudinal speed based on vehicle coordinate system; i=1 is the first quadrant; i=2 is the second quadrant; i=3 is the third quadrant; and i=4 is the fourth quadrant. Then, the coordinate equations for the circular bends are expressed as follows:(4)Xd=X0+Rcosθi,Yd=Y0+Rsinθi,where X0,Y0 is the circle center coordinate; Yd is the desired lateral displacement based on the ground inertial coordinate system; and Xd is the desired longitudinal displacement based on the ground inertial coordinate system.After the derivative of the function ofYd with respect to Xd, the function of the desired yaw angle Ωd with respect to Xd is expressed as follows:(5)ΩdXd=arctanYd′Xd,where Ωd represents the desired yaw angle based on the ground inertial coordinate system.After the derivative of the function ofΩd with respect to t, the function of the desired yaw rate rd with respect to t is formulated as follows:(6)rdt=Ωd′t.Taking the restriction of traveling clockwise as an example, Figure1 shows the right-angle circular bends under four different quadrants, where α is the corner angle and 0≤α≤180 (deg).Figure 1
Right-angle circular bends under four different quadrants: (a) fourth quadrant; (b) first quadrant; (c) third quadrant; (d) second quadrant.Taking the restrictions of traveling anticlockwise and the second quadrant (i=2) as an example, Figure 2 shows the four circular bend diagrams with acute, right, obtuse, and U-shaped angles that meet the restrictions.Figure 2
Four circular bend diagrams under the restrictions of travel anticlockwise and second quadrant (i=2): (a) acute angle; (b) right angle; (c) obtuse angle; (d) U-shaped angle.
(a)
(b)
(c)
(d)Finally, the target paths can be converted into the desired front and rear wheel steering angles and used as initial inputs of autonomous vehicles. For a 4WS vehicle, the yaw rate response is generally set to a first-order lag response and the sideslip angle response is generally set to 0 [27–31]. Its initial inputs can be formulated as follows:(7)δfdsδsws=Izs+mbuc+laCaf/uclCaf⋅Grd1+Trs,δrdsδsws=−Izs−mauc+lbCar/uclCar⋅Grd1+Trs,where s is a Laplace operator; δsw=rd/Grd; Grd is the desired steering sensitivity coefficient and Grd=uc/i0l+i0lKuc2; K is the stability factor and K=mbCar−aCaf/l2CafCar; δfd and δrd are initial inputs of front and rear wheel steering angles for 4WS vehicles, respectively; and Tr=1/l2CafCar1+Kuc2/uc2mIz.


## 2.2. Vehicle Model

A vehicle model (D-class Sedan) of CarSim® is used in the following simulations. The vehicle model is front-drive (Generic 200 KW powertrain with an automatic transmission). A tire model (215/55 R17) used in the vehicle model is a standard lookup table tire model. The main parameters for the vehicle model are shown in Table1.Table 1
Parameters of D-class sedan.
Model parameter
Notation
Unit
Value
Yaw moment of inertia
I
z
kg
⋅
m
2
4607.47
Vehicle mass
m
kg
1530
Distances from vehicle C.G. to front axle
a
m
1.11
Distances from vehicle C.G. to rear axle
b
m
1.66622
Front-wheel track width
T
wf
m
1.55
Rear wheel track width
T
wr
m
1.55
Vehicle wheelbase
l
m
2.77622
Vehicle height
h
0
m
1.471
Vehicle length
l
0
m
4.52
Vehicle width
w
0
m
1.8
Cornering stiffness of front axle
C
af
N
/
rad
97937∗ 2
Cornering stiffness of rear axle
C
ar
N
/
rad
70287∗ 2
Longitudinal stiffness of front axle
C
lf
N
/
rad
85740∗ 2
Longitudinal stiffness of rear axle
C
lr
N
/
rad
57110∗ 2
Initial transmission ratio ofδsw to δf
i
0
23
Front-wheel optimal slip ratio
s
f
0.2
Rear wheel optimal slip ratio
s
r
0.2


## 3. Estimator Based on EKF

### 3.1. Vehicle State Parameter Estimation

Ignoring the effects of the suspension system, road roughness, and air resistance, a three-degree-of-freedom (3DOF) vehicle dynamics model considering lateral, longitudinal, and yaw motions is used to estimate the longitudinal speed, sideslip angle, and yaw rate of vehicles [32]. The state space equation of the vehicle dynamics model is shown as follows:(8)r˙=−a2Caf+b2CarIzucr−aCaf−bCarIzβ+aCafIzδf,β˙=bCar−aCafIz−1r−Caf+Carmucβ+Cafmucδf,u˙c=rβuc+ax.The measurement equation of the vehicle model is shown as follows:(9)ay=bCar−aCafmucr−Caf+Carmβ+Cafmδf,where δf is the front-wheel steering angle, r is the yaw rate based on the vehicle coordinate system, β is the sideslip angle based on the vehicle coordinate system, ax is the longitudinal acceleration based on the vehicle coordinate system, ay is the lateral acceleration based on the vehicle coordinate system, the input vector is u=δfaxT, the state vector is x=rβucT, and the measurement vector is y=ay. After the partial differential of equations (8) and (9) with respect to x, the Jacobian matrices F and H in Figure 3 are calculated as follows:(10)F=−Carb2−Cafa2Izuc−Cafa+bCarIza2Caf+b2CarIzuc2rbCar−aCafmuc2−1−Caf−Carmuc2aCaf−bCarrmuc3+Caf+Carβ−Cafδfmuc2βucrucrβ,H=bCar−aCafmuc−Caf−CarmaCaf−bCarmuc2r.Figure 3
Flow chart of EKF algorithm.Then, the flow chart of EKF algorithm is shown in Figure3, where Φ=eF∗Δt≈I+FΔt; Δt=0.01 is the sampling time; Q1 is the covariance matrix of process noise; R1 is the covariance matrix of measurement noise; f is the state space equation; x⌢k is the state estimation vector at time k; Pk is the state error covariance matrix at time k; x⌢k−1 is the state estimation vector at time k−1; and Pk−1 is the state error covariance matrix at time k−1. Given the initial values x⌢k−1 and Pk−1 and x⌢k and Pk can be updated cyclically. The above estimation process is continually looped until the number of set steps.


### 3.2. Tire Longitudinal and Lateral Forces Estimation

Professor Pacejka put forward a “magic formula” tire model in 1987, which estimates the tire longitudinal and lateral forces by constructing the trigonometric function group. The tire longitudinal and lateral force estimation equations under a single working condition are shown as follows:(11)Fx0=μ0DxsinCxarctanBxλ−EBxλ−arctanBxλx,Fy0=μ0DysinCyarctanByβ−EyByβ−arctanByβ,where μ0=1; λ is the slip ratio; Bx and By are the stiffness factors of the tire longitudinal and lateral forces, respectively; Cx and Cy are the curve shape factors of the tire longitudinal and lateral forces, respectively; Dx and Dy are the curve peak factors of the tire longitudinal and lateral forces, respectively; and Ex and Ey are the curvature factors at the curve peak of the tire longitudinal and lateral forces, respectively. These factors are represented by equations containing several different coefficients [33].Although the “magic formula” has a high fitting accuracy for the tire longitudinal and lateral forces, these coefficients in the above factors are difficult to fit [34, 35]. In order to estimate the longitudinal and lateral forces of the tire model (215/55 R17) conveniently and accurately, we use the combined slip theory, trigonometric function group fitting, and cubic spline interpolation to estimate longitudinal and lateral forces of the tire model.Firstly, combined with the longitudinal and lateral force tabular data of the tire model (the tabular data are measured on a laboratory or on-road tester, and they are all measured under the condition that the road adhesion coefficientμ is 1), tabular data under 10 different tire loads are selected, which are represented by the purple, pink, green, cyan, and blue lines in Figure 4. Using the trigonometric function group to fit these tabular data, the fitting results are represented by the five types of black lines in Figure 4. These factors after fitting are shown as follows (RMSE is the root mean square error):(1)
Tire load 1725:By1 is 9.342; Cy1 is 2.753; Dy1 is 1891.4; Ey1 is 1.123; and RMSE is 0.0096(2)
Tire load 3500:By2 is 9.909; Cy2 is 2.694; Dy2 is 3698.8; Ey2 is 1.114; and RMSE is 0.0089(3)
Tire load 6100:By3 is 10.17; Cy3 is 2.626; Dy3 is 5382.1; Ey3 is 1.109; and RMSE is 0.0077(4)
Tire load 6950:By4 is 9.943; Cy4 is 2.573; Dy4 is 6971.6; Ey4 is 1.112; and RMSE is 0.0066(5)
Tire load 9005:By5 is 9.029; Cy5 is 2.565; Dy5 is 8564.5; Ey5 is 1.126; and RMSE is 0.0055(6)
Tire load 2105:Bx1 is 11.37; Cx1 is 1.528; Dx1 is 2729.5; Ex1 is −0.5744; and RMSE is 0.0016(7)
Tire load 3995:Bx2 is 11.37; Cx2 is 1.528; Dx2 is 5459.1; Ex2 is −0.5744; and RMSE is 0.0016.(8)
Tire load 6120:Bx3 is 11.37; Cx3 is 1.528; Dx3 is 8188.6; Ex3 is −0.5744; and RMSE is 0.0016(9)
Tire load 7900:Bx4 is 11.37; Cx4 is 1.528; Dx4 is 10918; Ex4 is −0.5744; and RMSE is 0.0016(10)
Tire load 10100:Bx5 is 11.37; Cx5 is 1.528; Dx5 is 13648; Ex5 is −0.5745; and RMSE is 0.0016Figure 4
Longitudinal and lateral force tabular data underμ=1: (a) the lateral tire forces under pure steering conditions; (b) the longitudinal tire forces under pure driving/braking conditions.
(a)
(b)Secondly, the cubic spline interpolation is used to estimate these factors under other tire loads. Taking theBy estimation under a certain tire load as an example, the circular calling format of the cubic spline interpolation is shown as follows:(12)yz=interp1z0,yi,:,z,“spline”,where z0=1725,3500,6100,6950,9005; z is a certain tire load; yi,:=Byi,i=1:5; yz is the interpolation result; and “spline” is the function about cubic spline interpolation. Similarly, the remaining factors can be estimated. Then, the longitudinal force estimation equation (under pure braking/driving conditions) and lateral force estimation equation (under pure steering conditions) under any tire load can be obtained.Finally, using the combined slip theory [36, 37], the tire longitudinal and lateral forces under joint conditions can be described as follows:(13)Fxij=μ0Fx0Fz,λ⋅σxσ,Fyij=μ0Fy0Fz,β⋅σyσ,where μ0=1; σx=−λ/1+λ; σy=−tanβ/1+λ; σ=σx2+σy2; Fz is the tire load; i=f,r; j=l,r; Fxfl, Fxrl, Fxfr, and Fxrr are the left front, left rear, right front, and right rear wheel longitudinal forces under joint conditions, respectively; and Fyfl, Fyrl, Fyfr, and Fyrr are the left front, left rear, right front, and right rear wheel lateral forces under joint conditions, respectively.


### 3.3. Road Adhesion Coefficient Estimation

A 3DOF vehicle dynamics model involving tire lateral and longitudinal forces under joint conditions is used to estimate the road adhesion coefficient [38]. The state space equation of the vehicle dynamics model is shown as follows:(14)ay=uflFxflδf+Fyfl+ufrFxfrδf+Fyfr+μrlFyrl+μrrFyrrm,ax=μflFxfl−Fyflδf+μfrFxfr−Fyfrδf+μrlFxrl+μrrFxrrm,r˙=uflaFxflδf+Fyfl+Twf/2−Fxfl+Fyflδf+ufraFxfrδf+Fyfr+Twf/2Fxfr−Fyfrδf+urf−bFyrl−Twr/2Fxrl+urr−bFyrr+Twr/2FxrrIz,where the input vector u=δfFxflFxfrFxrlFxrrFyflFyfrFyrlFyrrT; the state vector x=μflμfrμrlμrrT, μfl, μfr, μrl, and μrr are the left front, right front, left rear, and right rear wheel adhesion coefficients, respectively; and the measurement vector y=axayr˙T. After the partial differential of equation (14) with respect to x, the Jacobian matrices F and H in Figure 3 are calculated as follows:(15)F=1000010000100001,H=Fxfl−FyflδfmFxfr−FyfrδfmFxrlmFxrrmFxflδf+FyflFxfrδf+FyfrFyrlFyrrH3,1H3,2H3,3H3,4,H3,1=aFxflδf+Fyfl+Twf−Fxfl+Fyflδf/2Iz,H3,3=−bFyrl−TwrFxrl/2Iz,H3,2=aFxfrδf+Fyfr+TwfFxfr−Fyfrδf/2Iz,H3,4=−bFyrr+TwrFxrr/2Iz.The flow chart of the road adhesion coefficient estimation algorithm is similar to that of the vehicle state parameter estimation, as shown in Figure3.


### 3.4. Verification of the Estimator

The overall block diagram of the estimator for vehicle state parameter and road adhesion coefficient is shown in Figure5.Figure 5
Overall block diagram of the estimator.A joint adhesion road is used in the verification. The features of the joint adhesion road are detailed as follows: a road of a total length of about 321 m, with the first 28 m a high-μ section (a road adhesion coefficient of 0.85), and the rest a low-μ section (a road adhesion coefficient of 0.5). In addition, a two-way single-lane road is used in the estimation, and it is a right-angled circular bend (the corner angle α is 90 deg), as shown in Figure 6(a), where the width of each lane is 6 m, the blue dashed line is the target path, and the turning radius is 37.5 m. The estimated road adhesion coefficient and vehicle state parameters at uc=50km/h are shown in Figures 6(b)–6(e), where a uniformly distributed random signal (the minimum value is 0 and the maximum value is 0.02) is added to the actual outputs of CarSim® vehicle considering uncertain disturbances. The parameter settings of vehicle state parameter estimation are shown as follows: Q1=I3×3; x⌢k−1=0050/3.6; Pk−1=I3×3; R1=0.01; and Pk can be updated cyclically in the form of 3×3 matrix. The parameter settings of road adhesion coefficient estimation are shown as follows: Q1=I4×4; x⌢k−1=0.50.50.50.5; Pk−1=I4×4; R1=I3×3; and Pk can be updated cyclically in the form of 4×4 matrix. Figure 6 shows that even in the case of a small turning radius and low-adhesion road, the estimator can quickly and accurately estimate the vehicle state parameters and road adhesion coefficient.Figure 6
Estimation results: (a) animation effect; (b) road adhesion coefficient; (c) sideslip angle; (d) yaw rate; (e) vehicle speed.
(a)
(b)
(c)
(d)
(e)


## 3.1. Vehicle State Parameter Estimation

Ignoring the effects of the suspension system, road roughness, and air resistance, a three-degree-of-freedom (3DOF) vehicle dynamics model considering lateral, longitudinal, and yaw motions is used to estimate the longitudinal speed, sideslip angle, and yaw rate of vehicles [32]. The state space equation of the vehicle dynamics model is shown as follows:(8)r˙=−a2Caf+b2CarIzucr−aCaf−bCarIzβ+aCafIzδf,β˙=bCar−aCafIz−1r−Caf+Carmucβ+Cafmucδf,u˙c=rβuc+ax.The measurement equation of the vehicle model is shown as follows:(9)ay=bCar−aCafmucr−Caf+Carmβ+Cafmδf,where δf is the front-wheel steering angle, r is the yaw rate based on the vehicle coordinate system, β is the sideslip angle based on the vehicle coordinate system, ax is the longitudinal acceleration based on the vehicle coordinate system, ay is the lateral acceleration based on the vehicle coordinate system, the input vector is u=δfaxT, the state vector is x=rβucT, and the measurement vector is y=ay. After the partial differential of equations (8) and (9) with respect to x, the Jacobian matrices F and H in Figure 3 are calculated as follows:(10)F=−Carb2−Cafa2Izuc−Cafa+bCarIza2Caf+b2CarIzuc2rbCar−aCafmuc2−1−Caf−Carmuc2aCaf−bCarrmuc3+Caf+Carβ−Cafδfmuc2βucrucrβ,H=bCar−aCafmuc−Caf−CarmaCaf−bCarmuc2r.Figure 3
Flow chart of EKF algorithm.Then, the flow chart of EKF algorithm is shown in Figure3, where Φ=eF∗Δt≈I+FΔt; Δt=0.01 is the sampling time; Q1 is the covariance matrix of process noise; R1 is the covariance matrix of measurement noise; f is the state space equation; x⌢k is the state estimation vector at time k; Pk is the state error covariance matrix at time k; x⌢k−1 is the state estimation vector at time k−1; and Pk−1 is the state error covariance matrix at time k−1. Given the initial values x⌢k−1 and Pk−1 and x⌢k and Pk can be updated cyclically. The above estimation process is continually looped until the number of set steps.


## 3.2. Tire Longitudinal and Lateral Forces Estimation

Professor Pacejka put forward a “magic formula” tire model in 1987, which estimates the tire longitudinal and lateral forces by constructing the trigonometric function group. The tire longitudinal and lateral force estimation equations under a single working condition are shown as follows:(11)Fx0=μ0DxsinCxarctanBxλ−EBxλ−arctanBxλx,Fy0=μ0DysinCyarctanByβ−EyByβ−arctanByβ,where μ0=1; λ is the slip ratio; Bx and By are the stiffness factors of the tire longitudinal and lateral forces, respectively; Cx and Cy are the curve shape factors of the tire longitudinal and lateral forces, respectively; Dx and Dy are the curve peak factors of the tire longitudinal and lateral forces, respectively; and Ex and Ey are the curvature factors at the curve peak of the tire longitudinal and lateral forces, respectively. These factors are represented by equations containing several different coefficients [33].Although the “magic formula” has a high fitting accuracy for the tire longitudinal and lateral forces, these coefficients in the above factors are difficult to fit [34, 35]. In order to estimate the longitudinal and lateral forces of the tire model (215/55 R17) conveniently and accurately, we use the combined slip theory, trigonometric function group fitting, and cubic spline interpolation to estimate longitudinal and lateral forces of the tire model.Firstly, combined with the longitudinal and lateral force tabular data of the tire model (the tabular data are measured on a laboratory or on-road tester, and they are all measured under the condition that the road adhesion coefficientμ is 1), tabular data under 10 different tire loads are selected, which are represented by the purple, pink, green, cyan, and blue lines in Figure 4. Using the trigonometric function group to fit these tabular data, the fitting results are represented by the five types of black lines in Figure 4. These factors after fitting are shown as follows (RMSE is the root mean square error):(1)
Tire load 1725:By1 is 9.342; Cy1 is 2.753; Dy1 is 1891.4; Ey1 is 1.123; and RMSE is 0.0096(2)
Tire load 3500:By2 is 9.909; Cy2 is 2.694; Dy2 is 3698.8; Ey2 is 1.114; and RMSE is 0.0089(3)
Tire load 6100:By3 is 10.17; Cy3 is 2.626; Dy3 is 5382.1; Ey3 is 1.109; and RMSE is 0.0077(4)
Tire load 6950:By4 is 9.943; Cy4 is 2.573; Dy4 is 6971.6; Ey4 is 1.112; and RMSE is 0.0066(5)
Tire load 9005:By5 is 9.029; Cy5 is 2.565; Dy5 is 8564.5; Ey5 is 1.126; and RMSE is 0.0055(6)
Tire load 2105:Bx1 is 11.37; Cx1 is 1.528; Dx1 is 2729.5; Ex1 is −0.5744; and RMSE is 0.0016(7)
Tire load 3995:Bx2 is 11.37; Cx2 is 1.528; Dx2 is 5459.1; Ex2 is −0.5744; and RMSE is 0.0016.(8)
Tire load 6120:Bx3 is 11.37; Cx3 is 1.528; Dx3 is 8188.6; Ex3 is −0.5744; and RMSE is 0.0016(9)
Tire load 7900:Bx4 is 11.37; Cx4 is 1.528; Dx4 is 10918; Ex4 is −0.5744; and RMSE is 0.0016(10)
Tire load 10100:Bx5 is 11.37; Cx5 is 1.528; Dx5 is 13648; Ex5 is −0.5745; and RMSE is 0.0016Figure 4
Longitudinal and lateral force tabular data underμ=1: (a) the lateral tire forces under pure steering conditions; (b) the longitudinal tire forces under pure driving/braking conditions.
(a)
(b)Secondly, the cubic spline interpolation is used to estimate these factors under other tire loads. Taking theBy estimation under a certain tire load as an example, the circular calling format of the cubic spline interpolation is shown as follows:(12)yz=interp1z0,yi,:,z,“spline”,where z0=1725,3500,6100,6950,9005; z is a certain tire load; yi,:=Byi,i=1:5; yz is the interpolation result; and “spline” is the function about cubic spline interpolation. Similarly, the remaining factors can be estimated. Then, the longitudinal force estimation equation (under pure braking/driving conditions) and lateral force estimation equation (under pure steering conditions) under any tire load can be obtained.Finally, using the combined slip theory [36, 37], the tire longitudinal and lateral forces under joint conditions can be described as follows:(13)Fxij=μ0Fx0Fz,λ⋅σxσ,Fyij=μ0Fy0Fz,β⋅σyσ,where μ0=1; σx=−λ/1+λ; σy=−tanβ/1+λ; σ=σx2+σy2; Fz is the tire load; i=f,r; j=l,r; Fxfl, Fxrl, Fxfr, and Fxrr are the left front, left rear, right front, and right rear wheel longitudinal forces under joint conditions, respectively; and Fyfl, Fyrl, Fyfr, and Fyrr are the left front, left rear, right front, and right rear wheel lateral forces under joint conditions, respectively.


## 3.3. Road Adhesion Coefficient Estimation

A 3DOF vehicle dynamics model involving tire lateral and longitudinal forces under joint conditions is used to estimate the road adhesion coefficient [38]. The state space equation of the vehicle dynamics model is shown as follows:(14)ay=uflFxflδf+Fyfl+ufrFxfrδf+Fyfr+μrlFyrl+μrrFyrrm,ax=μflFxfl−Fyflδf+μfrFxfr−Fyfrδf+μrlFxrl+μrrFxrrm,r˙=uflaFxflδf+Fyfl+Twf/2−Fxfl+Fyflδf+ufraFxfrδf+Fyfr+Twf/2Fxfr−Fyfrδf+urf−bFyrl−Twr/2Fxrl+urr−bFyrr+Twr/2FxrrIz,where the input vector u=δfFxflFxfrFxrlFxrrFyflFyfrFyrlFyrrT; the state vector x=μflμfrμrlμrrT, μfl, μfr, μrl, and μrr are the left front, right front, left rear, and right rear wheel adhesion coefficients, respectively; and the measurement vector y=axayr˙T. After the partial differential of equation (14) with respect to x, the Jacobian matrices F and H in Figure 3 are calculated as follows:(15)F=1000010000100001,H=Fxfl−FyflδfmFxfr−FyfrδfmFxrlmFxrrmFxflδf+FyflFxfrδf+FyfrFyrlFyrrH3,1H3,2H3,3H3,4,H3,1=aFxflδf+Fyfl+Twf−Fxfl+Fyflδf/2Iz,H3,3=−bFyrl−TwrFxrl/2Iz,H3,2=aFxfrδf+Fyfr+TwfFxfr−Fyfrδf/2Iz,H3,4=−bFyrr+TwrFxrr/2Iz.The flow chart of the road adhesion coefficient estimation algorithm is similar to that of the vehicle state parameter estimation, as shown in Figure3.


## 3.4. Verification of the Estimator

The overall block diagram of the estimator for vehicle state parameter and road adhesion coefficient is shown in Figure5.Figure 5
Overall block diagram of the estimator.A joint adhesion road is used in the verification. The features of the joint adhesion road are detailed as follows: a road of a total length of about 321 m, with the first 28 m a high-μ section (a road adhesion coefficient of 0.85), and the rest a low-μ section (a road adhesion coefficient of 0.5). In addition, a two-way single-lane road is used in the estimation, and it is a right-angled circular bend (the corner angle α is 90 deg), as shown in Figure 6(a), where the width of each lane is 6 m, the blue dashed line is the target path, and the turning radius is 37.5 m. The estimated road adhesion coefficient and vehicle state parameters at uc=50km/h are shown in Figures 6(b)–6(e), where a uniformly distributed random signal (the minimum value is 0 and the maximum value is 0.02) is added to the actual outputs of CarSim® vehicle considering uncertain disturbances. The parameter settings of vehicle state parameter estimation are shown as follows: Q1=I3×3; x⌢k−1=0050/3.6; Pk−1=I3×3; R1=0.01; and Pk can be updated cyclically in the form of 3×3 matrix. The parameter settings of road adhesion coefficient estimation are shown as follows: Q1=I4×4; x⌢k−1=0.50.50.50.5; Pk−1=I4×4; R1=I3×3; and Pk can be updated cyclically in the form of 4×4 matrix. Figure 6 shows that even in the case of a small turning radius and low-adhesion road, the estimator can quickly and accurately estimate the vehicle state parameters and road adhesion coefficient.Figure 6
Estimation results: (a) animation effect; (b) road adhesion coefficient; (c) sideslip angle; (d) yaw rate; (e) vehicle speed.
(a)
(b)
(c)
(d)
(e)


## 4. MPC Feed-Back Controller

MPC generally includes predictive model, rolling optimization, and feed-back correction. It mainly compensates for system errors through feed-back correction and compensates for the effects caused by time variance and interference through rolling optimization [39, 40]. Since the effects of the tire nonlinearity, tire load transfer, road low-adhesion, and crosswind, the actual trajectories often deviate from the desired trajectories in path tracking [20]. To minimize these deviations, a MPC feed-back controller is used to follow the desired lateral displacement and the desired yaw angle and minimize the tracking errors between the desired values and the actual values. First, a nonlinear dynamic model considering the transformation between the vehicle coordinate system and the ground inertial coordinate system is described as follows:(16)v˙y=−Caf+Carmucvy+bCar−aCafmuc−ucr+Cafmδf+Carmδr+Clfsfmδf+Clrsrmδr,r˙=aClfsfIzδf−bClrsrIzδr−Cafa2+Carb2Izucr+bCar−aCafIzucvy+aCafIzδf−bCarIzδr,u˙c=vyr+Clfsfm+Clrsrm+Cafvy+ar/uc−δfmδf+Cvy−br/uc−δrarmδr,Y˙=ucsinΩ+vycosΩ,X˙=uccosΩ−vysinΩ,where δr is the rear-wheel steering angle; vy is the lateral speed based on the vehicle coordinate system; Y˙ is the lateral speed based on the ground inertial coordinate system; X˙ is the longitudinal speed based on the ground inertial coordinate system; and Ω is the yaw angle based on the ground inertial coordinate system. Then, it is transformed into(17)x˜˙=Ax˜+Bu˜,where(18)x˜=x−xd,u˜=u−ud,u=δf,ud=δfd,x=vyucΩrYXΤ,xd=vyducdΩdrdYdXdΤ,A=−Caf+Carmuc∂fvy∂uc0bCar−aCafmuc−uc00r+Cafδf,t−1muc+Carδr,t−1muc∂fuc∂uc0vy+Cafaδf,t−1muc−Carbδr,t−1muc00000100bCar−aCafIzuc∂fr∂uc0−Cafa2+Carb2Izuc00cosΩsinΩuccosΩ−vysinΩ000−sinΩcosΩ−vycosΩ−ucsinΩ000,B=Cafm+Clfsfm−2Cafδf,t−1+Cafvy+ar/ucm0aClfsfIz+aCafIz00,∂fuc∂uc=−Cafδf,t−1vy+armuc2−Carδr,t−1vy−brmuc2,∂fvy∂uc=Caf+Carmuc2vy+aCaf−bCarmuc2r−r,∂fr∂uc=Cafa2+Carb2Izuc2r−bCar−aCafIzuc2vy.To apply equation (17) to the MPC system, it is discretized as follows:(19)x˜k+1=Akx˜k+Bku˜k,Ak=I+TA,Bk=TB,where T is the sampling time. It is further converted into(20)ξk+1=A˜ξk+B˜Δu˜k,ηk=C˜ξk,where the state vector ξk=x˜ku˜k−1; the input vector Δu˜k=u˜k−u˜k−1; the output vector η=XYT; A˜=An×nBn×m0m×nIm×m and B˜=Bn×mIm×m; and n and m are the dimensions of the state vector and output vector, respectively.Next is a rolling optimization process, we can convert the MPC problem into a quadratic programming problem [41]. Its mathematical model is expressed as follows:(21)minΔU,ε12ΔUεTH2ΔUε+f2TΔUεs.t.AconsΔUεT≤bcons,lb≤ΔUεT≤ub,where ΔU is the input increment; H2 is the quadratic matrix of quadratic programming; f2 is the one item vector of quadratic programming; Acons is the coefficient matrix in linear inequality constraint; bcons is the right-end vector in linear inequality constraint; lb=ΔUmin0T; ub=ΔUmaxεT; and ε is the relaxation factor. Its optimization objective function is set to(22)Jξk,u˜k−1,ΔUk=∑i=1Npηk+ik−ηdk+ikQ22+∑i=1Nc−1Δu˜k+ikR22+ρε2,where Q2 and R2 are weighting matrices; ρ is the weight coefficient; Np is the prediction time domain; and Nc is the control time domain.To avoid a sudden change in the input vector, the problem of constrained optimization for each step is shown as follows [42, 43]:(23)minΔUt,εJξk,u˜k−1,ΔUk,s.t.Δu˜mini+k≤Δu˜i+k≤Δu˜maxi+k,i=0,1,…,Nc−1,u˜mini+k≤u˜i+k≤u˜maxi+k,i=0,1,…,Nc−1,η˜mini+k−ε≤η˜i+k≤η˜maxi+k+ε,i=0,1,…,Nc−1.Calling the quadratic programming function (quadprog), which is expressed as follows:(24)ΔU,fval,exitflag=quadprogH2,f2,Acons,bcons,,,lb,ub,ΔU0,options,where ΔU0=zerosNc+1,1 and options=optimset“Algorithm”,“active−set”.After solving each control cycle, a series of input increments in the control time domain are obtained:(25)ΔU=Δu˜kk,Δu˜k+1k,…,Δu˜k+Nc−1kT.Finally, there is a process of feed-back correction. The first element of the control sequence acts on the MPC system:(26)Uek=uk=uk−1+Δu˜kk,where Ue represents the feed-back control (front-wheel steering angle). The above control process is continually looped until path tracking is completed.


## 5. Curved Path Tracking Simulations

### 5.1. Two-Path Tracking Controllers

Two-path tracking controllers for autonomous vehicles in the following simulations are shown as follows:(1)
The optimal preview feed-back controller for an FWS autonomous vehicle is defined as controller 1, which is an FWS path tracking controller. CarSim® has its own optimal preview path tracking control model located in “Steering by the Closed-loop Driver Model.” The red vehicle in the following simulations represents controller 1.(2)
The MPC feed-back controller for a 4WS autonomous vehicle is defined as controller 2, which is a 4WS path tracking controller. The blue vehicle in the following simulations represents controller 2. The block diagram of controller 2 is shown in Figure7, where δf=δfd+Ue and δr=δrd.Figure 7
Block diagram of controller 2.


### 5.2. Circular Bend Simulations at Low Speeds and on High-Adhesion Roads

For the circular bend simulation, the coordinatesXdandYd of the target path can be updated continuously by polar coordinate equations depending on R, uc, α, and X0,Y0 variables.The target path of this simulation is a right-angled circular bend (travel anticlockwise and second quadrant), which is shown in Figure2(b) (the corner angle α is 90 deg). A one-way four-lane road is applied to the simulation, where the width of each lane is 3.5 m, the turning radius R is 37.5 m, the initial coordinates of the autonomous vehicle are 0,0, the ending coordinates of the autonomous vehicle are 300,300, the circle center coordinates are X0,Y0=262.5,37.5, and the road adhesion coefficient of the dry asphalt road is set to 0.85. The road conditions meet the SAE standards J3087_201710 [44]. The path tracking animation effects and steering angles for the red and blue vehicles at uc=50km/h and μ=0.85 are shown in Figure 8. The longitudinal and lateral displacements, lateral accelerations, yaw angles, and yaw rates for the red and blue vehicles at uc=50km/h and μ=0.85 are shown in Figure 9. The detailed constraints of the MPC are shown as follows: the input increment constraints, −1rad/s≤Δu=Δδf≤1rad/s; the input vector constraints, −0.25rad≤u=δf≤0.25rad; the output vector constraints: 0m0m≤η=XY≤300m300m; Q2=8000008000080; Np=30; Nc=10; T=0.01; ε=1000; ρ=1000; and R2=5∗105. The detailed optimal preview parameter settings are shown as follows: the driver preview time is 0.5 s; the driver time lag is 0 s; the low-speed dynamic limit is 10 km/h; the maximum steering wheel angle is 720 deg; and the maximum steering wheel angle rate is 1200 deg/s.Figure 8
Curved path tracking simulation results for red and blue vehicles atuc=50km/h and μ=0.85: (a) animation effect; (b) front-wheel steering angle; (c) rear-wheel steering angle.
(a)
(b)
(c)Figure 9
Simulation results for red and blue vehicles atuc=50km/h and μ=0.85: (a) Y with respect to X; (b) r with respect to t; (c) Ω with respect to t; (d) ay with respect to t.
(a)
(b)
(c)
(d)It can be seen from Figures8(b) and 8(c) that the red vehicle is an FWS vehicle and the blue vehicle is a 4WS vehicle. Even if the two vehicles have different front and rear wheel steering angles, Figures 9(a) and 9(c) show that they can still achieve effective path tracking without pushing ahead (understeer) or tail flicking (oversteer) at low speeds and on high-adhesion roads. It can be seen that there is little difference in the control effect under steady-state conditions in Figures 9(b) and 9(d).


### 5.3. Circular Bend Simulations at Low Speeds and on Low-Adhesion Roads

The target path of this simulation is the same as in Section5.2. The road adhesion coefficient of the wet and slippery road is set to 0.5. The curved path tracking animation effects and steering angles for red and blue vehicles at uc=50km/h and μ=0.5 are shown in Figure 10. The longitudinal and lateral displacements, lateral accelerations, yaw angles, and yaw rates for red and blue vehicles at uc=50km/h and μ=0.5 are shown in Figure 11. The constraints of the MPC and the optimal preview parameter settings are the same as in Section 5.2.Figure 10
Curved path tracking simulation results for red and blue vehicles atuc=50km/h and μ=0.5: (a) animation effect; (b) front-wheel steering angle; (c) rear-wheel steering angle.
(a)
(b)
(c)Figure 11
Simulation results for red and blue vehicles atuc=50km/h and μ=0.5: (a) Y with respect to X; (b) r with respect to t; (c) Ω with respect to t; (d) ay with respect to t.
(a)
(b)
(c)
(d)Whenμ drops to 0.5, Figure 10(a) shows that the red vehicle has obvious sideslip (the maximum deviation is about 4 m). Figure 10(a) shows that the blue vehicle has smaller sideslip (the maximum deviation is about 1 m), and its path tracking performance is significantly better than that of the red vehicle on the low-adhesion road. The reason for this phenomenon is that the rear wheels of the blue vehicle are steered in the opposite direction to the front ones when entering the corner, and its rear wheels are steered in the same direction as the front ones when driving out of the corner, as shown in Figures 10(b) and 10(c). This situation greatly improves the safety and stability of the 4WS vehicle in path tracking on a small turning radius and low-adhesion curved road [45]. Figure 11 shows that the yaw angle, yaw rate, and lateral acceleration values of the blue vehicle are closer to the desired values, and the vehicle maintains good path tracking performance even on low-adhesion roads with small turning radius.


### 5.4. True Variable Curvature Curved Road Simulations at High Speeds and on Low-Adhesion Roads

For the curved path tracking on real roads, the polar coordinate equation (4) are replaced by the coordinates of the electronic map (E-map), and the coordinates Xd,Yd of the target path can be updated continuously by high-precision vehicle positioning system combining GPS with INS.A true variable curvature curved road is shown as follows: the road adhesion coefficient of the wet and slippery road is set to 0.5; the minimum turning radius of the bend ① isRmin≈167m; the minimum turning radius of the bend ② is Rmin≈143m; and the minimum turning radius of the bend ③ is Rmin≈125m. The coordinates Xd,Yd of the true variable curvature curved road are shown in Figure 12(a). The cross slopes of the true variable curvature curved road are shown in Figure 12(b).Figure 12
True variable curvature curved road: (a) coordinatesXd,Yd and (b) cross slopes.
(a)
(b)The constraints of the MPC are shown as follows: the output vector constraints is0m−350m≤η=XY≤900m0m, and other constraints are the same as in Section 5.2. The optimal preview parameter settings are the same as in Section 5.2. The curved path tracking effects for red and blue vehicles at uc=100km/h and μ=0.5 are shown in Figure 13. Figure 13 shows that the red vehicle loses control and rushes out of the runway when passing the bend ③, but the blue vehicle maintains good path tracking performance even at high speeds and on low-adhesion roads.Figure 13
Curved path tracking effects for red and blue vehicles: (a) bend ①; (b) bend ②; (c) bend ③; (d)Y with respect to X.
(a)
(b)
(c)
(d)


## 5.1. Two-Path Tracking Controllers

Two-path tracking controllers for autonomous vehicles in the following simulations are shown as follows:(1)
The optimal preview feed-back controller for an FWS autonomous vehicle is defined as controller 1, which is an FWS path tracking controller. CarSim® has its own optimal preview path tracking control model located in “Steering by the Closed-loop Driver Model.” The red vehicle in the following simulations represents controller 1.(2)
The MPC feed-back controller for a 4WS autonomous vehicle is defined as controller 2, which is a 4WS path tracking controller. The blue vehicle in the following simulations represents controller 2. The block diagram of controller 2 is shown in Figure7, where δf=δfd+Ue and δr=δrd.Figure 7
Block diagram of controller 2.


## 5.2. Circular Bend Simulations at Low Speeds and on High-Adhesion Roads

For the circular bend simulation, the coordinatesXdandYd of the target path can be updated continuously by polar coordinate equations depending on R, uc, α, and X0,Y0 variables.The target path of this simulation is a right-angled circular bend (travel anticlockwise and second quadrant), which is shown in Figure2(b) (the corner angle α is 90 deg). A one-way four-lane road is applied to the simulation, where the width of each lane is 3.5 m, the turning radius R is 37.5 m, the initial coordinates of the autonomous vehicle are 0,0, the ending coordinates of the autonomous vehicle are 300,300, the circle center coordinates are X0,Y0=262.5,37.5, and the road adhesion coefficient of the dry asphalt road is set to 0.85. The road conditions meet the SAE standards J3087_201710 [44]. The path tracking animation effects and steering angles for the red and blue vehicles at uc=50km/h and μ=0.85 are shown in Figure 8. The longitudinal and lateral displacements, lateral accelerations, yaw angles, and yaw rates for the red and blue vehicles at uc=50km/h and μ=0.85 are shown in Figure 9. The detailed constraints of the MPC are shown as follows: the input increment constraints, −1rad/s≤Δu=Δδf≤1rad/s; the input vector constraints, −0.25rad≤u=δf≤0.25rad; the output vector constraints: 0m0m≤η=XY≤300m300m; Q2=8000008000080; Np=30; Nc=10; T=0.01; ε=1000; ρ=1000; and R2=5∗105. The detailed optimal preview parameter settings are shown as follows: the driver preview time is 0.5 s; the driver time lag is 0 s; the low-speed dynamic limit is 10 km/h; the maximum steering wheel angle is 720 deg; and the maximum steering wheel angle rate is 1200 deg/s.Figure 8
Curved path tracking simulation results for red and blue vehicles atuc=50km/h and μ=0.85: (a) animation effect; (b) front-wheel steering angle; (c) rear-wheel steering angle.
(a)
(b)
(c)Figure 9
Simulation results for red and blue vehicles atuc=50km/h and μ=0.85: (a) Y with respect to X; (b) r with respect to t; (c) Ω with respect to t; (d) ay with respect to t.
(a)
(b)
(c)
(d)It can be seen from Figures8(b) and 8(c) that the red vehicle is an FWS vehicle and the blue vehicle is a 4WS vehicle. Even if the two vehicles have different front and rear wheel steering angles, Figures 9(a) and 9(c) show that they can still achieve effective path tracking without pushing ahead (understeer) or tail flicking (oversteer) at low speeds and on high-adhesion roads. It can be seen that there is little difference in the control effect under steady-state conditions in Figures 9(b) and 9(d).


## 5.3. Circular Bend Simulations at Low Speeds and on Low-Adhesion Roads

The target path of this simulation is the same as in Section5.2. The road adhesion coefficient of the wet and slippery road is set to 0.5. The curved path tracking animation effects and steering angles for red and blue vehicles at uc=50km/h and μ=0.5 are shown in Figure 10. The longitudinal and lateral displacements, lateral accelerations, yaw angles, and yaw rates for red and blue vehicles at uc=50km/h and μ=0.5 are shown in Figure 11. The constraints of the MPC and the optimal preview parameter settings are the same as in Section 5.2.Figure 10
Curved path tracking simulation results for red and blue vehicles atuc=50km/h and μ=0.5: (a) animation effect; (b) front-wheel steering angle; (c) rear-wheel steering angle.
(a)
(b)
(c)Figure 11
Simulation results for red and blue vehicles atuc=50km/h and μ=0.5: (a) Y with respect to X; (b) r with respect to t; (c) Ω with respect to t; (d) ay with respect to t.
(a)
(b)
(c)
(d)Whenμ drops to 0.5, Figure 10(a) shows that the red vehicle has obvious sideslip (the maximum deviation is about 4 m). Figure 10(a) shows that the blue vehicle has smaller sideslip (the maximum deviation is about 1 m), and its path tracking performance is significantly better than that of the red vehicle on the low-adhesion road. The reason for this phenomenon is that the rear wheels of the blue vehicle are steered in the opposite direction to the front ones when entering the corner, and its rear wheels are steered in the same direction as the front ones when driving out of the corner, as shown in Figures 10(b) and 10(c). This situation greatly improves the safety and stability of the 4WS vehicle in path tracking on a small turning radius and low-adhesion curved road [45]. Figure 11 shows that the yaw angle, yaw rate, and lateral acceleration values of the blue vehicle are closer to the desired values, and the vehicle maintains good path tracking performance even on low-adhesion roads with small turning radius.


## 5.4. True Variable Curvature Curved Road Simulations at High Speeds and on Low-Adhesion Roads

For the curved path tracking on real roads, the polar coordinate equation (4) are replaced by the coordinates of the electronic map (E-map), and the coordinates Xd,Yd of the target path can be updated continuously by high-precision vehicle positioning system combining GPS with INS.A true variable curvature curved road is shown as follows: the road adhesion coefficient of the wet and slippery road is set to 0.5; the minimum turning radius of the bend ① isRmin≈167m; the minimum turning radius of the bend ② is Rmin≈143m; and the minimum turning radius of the bend ③ is Rmin≈125m. The coordinates Xd,Yd of the true variable curvature curved road are shown in Figure 12(a). The cross slopes of the true variable curvature curved road are shown in Figure 12(b).Figure 12
True variable curvature curved road: (a) coordinatesXd,Yd and (b) cross slopes.
(a)
(b)The constraints of the MPC are shown as follows: the output vector constraints is0m−350m≤η=XY≤900m0m, and other constraints are the same as in Section 5.2. The optimal preview parameter settings are the same as in Section 5.2. The curved path tracking effects for red and blue vehicles at uc=100km/h and μ=0.5 are shown in Figure 13. Figure 13 shows that the red vehicle loses control and rushes out of the runway when passing the bend ③, but the blue vehicle maintains good path tracking performance even at high speeds and on low-adhesion roads.Figure 13
Curved path tracking effects for red and blue vehicles: (a) bend ①; (b) bend ②; (c) bend ③; (d)Y with respect to X.
(a)
(b)
(c)
(d)


## 6. Conclusions

The results of this paper can be summarized as follows.Since the combined slip theory, trigonometric function group fitting, and cubic spline interpolation are used, it is more convenient and accurate in estimating longitudinal and lateral forces and road adhesion coefficients. Compared to the FWS vehicle with an optimal preview path tracking controller, the 4WS vehicle with a MPC path tracking controller can more effectively improve safety and stability at high speeds and on low-adhesion roads. In the future, the path tracking test in the real 4WS autonomous vehicles needs to be implemented and the suspension and braking control systems for 4WS autonomous vehicles in curved path tracking need to be studied further.



---

*Source: 3108589-2020-01-20.xml*