# Design of an Anti-Windup PID Algorithm for Differential Torque Steering Systems

Authors: Dengzhi Peng; Bin Huang; Hao Huang
Journal: Shock and Vibration
 (2022)

---

## Abstract

The EV (electric vehicle) with a wheel hub motor has the advantage of independent driving torque distribution for each wheel, which allows the vehicle performance to be improved. Therefore, a lot of work has been done to investigate the torque allocation algorithm for the mechanical and differential torque integration steering system. To investigate the differential steering process, the 2 DOF (degree of freedom) dual-track reference models with linear and nonlinear tire models are established, and based on the steering process analysis and yaw rate gain calculation, a BP-NN (backpropagation neural network) model is initiated to maintain the accuracy of the calculated yaw rate gain. The limitation of DTSS (differential torque steering system) and the difference of reference models with linear and nonlinear tires are drawn. In addition, an anti-windup variable PID (proportion integration differentiation) controller is designed for the torque distribution. Based on the built 8 DOF model, the vehicle performance indicators are calculated and compared, the gap between the models with linear and nonlinear tires is non-negligible, and a reference model with a nonlinear tire model is recommended for the relevant research. The anti-windup variable PID controller has a better performance than the normal PID controller except for the stability phase plane indicator.


---

## Body

## 1. Introduction

The EV (electric vehicle) has attracted the researcher’s attention for some time due to its advantage in energy conservation and environment protection. Among the different kinds of EVs, the EV with wheel-motor shows superior performance in advanced chassis control algorithms because of its simplified chassis and independent driving torque.The independent driving torque on the vehicle wheels allows the vehicle to have a better ability on torque vectoring (TV) control. The TV control can improve the vehicle’s dynamic performance according to the driver’s attention without significantly damaging the vehicle’s speed. In addition, it provides another steering method, differential torque steering system (DTSS) [1].Most of the relevant research takes the DTSS as a supplementary method of the Ackerman or articulated steering system to improve the vehicle performance [2, 3]. Based on the coordination of active front steering (AFS), direct yaw control (DYC), and TV control, an integrated control algorithm is provided to maintain the vehicle stability during extreme work conditions [4]. In the work of Zheng Hongyu et al., the TV control is combined with the AFS to enhance the yaw and roll stability of a coach [5]. Normally, the torque control distribution algorithm is made up of 2 or 3 layers [6]. The number of layers depends on the preference of the researchers, but the function of the distribution algorithms is similar, including the stability judgment, slip rate calculation, torque allocation, and torque distribution. The phase plane method can represent the nonlinear characteristics; therefore, it is widely used in vehicle stability judgment. The phase plane method can be classified according to the state variable used in the phase plane, including sideslip angle-sideslip angular velocity [7], yaw rate-sideslip angle [8], and front-rear tire sideslip angle [9]. Different torque control algorithms, such as PID (proportion integration differentiation) [10], LQR (linear quadratic regulator control) [11], MPC (model predictive control) [12], and sliding-mode control [13], are initiated to allocate the needed torque during the steering process. For the torque distribution process, even distribution, generalized inverse matrix, and least square are the commonly used algorithms [14]. In addition, the torque distribution inside the motor also attracts the researchers’ attention. Kada Hartani et al. analyzed the impact factor of stator flux, stator resistance, and temperature on the electromagnetic torque of permanent magnet synchronous motor (PMSM) [15].According to the current research, the TV control is combined with the Ackerman steering system, articulated steering system to reduce the steering radius and improve the stability during the steering process. Rarely work has been done to look into the process and limitations of the DTSS itself. In this paper, the steering process of the DTSS is analyzed based on a 2 DOF (degree of freedom) dual-track model, and an anti-windup PID method is used to determine the motor torque. The distribution algorithm is verified by comparing the vehicle performance indicators based on the established 8 DOF vehicle dynamic model.


## 2. Vehicle Dynamic Model

To investigate the differential torque steering control algorithms, the vehicle dynamic model is established based on the vehicle dynamic theory. For an operating vehicle, lots of freedom exists in the whole system, and the simplification must be done based on the researchers’ focus. Normally, a 2 DOF model is established as the reference model and the 8 DOF model is established to verify the torque distribution algorithms in the relevant research [16, 17]. Since the 8 DOF model is much more complicated than the 2 DOF model, the 8 DOF model is introduced first. The 8 DOF model is made up of the vehicle body, wheel, motor, and tire model.


### 2.1. Vehicle Body Model

The 8 DOF model has the motion of longitudinal, lateral, yaw, roll, and rotation of the four wheels, as shown in Figure1. The four motions can be identified as the following equations:(1)mv˙x−vyγ+mshsγϕ˙=∑i=14Fxi−12CdAρvx2−mgf,(2)mv˙y+vxγ−mshsϕ¨=∑i=14Fyi,(3)Ixϕ¨−mshsv˙y+γvx−Ixzγ˙=msghsϕ−Kf+Krϕ−Cf+Crϕ˙,(4)Izγ˙−Ixzϕ¨=M+Fy1+Fy2Lf−Fy3+Fy4Lr,(5)M=Lw2Fx2+Fx4−Fx1−Fx3,where m means the vehicle mass; ms is the sprung mass; vx and vy present the longitudinal and lateral speed; γ is the yaw rate; ϕ is the roll angle; Fxi and Fyi are the longitudinal and lateral force of the four wheels, i=1,2,3,4; Cd presents the air drag coefficient; A is the front area of the vehicle; ρ means the air density; f is the rolling resistance coefficient; Ix and Iz present the roll and yaw inertial of the vehicle; Ixz means the product of inertia of x, z, and axes; hs is the height of the sprung mass; Lf, Lr, and Lw mean the distance between the vehicle gravity center to the front axle, rear axle, and the wheel distance of axle; Kf and Kr are the rolling stiffness of the front and rear axle; and Cf and Cr are the rolling damper of the front and rear axle.Figure 1
Diagram of the 8 DOF model.


### 2.2. Wheel Model

The rotation of the four wheels is taken into consideration in the wheel dynamic model. Each wheel can be demonstrated as follows:(6)Iwω˙i=Ti−FxiRw,where Iw is inertial of the wheel, Ti means the driving torque on the wheel, and Rw is the rolling radius of the wheel. The slip rate of the wheel, si, can be acquired by (7)si=ωiRw−vwiωiRw,where vwi is the speed of the wheel. For the pure DTSS, there is no steering angle for the front and rear axle. For the wheels on the left side, the wheel speeds, vw1 and vw3, can be expressed as(8)vw1,3=vx−Lw2γ.For the wheels on the right side, the wheel speeds,vw2 and vw4, can be expressed as(9)vw2,4=vx+Lw2γ.


### 2.3. Motor Model

For the electrical vehicle, as the driving unit, the motor is quite critical for the vehicle dynamic model. There are kinds of methods to establish the motor model [18]. In this case, the torque distribution algorithm is the concern of our work. Therefore, the model is built with a simple method that can represent the torque characteristic of the motor.The motor model calculates the maximum torque in the current rotation speed, based on the dynamic response characteristic simulated by the first-order inertial response unit, and the output torque of the motor,T, can be acquired.(10)T=11+τtmaxTmmax,Td,where τt is the constant time in a first-order system, Tmmax means the maximum output of the motor in the current rotation speed, and Td represents the demand torque. In this case, the τt is 0.002 s. The external characteristic curve of the motor is shown in Figure 2.Figure 2
Motor external characteristic curve.


### 2.4. Tire Model

The tire system contacts the vehicle to the road surface, and all of the vibration and forces are translated by it. The accuracy of the tire model has a great influence on vehicle performance. A lot of work has been done to establish the tire model, the Fila tire model [19], UA tire model [20], Gim tire model [21], Dugoff tire model [22], HRSI tire model [23], UniTire model [24], and Magic Formula tire model (MF tire) [25] are widely used in the vehicle dynamic model. In this case, two tire models are used: linear tire and MF tire. The linear tire can be expressed as follows:(11)Fy1,2=−kfαf,Fy3,4=−krαr,(12)αf=vy+Lfγvx,αr=vy−Lrγvx,where kf and kr are the cornering stiffness of the front and rear tire and αf and αr are the sideslip angle of the front and rear tires, respectively.For the MF tire model, its parameters are fit from the analysis of the test. It can be expressed as follows:(13)YX=DsinCarctanBX−EBX−arctanBX.WhenY is the longitudinal force, Fx, the variable, X, is the slip rate of the wheel, s. When Y is the lateral force, Fy, the corresponding variable, X, is the sideslip angle the wheel, α. The B, C, D, and E are parameters fitted based on the test on different road types, wheel load, camber angle, temperature, inflation, and tread wear.Based on the vehicle body, wheel, motor, and tire model, the 8 DOF model is established. The 2 DOF model can be built based on the simplification of equations (2) and (4) and the tire model. Equations (2) and (4) can be simplified as(14)mv˙y+vxγ=∑I=14Fyi,Izγ˙=M+Fy1+Fy2Lf−Fy3+Fy4Lr.


## 2.1. Vehicle Body Model

The 8 DOF model has the motion of longitudinal, lateral, yaw, roll, and rotation of the four wheels, as shown in Figure1. The four motions can be identified as the following equations:(1)mv˙x−vyγ+mshsγϕ˙=∑i=14Fxi−12CdAρvx2−mgf,(2)mv˙y+vxγ−mshsϕ¨=∑i=14Fyi,(3)Ixϕ¨−mshsv˙y+γvx−Ixzγ˙=msghsϕ−Kf+Krϕ−Cf+Crϕ˙,(4)Izγ˙−Ixzϕ¨=M+Fy1+Fy2Lf−Fy3+Fy4Lr,(5)M=Lw2Fx2+Fx4−Fx1−Fx3,where m means the vehicle mass; ms is the sprung mass; vx and vy present the longitudinal and lateral speed; γ is the yaw rate; ϕ is the roll angle; Fxi and Fyi are the longitudinal and lateral force of the four wheels, i=1,2,3,4; Cd presents the air drag coefficient; A is the front area of the vehicle; ρ means the air density; f is the rolling resistance coefficient; Ix and Iz present the roll and yaw inertial of the vehicle; Ixz means the product of inertia of x, z, and axes; hs is the height of the sprung mass; Lf, Lr, and Lw mean the distance between the vehicle gravity center to the front axle, rear axle, and the wheel distance of axle; Kf and Kr are the rolling stiffness of the front and rear axle; and Cf and Cr are the rolling damper of the front and rear axle.Figure 1
Diagram of the 8 DOF model.


## 2.2. Wheel Model

The rotation of the four wheels is taken into consideration in the wheel dynamic model. Each wheel can be demonstrated as follows:(6)Iwω˙i=Ti−FxiRw,where Iw is inertial of the wheel, Ti means the driving torque on the wheel, and Rw is the rolling radius of the wheel. The slip rate of the wheel, si, can be acquired by (7)si=ωiRw−vwiωiRw,where vwi is the speed of the wheel. For the pure DTSS, there is no steering angle for the front and rear axle. For the wheels on the left side, the wheel speeds, vw1 and vw3, can be expressed as(8)vw1,3=vx−Lw2γ.For the wheels on the right side, the wheel speeds,vw2 and vw4, can be expressed as(9)vw2,4=vx+Lw2γ.


## 2.3. Motor Model

For the electrical vehicle, as the driving unit, the motor is quite critical for the vehicle dynamic model. There are kinds of methods to establish the motor model [18]. In this case, the torque distribution algorithm is the concern of our work. Therefore, the model is built with a simple method that can represent the torque characteristic of the motor.The motor model calculates the maximum torque in the current rotation speed, based on the dynamic response characteristic simulated by the first-order inertial response unit, and the output torque of the motor,T, can be acquired.(10)T=11+τtmaxTmmax,Td,where τt is the constant time in a first-order system, Tmmax means the maximum output of the motor in the current rotation speed, and Td represents the demand torque. In this case, the τt is 0.002 s. The external characteristic curve of the motor is shown in Figure 2.Figure 2
Motor external characteristic curve.


## 2.4. Tire Model

The tire system contacts the vehicle to the road surface, and all of the vibration and forces are translated by it. The accuracy of the tire model has a great influence on vehicle performance. A lot of work has been done to establish the tire model, the Fila tire model [19], UA tire model [20], Gim tire model [21], Dugoff tire model [22], HRSI tire model [23], UniTire model [24], and Magic Formula tire model (MF tire) [25] are widely used in the vehicle dynamic model. In this case, two tire models are used: linear tire and MF tire. The linear tire can be expressed as follows:(11)Fy1,2=−kfαf,Fy3,4=−krαr,(12)αf=vy+Lfγvx,αr=vy−Lrγvx,where kf and kr are the cornering stiffness of the front and rear tire and αf and αr are the sideslip angle of the front and rear tires, respectively.For the MF tire model, its parameters are fit from the analysis of the test. It can be expressed as follows:(13)YX=DsinCarctanBX−EBX−arctanBX.WhenY is the longitudinal force, Fx, the variable, X, is the slip rate of the wheel, s. When Y is the lateral force, Fy, the corresponding variable, X, is the sideslip angle the wheel, α. The B, C, D, and E are parameters fitted based on the test on different road types, wheel load, camber angle, temperature, inflation, and tread wear.Based on the vehicle body, wheel, motor, and tire model, the 8 DOF model is established. The 2 DOF model can be built based on the simplification of equations (2) and (4) and the tire model. Equations (2) and (4) can be simplified as(14)mv˙y+vxγ=∑I=14Fyi,Izγ˙=M+Fy1+Fy2Lf−Fy3+Fy4Lr.


## 3. Differential Torque Steering Process Analysis

To simplify the analysis, the following assumptions are made in this case. There is no side slip during the differential steering process. The sideslip angle of the tires is the same. The linear tire is used in this section.


### 3.1. Differential Torque Steering Process

Based on the above-mentioned assumption, the state function of the differential torque system can be expressed as follows:(15)x˙=Ax+Bu,(16)x=γvyvx,A=a11a12a21a22=2vxIzLf2kf+Lr2kr2Lfkf−LrkrIz2Lfkf−Lrkrmvx2−12kf+krmvx,B=1Iz0T,u=M.To access the stability of this system, equation (17) is established as follows:(17)detsI−A=s−a11−a12−a21s−a22=s2−a11+a22s+a11a22−a12a21,(18)Δ=a11+a222−4a11a22−a12a21,where s can be gained as follows:(19)s=a11+a222±a11+a222−4a11a22−a12a212,Δ≥0,a11+a222±a11+a222−4a11a22−a12a212i,Δ<0.The condition for a steady system is that all of the eigenvalues of the state matrix,A, have a negative real part. Since the values of kf and kr in a11 and a22 are negative, a11 and a22 always have a negative value. Therefore, if a11a22>a12a21, this system is stable.(20)a11a22−a12a21=4Lf+Lr2kfkrIzmvx2+2Lfkf−LrkrIz>0,in which kfkr >0; therefore, 4Lf+Lr2kfkr/Izmvx2 >0. Thus, when Lfkf−Lrkr >0, the system is always stable. If it is lower than zero, the system will be unstable when the velocity equals the critical speed, vxcr. Otherwise, the system is stable.(21)vxcr=2Lfkf−Lrkrm.


#### 3.1.1. Yaw Rate Gain during the Steady Steering Process

During the steady steering process, the gradient of yaw rate and lateral speed is zero. The equations (2) and (3) can be expressed as(22)2Lfkfαf−2Lrkrαr+M=0,2kfαf+2krαr=mvxγ.Based on equations (11) and (12), and (15), the yaw rate gain during the steady steering process can be acquired.(23)γMs=−kf+krvxLfkf−Lrkrmvx2+2kfkrLf+Lr2.


### 3.2. Torque Distribution Algorithm

According to equation (23), the torque needed to steer in the differential torque steering can be expressed as follows:(24)M=−Lfkf−Lrkrmvx2+2kfkrLf+Lr2kf+krvxγ,where Lf, Lr, and m are constant values. The relevant parameters of the vehicle are shown in Table 1 [26].Table 1
Vehicle parameters.
ParameterValueUnitVehicle mass1704.7kgYawing moment inertia3048.1kg·m2Front tire cornering stiffness−49412N/radRear tire cornering stiffness−47642N/radDistance between the front axle and CG1.035mDistance between the rear axle and CG1.675mRotation moment inertia of the wheel2.5kg·m2Wheel radius0.313mBased on the parameters in Table1, the impact factor of vehicle speed, tire cornering stiffness on differential torque can be analyzed.Equation (24) is established based on the assumption of the linear tire. However, the relationship between the lateral force and tire sideslip angle is nonlinear. When the tire sideslip angle exceeds the range of ±5°, the gap of the lateral force of the tire gets non-negligible. The comparison of the lateral force acquired from the linear tire model and MF tire model of the front and the rear tire is shown in Figures 3 and 4.Figure 3
Comparison of lateral force of the front tire.Figure 4
Comparison of lateral force of the rear tire.According to equation (24) and the variation of sideslip angle and tire cornering stiffness during the steering process, the needed torque can be calculated. In Figure 5, the dash lines represent the torque based on the linear tire, and the solid lines mean the tire cornering stiffness is nonlinear during the process.Figure 5
Comparison of need torque during the steering process.According to Figure5, the needed torque increases with the yaw rate and decreases with the vehicle velocity. The needed torque of taking the tire cornering stiffness as constant is higher than taking it as the variable. When the yaw rate is lower than 0.2 rad/s, the velocity is higher than 4 m/s, and the calculated torque in the two methods is similar.Based on the torque characteristic of the motor, the output torque decreases with the rotation speed. For the vehicle equipped with a wheel-hub motor, the output also decreases with the vehicle speed. According to equation (1), the driving torque needs to overcome the resistance force, including the wind, slope, and rolling resistance force. In this case, the relationship of the force in the longitudinal direction is shown in Figure 6.Figure 6
Variation of the force in the longitudinal direction.According to the research of Antanaitis and Brent [27], there is a different maximum speed for the vehicle on slopes with a different angle. The vehicle has the greatest resistance force when driving on a slope of 4 degrees; therefore, they are set as the boundary condition in Figure 5. This vehicle will have to decrease the speed during the steering process when the speed reaches 35 m/s. The maximum driving force drops sharply at the velocity of 20 m/s, the force that can be used for differential steering also drops dramatically. Thus, the steering process at the speed of 20 m/s is chosen as the work condition for further analysis.Different yaw rate and tire models are set in the 8 DOF model to investigate the impact factor of the tire model, and the result is shown in Figure4. The gap between the result of the linear tire and MF tire gets smaller with the increase of the vehicle speed; therefore, only the minimum gap condition, velocity is 20 m/s, is compared.As shown in Figure7, the yaw rate can reach the ideal value in all the conditions. The response of the yaw rate with the linear tire is worse than the MF tire, and the gap between them gets larger with the increase of the wanted yaw rate. In conclusion, the gap gets larger when the vehicle speed is lower. Therefore, to maintain the accuracy of the torque acquired by equation (24), the cornering stiffness of the front and rear tire, kf and kr, should be adjusted with the vehicle’s working condition. The sideslip angle of the tire is hard to measure when the vehicle is operating. Thus, a BP-NN (backpropagation neural network) model is established to acquire the tire cornering stiffness based on the vehicle velocity and yaw rate.Figure 7
Comparison of yaw rate with different tire models.


### 3.3. BP-NN Model

The BP-NN is the improved feedforward algorithm of ANN (artificial neural network), which was developed by Rumelhart et al. [28]. BP-NN is widely used in parameter fitting and optimization in many fields, such as manufacturing, aerospace, and vehicle design [29–32]. Figure 8 illustrates a typical configuration of the BP-NN.Figure 8
Diagram of BP-NN.In Figure8, there is one input layer, one output layer, and one or more hidden layers between them. In the input layer, there are n inputs; for the output layer, m outputs are included. In addition, s neurons lie in the hidden layer. Suppose the output of the hidden layer and output layer are bj and yk, the threshold value of the hidden layer and out layer is θj and θk, the transfer functions of the hidden layer and output layer are f1 and f2, and the weight factor between the input layer and hidden layer is wij, while the factor between the hidden layer and output layer is wjk. For the hidden layer, the output of jth neuro can be gained by(25)bj=f1∑i=1nwijxi−θj.The result for the output layer is as follows:(26)yk=f2∑j=1swjkbj−θk.In equations (25) and (26), i=1,2,…,n;j=1,2,…,s;k=1,2,…,m.In this case, there are two inputs, the yaw rate and vehicle speed, and two outputs, the cornering stiffness of the front and rear tire. The samples of the inputs and outputs used for the training process are gained by a 2 DOF reference model. The training function of the established BP-NN is set as “Levenberg–Marquardt,” and the adaption learning function is set as “Gradient descent with momentum weight and bias learning function.” The detailed setting of the layers is shown in Table2.Table 2
Parameter of the BP-NN model.
Properties for layerLayer 1Layer 2Layer 3Number of neurons201515Transfer functionsHyperbolic tangent sigmoid transfer functionThe accuracy of the established BP-NN model should be verified before further analysis. Normally, the coefficient of determination,R2; root mean square error, RMSE; and maximum absolute percentage error, MAPE, are used to assess the accuracy of the fitted model [33, 34]. These indicators can be gained according to the equations as follows:(27)R2=∑i=1ntesty⌢i−y¯2∑i=1ntestyi−y¯2,RMSE=1num∑i=1numy⌢i−yi2,(28)MAPE=maxy⌢i−yiyi×100%,where y⌢i and yi represent the calculated value of the fitted model and real value of the samples, num means the number of samples, and y¯ is the mean value of the test result. According to the above equations and the original and fitted values, the R2, RMSE, and MAPE of the two tire cornering stiffness fitted models can be acquired and listed in Table 3.Table 3
Error analysis of the approximation model.
Fitted modelR2RMSEMAPE (%)Front tire0.980.911.01Rear tire0.991.152.12ForR2, the value closer to 1 means better accuracy. For RMSE and MAPE, the value closer to 0 means better accuracy. The error analysis in Table 3 shows that the established BP-NN model meets the requirement for further analysis.


### 3.4. Anti-Windup Variable PID Controller

To keep the output torque of the motor from saturation, an anti-windup variable PID controller is established based on the traditional PID controller. The output torque of the motors,Ti, can be expressed as follows:(29)Ti=Ti′±MReff2Lw,Ti∈Tmax,Tmin,where the footprint, i=1,2,3,4, means the four motors of the vehicle, Ti′ is the drive torque in the normal work condition, ΔTi represents the torque needed to initiate the differential torque steering process, and Tmax and Tmin are the maximum and minimum torque of the motor in the current state.For a typical negative feedback PID control system, the controller can be expressed as(30)u=Kpe+KiKp∫edt+KdKpdedt,where the control parameters Kp, Ki, and Kd are tuned according to the object model or the output curve. u is the output of the controller, and e means the error between the input and output. For the vehicle system, it is a time-varying system with nonlinear characteristics, and a PID controller with constant parameters is not effective enough [35]. To deal with the situation, the Gaussian PID method is proposed, in which the Gaussian function is used to define the control parameter of the PID controller. The Gaussian function can be given by [36](31)fδ=k1−k1−k0exp−qδ2,where δ is the input signal; k0 and k1 are the bonding limit for the zero and infinity input, respectively; and the parameter q defines the concavity openness of the output curve, which can be gained based on the reference input δr and the other variable in equation (28).(32)q=−Infδ−k1/k0−k1δr2=−In1−λδr2,where −δr and δr are the startup and ending point of the upward and downward transient in the upward Gaussian function (as Figure 9 shows); at these two points, fδ reaches to the λ percent of the k0-to-k1 range.Figure 9
Gaussian function. (a) Downward. (b) Upward.
(a)(b)To determine the Gaussian function type of the PID controller parameters, the yaw rate with differentKp, Ki, and Kd is gained and compared in Figures 10, 11, and 12.Figure 10
Yaw rate variation with differentKp.Figure 11
Yaw rate variation with differentKi.Figure 12
Yaw rate variation with differentKd.In Figure10, the yaw rate response oscillation amplitude gets smaller with the increase of Kp, the fluctuation disappears when the value is 27000, and it keeps increasing, which will lead to a higher accommodation time. Therefore, a higher Kp when the response error is relatively big and a lower Kp when the error gets smaller are preferred, and the upward Gaussian function is chosen for Kp.In Figure11, a higher Ki means bigger oscillation amplitude and accommodation time, and the ideal Ki should increase with the response error; thus, the downward Gaussian function is set for Ki.In Figure12, a higher Kd means longer accommodation time and lower oscillation amplitude when the error is bigger, and the ideal Kd should decrease with the response error; thus, the upward Gaussian function is set for Kd. Based on the (30) and analysis of Figures 10, 11, and 12, the three control parameters are determined based on the following equations [37]:(33)Kp=Kp01+Kp11−exp−qe2,Ki=Ki01+Ki1exp−qe2,Kd=Kd01+Kd11−exp−qe2.


## 3.1. Differential Torque Steering Process

Based on the above-mentioned assumption, the state function of the differential torque system can be expressed as follows:(15)x˙=Ax+Bu,(16)x=γvyvx,A=a11a12a21a22=2vxIzLf2kf+Lr2kr2Lfkf−LrkrIz2Lfkf−Lrkrmvx2−12kf+krmvx,B=1Iz0T,u=M.To access the stability of this system, equation (17) is established as follows:(17)detsI−A=s−a11−a12−a21s−a22=s2−a11+a22s+a11a22−a12a21,(18)Δ=a11+a222−4a11a22−a12a21,where s can be gained as follows:(19)s=a11+a222±a11+a222−4a11a22−a12a212,Δ≥0,a11+a222±a11+a222−4a11a22−a12a212i,Δ<0.The condition for a steady system is that all of the eigenvalues of the state matrix,A, have a negative real part. Since the values of kf and kr in a11 and a22 are negative, a11 and a22 always have a negative value. Therefore, if a11a22>a12a21, this system is stable.(20)a11a22−a12a21=4Lf+Lr2kfkrIzmvx2+2Lfkf−LrkrIz>0,in which kfkr >0; therefore, 4Lf+Lr2kfkr/Izmvx2 >0. Thus, when Lfkf−Lrkr >0, the system is always stable. If it is lower than zero, the system will be unstable when the velocity equals the critical speed, vxcr. Otherwise, the system is stable.(21)vxcr=2Lfkf−Lrkrm.


### 3.1.1. Yaw Rate Gain during the Steady Steering Process

During the steady steering process, the gradient of yaw rate and lateral speed is zero. The equations (2) and (3) can be expressed as(22)2Lfkfαf−2Lrkrαr+M=0,2kfαf+2krαr=mvxγ.Based on equations (11) and (12), and (15), the yaw rate gain during the steady steering process can be acquired.(23)γMs=−kf+krvxLfkf−Lrkrmvx2+2kfkrLf+Lr2.


## 3.1.1. Yaw Rate Gain during the Steady Steering Process

During the steady steering process, the gradient of yaw rate and lateral speed is zero. The equations (2) and (3) can be expressed as(22)2Lfkfαf−2Lrkrαr+M=0,2kfαf+2krαr=mvxγ.Based on equations (11) and (12), and (15), the yaw rate gain during the steady steering process can be acquired.(23)γMs=−kf+krvxLfkf−Lrkrmvx2+2kfkrLf+Lr2.


## 3.2. Torque Distribution Algorithm

According to equation (23), the torque needed to steer in the differential torque steering can be expressed as follows:(24)M=−Lfkf−Lrkrmvx2+2kfkrLf+Lr2kf+krvxγ,where Lf, Lr, and m are constant values. The relevant parameters of the vehicle are shown in Table 1 [26].Table 1
Vehicle parameters.
ParameterValueUnitVehicle mass1704.7kgYawing moment inertia3048.1kg·m2Front tire cornering stiffness−49412N/radRear tire cornering stiffness−47642N/radDistance between the front axle and CG1.035mDistance between the rear axle and CG1.675mRotation moment inertia of the wheel2.5kg·m2Wheel radius0.313mBased on the parameters in Table1, the impact factor of vehicle speed, tire cornering stiffness on differential torque can be analyzed.Equation (24) is established based on the assumption of the linear tire. However, the relationship between the lateral force and tire sideslip angle is nonlinear. When the tire sideslip angle exceeds the range of ±5°, the gap of the lateral force of the tire gets non-negligible. The comparison of the lateral force acquired from the linear tire model and MF tire model of the front and the rear tire is shown in Figures 3 and 4.Figure 3
Comparison of lateral force of the front tire.Figure 4
Comparison of lateral force of the rear tire.According to equation (24) and the variation of sideslip angle and tire cornering stiffness during the steering process, the needed torque can be calculated. In Figure 5, the dash lines represent the torque based on the linear tire, and the solid lines mean the tire cornering stiffness is nonlinear during the process.Figure 5
Comparison of need torque during the steering process.According to Figure5, the needed torque increases with the yaw rate and decreases with the vehicle velocity. The needed torque of taking the tire cornering stiffness as constant is higher than taking it as the variable. When the yaw rate is lower than 0.2 rad/s, the velocity is higher than 4 m/s, and the calculated torque in the two methods is similar.Based on the torque characteristic of the motor, the output torque decreases with the rotation speed. For the vehicle equipped with a wheel-hub motor, the output also decreases with the vehicle speed. According to equation (1), the driving torque needs to overcome the resistance force, including the wind, slope, and rolling resistance force. In this case, the relationship of the force in the longitudinal direction is shown in Figure 6.Figure 6
Variation of the force in the longitudinal direction.According to the research of Antanaitis and Brent [27], there is a different maximum speed for the vehicle on slopes with a different angle. The vehicle has the greatest resistance force when driving on a slope of 4 degrees; therefore, they are set as the boundary condition in Figure 5. This vehicle will have to decrease the speed during the steering process when the speed reaches 35 m/s. The maximum driving force drops sharply at the velocity of 20 m/s, the force that can be used for differential steering also drops dramatically. Thus, the steering process at the speed of 20 m/s is chosen as the work condition for further analysis.Different yaw rate and tire models are set in the 8 DOF model to investigate the impact factor of the tire model, and the result is shown in Figure4. The gap between the result of the linear tire and MF tire gets smaller with the increase of the vehicle speed; therefore, only the minimum gap condition, velocity is 20 m/s, is compared.As shown in Figure7, the yaw rate can reach the ideal value in all the conditions. The response of the yaw rate with the linear tire is worse than the MF tire, and the gap between them gets larger with the increase of the wanted yaw rate. In conclusion, the gap gets larger when the vehicle speed is lower. Therefore, to maintain the accuracy of the torque acquired by equation (24), the cornering stiffness of the front and rear tire, kf and kr, should be adjusted with the vehicle’s working condition. The sideslip angle of the tire is hard to measure when the vehicle is operating. Thus, a BP-NN (backpropagation neural network) model is established to acquire the tire cornering stiffness based on the vehicle velocity and yaw rate.Figure 7
Comparison of yaw rate with different tire models.


## 3.3. BP-NN Model

The BP-NN is the improved feedforward algorithm of ANN (artificial neural network), which was developed by Rumelhart et al. [28]. BP-NN is widely used in parameter fitting and optimization in many fields, such as manufacturing, aerospace, and vehicle design [29–32]. Figure 8 illustrates a typical configuration of the BP-NN.Figure 8
Diagram of BP-NN.In Figure8, there is one input layer, one output layer, and one or more hidden layers between them. In the input layer, there are n inputs; for the output layer, m outputs are included. In addition, s neurons lie in the hidden layer. Suppose the output of the hidden layer and output layer are bj and yk, the threshold value of the hidden layer and out layer is θj and θk, the transfer functions of the hidden layer and output layer are f1 and f2, and the weight factor between the input layer and hidden layer is wij, while the factor between the hidden layer and output layer is wjk. For the hidden layer, the output of jth neuro can be gained by(25)bj=f1∑i=1nwijxi−θj.The result for the output layer is as follows:(26)yk=f2∑j=1swjkbj−θk.In equations (25) and (26), i=1,2,…,n;j=1,2,…,s;k=1,2,…,m.In this case, there are two inputs, the yaw rate and vehicle speed, and two outputs, the cornering stiffness of the front and rear tire. The samples of the inputs and outputs used for the training process are gained by a 2 DOF reference model. The training function of the established BP-NN is set as “Levenberg–Marquardt,” and the adaption learning function is set as “Gradient descent with momentum weight and bias learning function.” The detailed setting of the layers is shown in Table2.Table 2
Parameter of the BP-NN model.
Properties for layerLayer 1Layer 2Layer 3Number of neurons201515Transfer functionsHyperbolic tangent sigmoid transfer functionThe accuracy of the established BP-NN model should be verified before further analysis. Normally, the coefficient of determination,R2; root mean square error, RMSE; and maximum absolute percentage error, MAPE, are used to assess the accuracy of the fitted model [33, 34]. These indicators can be gained according to the equations as follows:(27)R2=∑i=1ntesty⌢i−y¯2∑i=1ntestyi−y¯2,RMSE=1num∑i=1numy⌢i−yi2,(28)MAPE=maxy⌢i−yiyi×100%,where y⌢i and yi represent the calculated value of the fitted model and real value of the samples, num means the number of samples, and y¯ is the mean value of the test result. According to the above equations and the original and fitted values, the R2, RMSE, and MAPE of the two tire cornering stiffness fitted models can be acquired and listed in Table 3.Table 3
Error analysis of the approximation model.
Fitted modelR2RMSEMAPE (%)Front tire0.980.911.01Rear tire0.991.152.12ForR2, the value closer to 1 means better accuracy. For RMSE and MAPE, the value closer to 0 means better accuracy. The error analysis in Table 3 shows that the established BP-NN model meets the requirement for further analysis.


## 3.4. Anti-Windup Variable PID Controller

To keep the output torque of the motor from saturation, an anti-windup variable PID controller is established based on the traditional PID controller. The output torque of the motors,Ti, can be expressed as follows:(29)Ti=Ti′±MReff2Lw,Ti∈Tmax,Tmin,where the footprint, i=1,2,3,4, means the four motors of the vehicle, Ti′ is the drive torque in the normal work condition, ΔTi represents the torque needed to initiate the differential torque steering process, and Tmax and Tmin are the maximum and minimum torque of the motor in the current state.For a typical negative feedback PID control system, the controller can be expressed as(30)u=Kpe+KiKp∫edt+KdKpdedt,where the control parameters Kp, Ki, and Kd are tuned according to the object model or the output curve. u is the output of the controller, and e means the error between the input and output. For the vehicle system, it is a time-varying system with nonlinear characteristics, and a PID controller with constant parameters is not effective enough [35]. To deal with the situation, the Gaussian PID method is proposed, in which the Gaussian function is used to define the control parameter of the PID controller. The Gaussian function can be given by [36](31)fδ=k1−k1−k0exp−qδ2,where δ is the input signal; k0 and k1 are the bonding limit for the zero and infinity input, respectively; and the parameter q defines the concavity openness of the output curve, which can be gained based on the reference input δr and the other variable in equation (28).(32)q=−Infδ−k1/k0−k1δr2=−In1−λδr2,where −δr and δr are the startup and ending point of the upward and downward transient in the upward Gaussian function (as Figure 9 shows); at these two points, fδ reaches to the λ percent of the k0-to-k1 range.Figure 9
Gaussian function. (a) Downward. (b) Upward.
(a)(b)To determine the Gaussian function type of the PID controller parameters, the yaw rate with differentKp, Ki, and Kd is gained and compared in Figures 10, 11, and 12.Figure 10
Yaw rate variation with differentKp.Figure 11
Yaw rate variation with differentKi.Figure 12
Yaw rate variation with differentKd.In Figure10, the yaw rate response oscillation amplitude gets smaller with the increase of Kp, the fluctuation disappears when the value is 27000, and it keeps increasing, which will lead to a higher accommodation time. Therefore, a higher Kp when the response error is relatively big and a lower Kp when the error gets smaller are preferred, and the upward Gaussian function is chosen for Kp.In Figure11, a higher Ki means bigger oscillation amplitude and accommodation time, and the ideal Ki should increase with the response error; thus, the downward Gaussian function is set for Ki.In Figure12, a higher Kd means longer accommodation time and lower oscillation amplitude when the error is bigger, and the ideal Kd should decrease with the response error; thus, the upward Gaussian function is set for Kd. Based on the (30) and analysis of Figures 10, 11, and 12, the three control parameters are determined based on the following equations [37]:(33)Kp=Kp01+Kp11−exp−qe2,Ki=Ki01+Ki1exp−qe2,Kd=Kd01+Kd11−exp−qe2.


## 4. Simulation and Discussion

The torque distributed by reference model with a linear tire with PID controller, a nonlinear tire with constant PID, and Gaussian PID controller during the steady-state steering process is initiated to the 8 DOF model. The yaw rate (Figure13), roll angle (Figure 14), sideslip angle and sideslip angular velocity (Figure 15), routine (Figure 16) of the vehicle, and tire sideslip angle (Figure 17), and output torque of the motors (Figure 18) are simulated and discussed in this part.Figure 13
Comparison of yaw rate response.Figure 14
Comparison of the vehicle roll angle.Figure 15
Comparison of sideslip angle and side slip angular velocity phase plane.Figure 16
Comparison of vehicle routine.Figure 17
Comparison of tire sideslip angle. (a) Front left tire. (b) Front right tire. (c) Rear left tire. (d) Rear right tire.
(a)(b)(c)(d)Figure 18
Comparison of motor output torque. (a) Front left tire. (b) Front right tire. (c) Rear left tire. (d) Rear right tire.
(a)(b)(c)(d)According to the comparison in Figure13, the linear tire with PID controller has the maximum peak value, longest accommodation time, but minimum delay time. The MF tire model with constant PID and Gaussian PID has a relatively close response, and the Gaussian PID controller has a lower peak value and accommodation time; in conclusion, it has the best response.In Figure14, the linear tire model has the maximum and minimum peak value, the MF tire model with the Gaussian PID algorithm has a higher peak value than the MF tire model with the constant PID algorithm. When the vehicle maintaining in a stable state (0∼2.5 s and after 15 s), the roll angle gap between these models can be ignored. During the variation period (2.5∼15 s), the roll angle of the Gaussian PID algorithm has the lowest fluctuation range.Figure15 shows the sideslip angle and sideslip angular velocity phase plane and the variation of these two indicators of the three models. The area between the purplish red dash lines is the stable area. It is obvious that the track of the linear tire model with PID exceeds the stable boundary, the track of nonlinear tires stays in the stable region, and the MF tire with constant PID has the smallest region; thus, it is more stable.Figure16 illustrates the vehicle routine of these three models. The routine of the two MF tires is close to each other, and the gap between the MF tire and linear tire is obvious.The sideslip angle of the tires is compared in Figure17. The trend of the tires is the same, and the sideslip angle decreases with the increase of yaw rate, starts to increase (about 7 s) after the yaw rate reached 0.4 rad/s (5 s), and becomes stable after some fluctuation. The linear tire has the minimum valley value, while the maximum peak value is calculated by the MF tire with a constant PID algorithm. The sideslip angle of the MF tire model with Gaussian PID control has the minimum fluctuation range.According to Figure18, the output torque of the left tires and right tires has the same trend. For the left tires, two valley values occur around 4 and 9 seconds, and a peak value emerges around 7 s. The torque of the linear tire model has the maximum peak value, while the MF tire with Gaussian PID has the minimum valley value. The torque variation range of the two models with MF tires is similar, which is lower than the linear tire model. Compared to the left tire, the output torque of the right tires has the opposite trend before the torque reaches a stable value, the torque increases with the yaw rate at 2.5 s, reaches a peak value around 5 s, and then decreases to 7 s and climb to a peak value at 9 s; at last, it decreases to a stable value. The linear tire has the minimum valley value, the first peak value of the three models is close to each other, and the linear tire model and MF tire model with constant PID have a higher 2nd peak value than the MF tire model with Gaussian PID algorithm.


## 5. Conclusion

The EV equipped with a wheel-hub motor allows each wheel to have different driving torque, which is a benefit to the vehicle’s dynamic control and performance improvement. Most of the current research focuses on the integrated control of DTSS and mechanic steering system. In this case, the DTSS itself is analyzed and an anti-windup PID control algorithm to distribute the motor torque is provided to improve the vehicle performance.Based on the vehicle dynamic theory, the 2 DOF reference models with linear and nonlinear tire models are established, the differential torque steering process is analyzed, and the yaw rate gain during the steady steering process is calculated. The DTSS has to kill the speed or acceleration ability to maintain the steering performance during critical conditions, such as running on the upward slope with a relatively high speed. The difference between the two 2 DOF models is compared. To maintain the accuracy of the torque distribution method, a BP-NN model is established, and an anti-windup PID controller is provided.According to the simulation of the distributed torque by linear tire model with PID, nonlinear tire model with PID, and nonlinear tire model with Gaussian PID, the vehicle performance indicators, such as the yaw rate, roll angle, sideslip angle, sideslip angular velocity, sideslip angle of tires, and motor torques are compared. The gap between the model with linear tire and nonlinear tire is non-negligible, and a reference model with a nonlinear tire model is recommended for the relevant research. Except for the sideslip angle and sideslip angular velocity phase plane track, the anti-windup variable PID controller has a better result than the normal PID control.In this case, only the steady steering scenario is taken as the input, and the vehicle performance in the other typical condition, such as double lane-change, fishhook, and snaking, should be taken into consideration. The PID parameter determination method is based on experience, and some adaptive or optimization algorithms, such as the genetic algorithm (GA) and particle swarm optimization (PSO), can be initiated in the parameter determination of the controller.



---

*Source: 9973379-2022-04-11.xml*