Lateral Stability Analysis of 4WID Electric Vehicle Based on Sliding Mode Control and Optimal Distribution Torque Strategy

Abstract

In this paper, we propose a lateral stability control strategy for four-wheel independent drive (4WID) electric vehicles. The control strategy adopts a hierarchical structure. First, a seven-degree-of-freedom (7DOF) 4WID electric vehicle model is established. Then, the upper controller adopts the integral sliding mode control (ISMC) method to obtain the desired yaw moment by controlling both the yaw rate and the sideslip angle. A new sliding mode reaching law (NSMRL) is designed to reduce chattering and make state variables converge faster, and the superiority of NSMRL is verified by theoretical analysis. The lower controller proposes a new optimal allocation algorithm, which selects the tire utilization rate and the standard deviation coefficient of the tire utilization rate as the objective function. The safety performance of vehicle is improved, and the instability caused by the significant difference in the stability margin between the four wheels under extreme road conditions is avoided. Finally, a simulation is carried out to verify the effectiveness of the proposed control strategy under single-lane-change and J-turn maneuvers.

1. Introduction

Electric vehicles (EVs) have recently become a research hotspot of new energy vehicles. They have significant advantages in reducing pollution, flexible control, and rapid response. Investment by automotive companies has developed various configurations of electric vehicles, among which EVs with 4WID are one of the most prominent configurations [1,2,3,4]. In 4WID electric vehicles, each in-wheel motor can control one wheel independently and adjust the wheel torque or speed through a more flexible control strategy, improving the vehicle’s handling and driving stability. The stability of EVs in an emergency is essential. Therefore, how to ensure the stability of 4WID electric vehicles under complex working conditions and to make full use of the advantages of independent control of four-wheel drive torque has become an important research topic.

Many studies show that direct yaw moment control (DYC) is one of the most effective control strategies to improve handling stability [5,6,7,8,9]. According to the existing literature, lateral stability control generally adopts a hierarchical structure to keep the control variables following the ideal value [10,11]. The DYC method includes two levels of control, where the upper controller provides the desired yaw moment, and the lower controller reasonably distributes the yaw moment to the four wheels.

The DYC system has nonlinear solid, uncertain, and coupling characteristics. The design of the upper controller is the first difficulty. In the early stage, the PID controller has the advantage of easy implementation and is the most extensively used scheme to generate yaw moment [12,13,14]. However, it is challenging to adjust the parameters of the PID controller, and the character is often limited. Some modern linear control methodologies are proposed to obtain a better control performance. The yaw moment that keeps the vehicle stable is calculated by using the linear quadratic regulator (LQR) method, and the vehicle stability is controlled [15,16,17]. However, the method cannot cope with the solid nonlinear characteristics, especially when the electric vehicle is traveling at high speed. The model predictive control (MPC) method has been adopted on 4WID electric vehicles. The disadvantage of the method is that the calculation of the iterative matrix is extensive, and the implementation of multiple variables is complex [18]. With the development of intelligent technology, fuzzy control is widely used to design the upper controller, and the method can calculate the yaw moment and keep the vehicle stable [19]. Although a fuzzy controller has the advantage of solving nonlinear problems, its control rules are complex, and the membership function is obtained according to mature experience, so it cannot be the best choice for vehicle stability control.

However, there are many interferences and uncertainties in the practice DYC system, which may come internally or externally. The above methods are complicated to use in solving these problems, significantly when the EVs run at a high speed. In recent decades, nonlinear variable structure SMC has been regarded as an effective control technique. Therefore, the SMC method has attracted extensive attention, some theories have been studied, and many achievements have been made [20,21]. The design of SMC includes two parts. The first part is the selection of the sliding mode surface, and the second part is the design of the sliding mode reaching law. Once the sliding surface is determined, the stability and dynamic character of the sliding motion are determined. The traditional linear sliding mode surface can meet the design requirements of the control system, and the parameter design is also easy [22,23,24,25]. Many studies for DYC systems have been based on a linear sliding mode surface, but it considers only yaw rate as a stability index [22,23,24]. The vehicle’s center of mass slip angle and yaw rate characterize the vehicle’s stability from different sides. Single control of any state parameter is inadequate. Therefore, many controllers consider both states simultaneously in the design process to improve stability [25,26,27]. However, its limitation is that the state tracking error of the system will not converge to zero within a finite time. Therefore, a linear sliding mode surface is unsuitable for nonlinear systems with high speed and accuracy requirements. The yaw rate and sideslip angle influence the integral SMC (ISMC). The results show that the controller significantly enhances the tracking performance and yaw damping [28,29]. The above SMC control strategies are based on the two-degree-of-freedom (2DOF) or three-degree-of-freedom (3DOF) bicycle model. The calculation of the yaw moment by adding the unknown additional torque has excellent limitations, which is equivalent to simplifying the actual force state when the vehicle is running, ignoring the influence caused by the change in tire force under extreme driving conditions. Due to the error between the 2DOF/3DOF model and the actual vehicle dynamic model, its control effect will be impacted. In addition, the SMC method suppresses external disturbances and parameter changes through switching functions, but the oscillation caused by chattering limits the advantages of sliding mode control [30]. Therefore, researchers have proposed many schemes to reduce chattering, such as reaching law [31,32], higher order sliding mode [33,34], and fractional sliding mode [35]. The concept of reaching law is introduced to reduce the chattering of the sliding mode. Through the design of the reaching law, the movement speed of the state quantity when approaching the sliding mode surface can be increased, and the chattering vibration can be effectively reduced. The isokinetic reaching law is applied in [22]. However, the coefficient of the reaching law is constant, and there is a considerable contradiction between reducing chattering and shortening the arrival time. Presently, the sliding mode controller design generally adopts the exponential reaching law method. Although it has suitable control performance for a specified system, the parameters of the exponential reaching law are typically fixed and have no self-adjustment function. A novel exponential reaching law is proposed to design the controller, in which the system state variables are used to suppress the chattering problem [36]. However, in the above reaching law, the discontinuity gain decreases rapidly due to the change in the sliding surface function, which reduces the robustness of the controller near the sliding surface and increases the approach time.

The purpose of the lower control of the DYC system is to reasonably distribute the yaw moment obtained by the upper control to the four wheels. Different control distribution methods will yield different results of motor torque. These methods mainly include a rule-based allocation algorithm and a torque allocation method based on optimization. The average distribution algorithm [37] and vertical load dynamic distribution algorithm [38] are rule-based allocation methods. These two methods do not take full advantage of the independent controllability of each wheel. Regarding the torque allocation method based on optimization, many optimization objective functions have been proposed to improve vehicle performance. For instance, reference [39] adopts a multi-objective driving/braking torque allocation algorithm to achieve minimum energy consumption while maintaining vehicle stability. Reference [40] proposes a method based on the minimum tire slip criterion to prevent sideslip due to tire force saturation. In this work, we add the steady-state margin standard deviation coefficient to the optimal objective function. The objective is to prevent the gap in the stability margin between the four tires from becoming too large to affect handling in extreme conditions.

Based on the above discussions, we propose a method combining novel SMC and the optimal distribution torque algorithm. An integral sliding mode surface method based on the 7DOF dynamics vehicle model is proposed in the upper sliding mode controller. The desired yaw moment is obtained by controlling both yaw rate and sideslip angle. A new sliding mode reaching law is designed to solve the problems of the traditional reaching law sliding mode’s surface-reaching time and system chattering, and its superiority is demonstrated. The lower controller adopts the utilization rate of four tires and its standard deviation coefficient as the objective function, which can keep the stability margin gap between tires from being too large.

This paper is divided into six parts. Section 1 introduces the motivation of this study. Section 2 presents the dynamic vehicle model. Section 3 is pertinent to designing the upper controller. The ISMC, based on the 7DOF model, controls the yaw rate and sideslip angle simultaneously to obtain the required yaw moment. Section 4 presents the lower control for optimal distribution. Section 5 evaluates the effectiveness of the proposed method under different road conditions, and the simulation results are given. Finally, the conclusion is given in Section 6.

2. Vehicle Dynamic Model and Problem Formulation

In this section, we first propose a 7DOF vehicle model, including the longitudinal motion, lateral motion, yaw motion, and rotational dynamics of four wheels, and a tire model is given. Then, a 2DOF vehicle model is established to calculate the ideal yaw rate and sideslip angle. Finally, the problem formulation is presented.

2.1. 7DOF Vehicle Model

The 7DOF vehicle model is shown in Figure 1. The vehicle coordinate system satisfies the right-sided requirement, and the positive direction of the X-axis is the vehicle’s forward direction. The following assumptions for establishing the 7DOF model are made for the vehicle:

① The vehicle always travels on a horizontal road, the displacement of the vehicle along the Z-axis is not considered, and the roll angle around the X-axis and the pitch angle around the Y-axis are both zero;

② The front tires of the left and right wheels have the same steering wheel angle, and the rear wheel angle is considered zero;

③ The mechanical properties of the four tires are the same; and

④ The air resistance is ignored.

In Figure 1, ${\omega }_{\gamma }$ is the yaw rate; $\beta$ is the sideslip angle; $L_f$ and $L_r$ are the distances from the front axle and rear axle to the center of centroid, respectively; $L$ is the total distance between the front axle and rear axle of the vehicle, $B_f$ and $B_r$ are the front track width and rear track width, respectively; $\delta$ is the front steering angle; $v_x$ and $v_y$ are the longitudinal velocity and lateral velocity at the center of centroid of vehicle, respectively; ${\alpha }_{\mathrm{ij}}$ is the slip angle; and $F_{xij}$ and $F_{yij}$ denote the longitudinal tire force and the lateral tire force, respectively (the first subscript of ${(\cdot )}_{ij}$ denotes front wheel and rear wheel, and the second subscript denotes left wheel and right wheel).

The longitudinal motion is modeled as

$$m\left({\dot{v}}_x-{\omega }_{\gamma }\cdot v_y\right)=\left(F_{xfl}+F_{xfr}\right)\cos \delta -\left(F_{yfl}+F_{yfr}\right)\sin \delta +F_{xrl}+F_{xrr}$$

(1)

where $m$ is the total mass of the car vehicle.

The lateral motion is expressed as

$$m\left({\dot{v}}_y+{\omega }_{\gamma }\cdot v_x\right)=\left(F_{xfl}+F_{xfr}\right)\sin \delta +\left(F_{yfl}+F_{yfr}\right)\cos \delta +F_{yrl}+F_{yrr}$$

(2)

The yaw motion is described as

$$I_z\cdot {\dot{\omega }}_{\gamma }^{}=\left[\left(F_{yfl}+F_{yfr}\right)\cos \delta \right]L_f+\left[\left(F_{yfl}-F_{yfr}\right)\sin \delta \right]\frac{B_f}{2}-\left(F_{yrl}+F_{yrr}\right)L_r+M_z$$

(3)

where $I_z$ is the yaw moment of inertia of vehicle, and $M_z$ is the yaw moment described by

$$M_z=L_f\left(F_{xfl}+F_{xfr}\right)\sin \delta +\frac{B_f}{2}\left(F_{xfr}-F_{xfl}\right)\cos \delta +\frac{B_r}{2}(F_{xrr}-F_{xrl})$$

(4)

The force state of the wheel during rotation is shown in Figure 2.

When the vehicle is running, the wheels are mainly subjected to the driving torque $T_{dij}$ output by the in-wheel motor, the braking torque $T_{bij}$ applied by the braking system, the rolling resistance torque $T_f$, and the friction force $F_{xij}$ and the vertical force $F_{zij}$ provided by the ground.

The wheel rotational dynamics are given as

$$J{\dot{\omega }}_{ij}=T_{dij}-T_{bij}-F_{xij}R-T_f$$

(5)

where $J$ is the moment of inertia of the wheel, ${\omega }_{ij}$ is the wheel angular velocity, and $R$ is the rolling radius of the tire.

Rolling resistance is mainly caused by the friction between the tire and the ground when the wheel is rolling. The following equation can obtain the rolling resistance torque.

$$T_f=f\cdot F_{zij}R$$

(6)

where $f$ is the tire rolling resistance coefficient.

The vertical load of each tire is given as

$$\left\{\begin{array}{l}{F}_{zfl}=\frac{mg{L}_r}{2L}-\frac{m{h}_g{a}_x}{2L}-\frac{m{h}_g{a}_y}{2B_f}\\ F_{zfr}=\frac{mg{L}_r}{2L}-\frac{m{h}_g{a}_x}{2L}+\frac{m{h}_g{a}_y}{2B_f}\\ F_{zrl}=\frac{mg{L}_f}{2L}+\frac{m{h}_g{a}_x}{2L}-\frac{m{h}_g{a}_y}{2B_r}\\ F_{zrr}=\frac{mg{L}_f}{2L}+\frac{m{h}_g{a}_x}{2L}+\frac{m{h}_g{a}_y}{2B_r}\end{array}\right.$$

(7)

where $g$ is the acceleration of gravity, $h_g$ is the distance from ground to centroid, and $a_x$ and $a_y$ are the longitudinal and lateral acceleration of the vehicle, respectively, which can be described as

$$a_x=\frac{1}{m}\left(F_{xfl}+F_{xfr}\right)\cos \delta +\frac{1}{m}\left(F_{xrl}+F_{xrr}\right)-\frac{1}{m}\left(F_{yfl}+F_{yfr}\right)\sin \delta$$

(8)

$$a_y=\frac{1}{m}\left(F_{xfl}+F_{xfr}\right)\sin \delta +\frac{1}{m}\left(F_{yrl}+F_{yrr}\right)+\frac{1}{m}\left(F_{yfl}+F_{yfr}\right)\cos \delta$$

(9)

2.2. Tire Model

The mechanical properties of tires play a vital role in the analysis of vehicle handling stability, so building a reasonable tire model is the basis for studying vehicle dynamics control. The magic formula tire model (MF Tire Model) commonly used in vehicle dynamics research is empirical. This model can better represent the nonlinear characteristics of tires, and its application is relatively mature. Some professional automobile simulation software can quickly obtain the relevant fitting parameters. Therefore, for the tire modeling in this work, we selected the MF Tire Model.

The MF Tire Model is obtained based on the recursive analysis of experimental data. The general expression is as follows:

$$y(x)=D\sin \{C\arctan \left[B\left(x+S_h\right)(1-E)\right. \left.\left.+E\arctan \left(B\left(x+S_h\right)\right)\right]\right\}+S_v$$

(10)

where $B,C,D,E,S_h,S_v$ are characteristic parameters.

The longitudinal force of the tire $F_{xij}$ is calculated as follows:

$$F_{xij}\left(S_{ij}\right)=D\sin \left\{C\arctan \left[B\left(S_{ij}+S_h\right)(1-E)\right.\right.\left.\left.+E\arctan \left(B\left(S_{ij}+S_h\right)\right)\right]\right\}+S_v$$

(11)

The lateral force of the tire $F_{yij}$ is calculated as follows:

$$F_{yij}\left({\alpha }_{ij}\right)=D\sin \{C\arctan \left[B\left({\alpha }_{ij}+S_h\right)(1-E)\right.\left.\left.+E\arctan \left(B\left({\alpha }_{ij}+S_h\right)\right)\right]\right\}+S_v$$

(12)

The inputs of the tire model consist of each tire slip ratio $S_{ij}$ and slip angle ${\alpha }_{ij}$. The slip rate $S_{ij}$ can be calculated by the wheel angular velocity ${\omega }_{ij}$ and the reference speed of the four wheels $v_{ij}$, as shown in the following equation:

$$S_{ij}=\frac{v_{ij}-{\omega }_{ij}R}{v_{ij}}$$

(13)

In addition, the input parameter slip angle ${\alpha }_{ij}$ is calculated as follows:

$$\begin{array}{l}{\alpha }_{fl}=\arctan \left(\frac{v_y+{\omega }_z{L}_f}{v_x-{\omega }_z{B}_f}\right)-\delta \\ {\alpha }_{fr}=\arctan \left(\frac{v_y+{\omega }_{\gamma }{L}_f}{v_x+{\omega }_{\gamma }{B}_f}\right)-\delta \\ {\alpha }_{rl}=\arctan \left(\frac{v_y-{\omega }_{\gamma }{L}_r}{v_x-{\omega }_{\gamma }{B}_f}\right)\\ {\alpha }_{rr}=\arctan \left(\frac{v_y-{\omega }_{\gamma }{L}_r}{v_x+{\omega }_{\gamma }{B}_f}\right)\end{array}$$

(14)

Under combined braking (driving) and steering conditions, the lateral and longitudinal forces should meet the following correction conditions:

$$\left\{\begin{array}{l}{F}_{xij}=\frac{\left|{\sigma }_{xij}\right|}{{\sigma }_{ij}}y(x),F_{yij}=\frac{\left|{\sigma }_{yij}\right|}{{\sigma }_{ij}}y(x)\\ {\sigma }_{ij}=\sqrt{{\sigma }_{xij}^2+{\sigma }_{yij}^2},{\sigma }_{xij}=\frac{S_{ij}}{1+S_{ij}},{\sigma }_{yij}=\frac{\tan {\alpha }_{ij}}{1+S_{ij}}\end{array}\right.$$

(15)

2.3. 2DOF Vehicle Reference Model

The 2DOF vehicle model shown in Figure 3, called the linear bicycle model, is employed here to calculate the ideal yaw rate and sideslip angle.

The dynamic equation of the lateral motion and yaw motion can be expressed as follows:

$$\left\{\begin{array}{c}m\left({\dot{v}}_y+m\cdot v_x\right)=\left(C_f+C_r\right)\beta +\frac{1}{v_y}\left(L_f{C}_f-L_r{C}_r\right){\omega }_{\gamma }-C_f\delta \\ I_z\cdot {\dot{\omega }}_{\gamma }=\left(L_f{C}_f-L_r{C}_r\right)\beta +\frac{1}{v_x}\left(L_f^2{C}_f-L_r^2{C}_r\right){\omega }_{\gamma }-L_f{C}_f\delta \end{array}\right.$$

(16)

where $C_f$ and $C_r$ are the front and rear tire cornering stiffness, respectively.

The 2DOF vehicle model does not consider the effects of tire nonlinearity, steering, or suspension. The car is driving in the linear region, which can better describe the stable state of the vehicle and the driver’s driving intention. The linear 2DOF vehicle model is always regarded as the reference model, and the ideal values of parameters such as the yaw rate and sideslip angle can be obtained. Therefore, the 2DOF vehicle model needs to be analyzed.

When the car is driving in a steady state, the yaw rate and the sideslip angle should satisfy the following equation:

$$\left\{\begin{array}{l}{\dot{v}}_y=0\\ {\dot{\omega }}_{\gamma }=0\end{array}\right.$$

(17)

According to vehicle dynamics, the arctangent value of the ratio of lateral velocity to longitudinal velocity is defined as the sideslip angle. Because the sideslip angle of vehicle motion is usually in a small range, the sideslip angle can be approximately expressed as

$$\beta =\frac{v_y}{v_x}$$

(18)

Substituting Equations (17) and (18) into (16), the ideal yaw rate ${\omega }_d$ and sideslip angle ${\beta }_d$ can be expressed as

$${\omega }_{\gamma d}=\frac{v_x/L}{1+K{v}_x^2}\delta$$

(19)

$${\beta }_d=\frac{v_x\delta }{L\left(1+K{v}_x^2\right)}\left(\frac{L_y}{v_x^2}+\frac{m{L}_f}{C_{\gamma }L}\right)v_x$$

(20)

where $K=\frac{m}{L^2}\left(\frac{L_p}{C_{\gamma }}-\frac{L_r}{C_f}\right)$, representing the stability factor of the vehicle’s steady-state response.

Substituting Equation (19) into (20), we can obtain the following equation:

$${\beta }_{\mathrm{d}}={\omega }_{\gamma }{}_{\mathrm{d}}\left(\frac{L_r}{v_x^2}+\frac{m{L}_f}{C_{r}L}\right)v_x$$

(21)

The lateral acceleration of the vehicle in a constant-speed circular motion is affected by the road adhesion coefficient $\mu$. The lateral acceleration in the limit steady state must meet

$$a_y=\frac{v_x^2}{R_r}=\frac{\left(R_r{\omega }_{\gamma }{}_d\right)v_x}{R_r}={\omega }_{\gamma }{}_d{v}_x\le \mu g$$

(22)

From Equation (22), the maximum expected yaw rate in the limit steady state is expressed as

$${\omega }_{\gamma }{}_d=\frac{\mu g}{v_x}$$

(23)

Combining Equations (19) and (23), we can calculate the ideal yaw rate:

$${\omega }_{\gamma }{}_d=\min \left\{\left|\frac{v_x/L}{\left(1+K{v}_x^2\right)}\delta \right|,\left|\frac{\mu g}{v_x}\right|\right\}sign(\delta )$$

(24)

Combining Equations (21) and (23), the maximum sideslip angle in the limit steady state is expressed as

$${\beta }_d=\mu g\left(\frac{L_r}{v_x{}^2}+\frac{m{L}_f}{C_{r}L}\right)$$

(25)

Combining Equations (21) and (25), we can obtain the ideal sideslip angle:

$${\beta }_d=\min \left\{\left|{\omega }_{\gamma }{}_d\left(\frac{L_r}{v_x^2}+\frac{m{L}_f}{C_{r}L}\right)v_x\right|,\left|\mu g\left(\frac{L_r}{v_x^2}+\frac{m{L}_f}{C_{r}L}\right)\right|\right\}sign(\delta )$$

(26)

2.4. Problem Formulation

The DYC control system designed in this paper is divided into upper and lower layers. The upper controller consists of a sliding mode controller and a speed tracking PI controller, as shown in Figure 4.

Based on the 7DOF vehicle model, the sliding mode controller was selected according to two control parameters: yaw rate and sideslip angle. The vehicle state parameters obtained from the 7DOF and 2DOF ideal reference models are compared. Then, the SMC algorithm controls the difference between the yaw rate and the sideslip angle. Finally, the magnitude of the additional yaw moment required to maintain the vehicle’s stability is calculated.

The speed tracking PI controller obtains the total demand force $F_x$ required to drive the vehicle at a constant speed through the PI control. Based on the difference $\Delta v_x$ between the live vehicle speed $v_x$ and the expected initial vehicle speed $v_{des}$, the total demand force $F_x$ for maintaining vehicle movement is calculated as follows:

$$\left\{\begin{array}{l}\Delta v_x=v_x-v_{des}\\ F_x=k_p\Delta v_x+k_i{\displaystyle \int \Delta v_x}dt\end{array}\right.$$

(27)

where $k_p$ and $k_i$ are the gain coefficient.

For the lower controller, the objective function of torque optimization distribution is first established, and then the constraints of the torque optimization distribution process are comprehensively considered. There are mainly additional yaw moment constraints, road adhesion conditions constraints, vehicle power demand constraints, and motor maximum torque constraints. Finally, the calculated additional yaw moment required to keep the vehicle stable is reasonably distributed on each wheel.

3. Upper Controller

Generally, the design of SMC includes the following two factors: ① the sliding surface function is designed so that the state trajectory of the system has suitable dynamic characteristics such as asymptotic stability after entering the sliding mode, and ② the sliding mode reaching law means the system state trajectories are driven onto the sliding surface for a finite time and are maintained in motion on it.

From the motion equation of the vehicle model, the lateral velocity, longitudinal velocity, and yaw angular velocity are selected to participate in the controller design. Assuming that the longitudinal vehicle velocity is constant, taking the derivative of this Equation (18), we can obtain the following equation:

$$\dot{\beta }=\frac{{\dot{v}}_y}{v_x}$$

(28)

Combining Equations (1)–(3) and (28), the differential equations of the yaw rate and sideslip angle can be obtained.

$$\dot{\beta }=\frac{\left(F_{xfl}+F_{xfr}\right)\sin \delta +\left(F_{yfl}+F_{y_{fr}}\right)\cos \delta +F_{yrl}+F_{yrr}}{m{v}_x}-{\omega }_{\gamma }$$

(29)

$${\dot{\omega }}_{\gamma }=\frac{L_f\left(F_{yfl}+F_{yfr}\right)\cos \delta +\frac{B_f}{2}\left(F_{yfl}-F_{y_{fr}}\right)\sin \delta -\left(F_{yrl}+F_{yrr}\right)L_r+M_z}{I_z}$$

(30)

3.1. Designing of Integral Sliding Mode Surface

The sliding surface of the linear structure can fully meet the design requirements of the control performance of the linear system, which makes the stability analysis simple and convenient when the system is in the sliding mode, and the parameter design is also easy. However, its limitation is that when the linear sliding surface is used, the state tracking error will not converge to zero in a finite time. Therefore, the linear sliding mode surface is suitable for nonlinear systems where the speed and accuracy requirements are not very high. Adding an integral term to the sliding mode surface can significantly improve the approaching process tracking accuracy.

The first step in the SMC design is the sliding mode surface selection. Once the sliding surface is determined, the sliding mode motion’s stability and dynamic quality are determined.

The integral sliding surface can be chosen as

$$s=c_1({\omega }_{\gamma }-{\omega }_{\gamma }{}_d)+{\displaystyle {\int }_0^t({\omega }_{\gamma }-{\omega }_{\gamma }{}_d)}dt+c_2(\beta -{\beta }_d)+{\displaystyle {\int }_0^t(\beta -{\beta }_d)}dt$$

(31)

where $c_1>0,c_2>0$, representing the weight coefficient, which can reflect the proportion of yaw rate and sideslip angle.

Taking the derivative of Equation (31), we can obtain the following equation:

$$\dot{s}=c_1\left({\dot{\omega }}_{\gamma }-{\dot{\omega }}_{\gamma d}\right)+\left({\omega }_{\gamma }-{\omega }_{\gamma d}\right)+c_2\left(\dot{\beta }-{\dot{\beta }}_d\right)+\left(\beta -{\beta }_d\right)$$

(32)

Substituting (29) and (30) into (32), the equation can be expressed as follows:

$$\begin{array}{l}\dot{s}=c_1\frac{L_f\left(F_{yfl}+F_{yfr}\right)\cos \delta +\frac{B}{2}\left(F_{yf}l-F_{yfr}\right)\sin \delta -\left(F_{yrl}+F_{yrr}\right)L_r}{I_z}+c_1\frac{M_z}{I_z}-c_1{\dot{\omega }}_{\gamma d}+{\omega }_{\gamma }-{\omega }_{\gamma d}\\ +c_2\frac{\left(F_{xfl}+F_{xfr}\right)\sin \delta +\left(F_{yfl}+F_{yfr}\right)\cos \delta +F_{yrl}+F_{yrr}}{m{v}_x}-{\omega }_{\gamma }-c_2{\dot{\beta }}_d+\beta -{\beta }_d\end{array}$$

(33)

3.2. Designing of Sliding Mode Reaching Law

A. Conventional Sliding Mode Reaching Law

The reaching law determines the system’s quality in the normal motion phase from outside the sliding surface into the sliding surface. By choosing different reaching laws, different dynamic quality characteristics can be obtained. The forms of the early reaching law include isokinetic reaching law, exponential reaching law, general reaching law, and power reaching law. Among them, the isokinetic and exponential reaching laws are the most commonly used.

The conventional sliding mode reaching law (CSMRL) can widely be expressed as

$$\dot{s}=-{\eta }_1\cdot \mathrm{sgn}(s)-{\eta }_{2}s$$

(34)

where ${\eta }_1>0, {\eta }_2>0$.

It is not difficult to prove that the approach rate (34) can meet the sliding mode reaching condition. The reaching law consists of two parts: $-{\eta }_1\cdot \mathrm{sgn}(s)$ is the constant rate term, and $-{\eta }_{2}s$ is the pure exponential term. When the system state point is far from the sliding mode surface, the reaching law is mainly determined by the pure exponential term. When the motion reaches the sliding mode surface, the pure exponential term is zero, and the constant rate term determines the reaching law. The constant rate reaching law is a part of the traditional exponential reaching law. It is known that the trajectory of approaching motion is banded under the action of symbolic function, and the system state variables cannot be stable at the origin but always chatter near the origin. The system states are unable to eventually converge to the equilibrium point, which indicates that yaw rate and sideslip angle cannot eventually reach the ideal value and the steady-state characteristics of the system are not ideal, so the isokinetic rate reaching law must be analyzed as follows:

$$\dot{s}=-{\eta }_1\cdot \mathrm{sgn}(s)$$

(35)

where $t_1$ is the required time for the system using the reaching law (35) to reach the sliding mode surface, and $s(t_1)=0$. Then, integrating Equation (35) from $0$ to $t_1$, the reaching time is derived as

$$t_1=\left\{\begin{array}{l}\frac{s(0)}{{\eta }_1} s(0)\ge 0\\ \frac{-s(0)}{{\eta }_1} s(0)<0\end{array}\right.$$

(36)

The expression of $t_1$ can be rewritten as follows:

$$t_1=\frac{\left|s\left(0\right)\right|}{{\eta }_1}$$

(37)

B. New Sliding Mode Reaching Law

From Equation (37), we can see that as ${\eta }_1$ increases, the time to reach the sliding mode surface becomes shorter, and the system’s robustness is enhanced. However, under the action of the sign function, increasing ${\eta }_1$ will bring high-frequency chattering. To solve the problem, a new sliding mode reaching law (NSMRL) is given as

$$\dot{s}=-{\eta }_{1}eq(x,s)\mathrm{sgn}(s)-{\eta }_{2}s$$

(38)

$$eq(x,s)=\frac{1}{\epsilon +\left(1+\frac{1}{{\tan }^n(|x_1|+\left|x_2\right|)}-\epsilon \right)e^{-\rho \left|s\right|}}$$

(39)

$$\left\{\begin{array}{l}{x}_1=e_{\omega }={\omega }_{\gamma }-{\omega }_{\gamma d}\\ x_2=e_{\beta }=\beta -{\beta }_d\end{array}\right.$$

(40)

where $\rho >1,0<\epsilon <1$, $n$ is even and $n>1,$ and $x={\left[\begin{array}{cc}{x}_1& x_2\end{array}\right]}^T$ is the state variable of the DYC system.

When the system is far from the sliding surface, and $\left|s\right|$ is large enough, $e^{-\rho \left|s\right|}$ will be very small, and the $eq(x,s)$ converges to $\frac{1}{\epsilon }$. The system will quickly approach the sliding mode surface under the joint action of variable speed term $-{\eta }_{1}eq(x,s)\mathrm{sgn}(s)$ and exponential term $-{\eta }_{2}s$. The speed is higher than the conventional reaching rate. When the system is close to the sliding surface, and the values of $\left|s\right|$ and $x$ are near 0, the $eq(x,s)$ converges to $\frac{{\tan }^n(|x_1|+\left|x_2\right|)}{1+{\tan }^n(|x_1|+\left|x_2\right|)}$. At this moment, the exponential term is small and close to 0. The speed change term $-{\eta }_{1}eq(x,s)\mathrm{sgn}(s)$ bears the leading role. Under the action of the sliding mode reaching law, the system approaches the origin. At the same time, due to the continuous reduction in the speed change term, its reduction speed is faster than the conventional reaching law, which realizes the smooth transition with the sliding mode surface and weakens the effects of chattering.

Taking the first term of the NSMRL, we can obtain the following equation:

$$\dot{s}=-\frac{{\eta }_1}{\epsilon +\left(1+\frac{1}{{\tan }^n(|x_1|+\left|x_2\right|)}-\epsilon \right)e^{-\rho \left|s\right|}}\mathrm{sgn}(s)$$

(41)

Conducting quantitative analysis and calculating the time $t_1$ to reach the sliding mode surface, Equation (41) can be written as follows:

$$\dot{s}\left[\epsilon +\left(1+\frac{1}{{\tan }^n(|x_1|+\left|x_2\right|)}-\epsilon \right)e^{-\rho \left|s\right|}\right]=-{\eta }_1\mathrm{sgn}(s)$$

(42)

By integrating Equation (42) from $0$ to $t_2$, we have

$${\displaystyle {\int }_{s(0)}^{s(t_2)}\frac{1}{\mathrm{sgn}(s)}\left[\epsilon +\left(1+\frac{1}{{\tan }^n(|x_1|+\left|x_2\right|)}-\epsilon \right)e^{-\rho \left|s\right|}\right]ds}={\displaystyle {\int }_0^{t_2}-{\eta }_1}dt$$

(43)

The reaching time $t_2$ is derived as

$$t_2=\left\{\begin{array}{l}\frac{1}{{\eta }_1}\left[\epsilon s(0)+\frac{1}{\rho }\left(1+\frac{1}{{\tan }^n(|x_1|+\left|x_2\right|)}-\epsilon \right)\left(1-e^{-\rho s(0)}\right)\right], s(0)\ge 0\\ \frac{1}{{\eta }_1}\left[-\epsilon s(0)+\frac{1}{\rho }\left(1+\frac{1}{{\tan }^n(|x_1|+\left|x_2\right|)}-\epsilon \right)\left(1-e^{\rho s(0)}\right)\right], s(0)<0\end{array}\right.$$

(44)

The expression of $t_2$ can be rewritten as follows:

$$t_2=\frac{1}{{\eta }_1}\left[\epsilon \left|s(0)\right|+\frac{1}{\rho }\left(1+\frac{1}{{\tan }^n(|x_1|+\left|x_2\right|)}-\epsilon \right)\left(1-e^{-\rho s\left|0\right|}\right)\right]$$

(45)

The time difference under the two reaching laws can be expressed as

$$\begin{array}{l}\Delta t=t_2-t_1=\frac{1}{{\eta }_1}\left[\epsilon \left|s(0)\right|+\frac{1}{\rho }\left(1+\frac{1}{{\tan }^n(|x_1|+\left|x_2\right|)}-\epsilon \right)\left(1-e^{-\rho \left|s(0)\right|}\right)\right]-\frac{\left|s(0)\right|}{{\eta }_1}\\ =\frac{1}{{\eta }_1}\left[(\epsilon -1)\left|s(0)\right|+\frac{1}{\rho }\left(1+\frac{1}{{\tan }^n(|x_1|+\left|x_2\right|)}-\epsilon \right)\left(1-e^{-\rho \left|s(0)\right|}\right)\right]\end{array}$$

(46)

Since $\left(1-e^{-\rho s\left|0\right|}\right)<1$, it can be known that

$$\Delta t<\frac{1}{{\eta }_1}\left[(\epsilon -1)\left|s(0)\right|+\frac{1}{\rho }\left(1+\frac{1}{{\tan }^n(|x_1|+\left|x_2\right|)}-\epsilon \right)\right]=Q$$

(47)

Let

$$\begin{array}{l}Q=\frac{1}{{\eta }_1}\left[(\epsilon -1)\left|s(0)\right|+\frac{1}{\rho }\left(1+\frac{1}{{\tan }^n(|x_1|+\left|x_2\right|)}-\epsilon \right)\right]\\ =\frac{1}{\rho {\eta }_1}\left[\rho (\epsilon -1)\left|s(0)\right|+\left(1+\frac{1}{{\tan }^n(|x_1|+\left|x_2\right|)}-\epsilon \right)\right]\end{array}$$

(48)

Here, $x_1,x_2\ne 0$ can be satisfied in that the value is always between 0 and $t$, if parameter $\rho$ is chosen as

$$\rho \left((\epsilon -1)\left|s(0)\right|\right)\gg \left(1+\frac{1}{{\tan }^n(|x_1|+\left|x_2\right|)}-\epsilon \right)$$

(49)

The equation $Q$ can be presented in a simplified form as

$$Q=(\epsilon -1)\frac{\left|s(0)\right|}{{\eta }_1}$$

(50)

Therefore, according to (49) and inequality (50), we can obtain

$$\Delta t<(\epsilon -1)\frac{\left|s(0)\right|}{{\eta }_1}$$

(51)

Note that the term $\frac{\left|s(0)\right|}{{\eta }_1}$ is strictly positive and $0<\epsilon <1$; then,$\Delta t<0$ can be guaranteed.

It can be seen that under the action of the new exponential term, the system state can reach the sliding surface in a shorter time than the constant rate term. In summary, when away from the sliding mode surface, the system state approach speed is composed of two parts: the new exponential term and the pure exponential term, and the speed is greater than that under the traditional exponential approach law. When approaching the sliding mode surface, the pure exponential term is zero. The rate of the new exponential term is gradually reduced, which can avoid the large inertia of the system state variable under the action of high speed and intense chattering on the sliding mode switching surface, ensuring the smoothness of the system control process.

The Lyapunov function proves the stability of the proposed NSMRL. First, the Lyapunov function is defined as

$$V=\frac{1}{2}{s}^2$$

(52)

The derivative of Equation (52) can be described as

$$\dot{V}=s\dot{s}$$

(53)

Substituting Equations (38)–(40) into (53), the equation can be expressed as follows:

$$\begin{array}{ll}\dot{V}& =s[-{\eta }_{1}eq(x,s)\mathrm{sgn}(s)-{\eta }_{2}s]\\ & =-{\eta }_1\left|s\right|\frac{1}{\epsilon +\left(1+\frac{1}{{\tan }^n(|x_1|+\left|x_2\right|)}-\epsilon \right)e^{-\rho \left|s\right|}}-{\eta }_2{s}^2\end{array}$$

(54)

In Equation (54), because ${\eta }_1>0,{\eta }_2>0$, and $1+1/{\tan }^n(|x_1|+\left|x_2\right|)-\epsilon >0$, $\dot{V}\le 0$ can be established. Therefore, according to the stability judgment conditions of the Lyapunov function, the proposed controller based on NSMRL can satisfy the sliding mode reaching condition.

According to Equations (33) and (38)-(40), the yaw moment controller is designed as

$$\begin{array}{ll}{M}_z& =-\frac{I_z}{c_1}\{\frac{{\eta }_1}{\epsilon +\left(1+\frac{1}{{\tan }^n(|x_1|+\left|x_2\right|)}-\epsilon \right)e^{-\rho \left|s\right|}}\mathrm{sgn}(s)+{\eta }_{2}s+\\ & c_1\frac{L_f\left(F_{y_{fl}}+F_{y_{fr}}\right)\cos \delta +\frac{B}{2}\left(F_{y_{fl}}-F_{y_{fr}}\right)\sin \delta -\left(F_{y_{rl}}+F_{y_{rr}}\right)L_r}{I_z}+\\ & c_1\frac{M_z}{I_z}-c_1{\dot{\omega }}_{\gamma d}+{\omega }_{\gamma }-{\omega }_{\gamma d}+c_2\frac{\left(F_{x_{fl}}+F_{x_{fr}}\right)\sin \delta +\left(F_{y_{fl}}+F_{y_{fr}}\right)\cos \delta +F_{y_{rl}}+F_{y_{rr}}}{m{v}_x}-{\omega }_{\gamma }-c_2{\dot{\beta }}_d+\beta -{\beta }_d\}\end{array}$$

(55)

Based on the above analysis, the designed ISMC can guarantee the DYC system to reach the sliding mode surface in a limited time and be stable on the sliding surface.

4. Lower Controller

The primary purpose of the lower controller is to reasonably distribute the additional yaw moment and total demand force obtained from the previous session to the four wheels. A comparison of the average distribution torque methods is also made here to verify the optimal distribution strategy designed.

4.1. Average Distribution Torque (ADT)

The ADT method is relatively simple so that the total longitudinal force $F_x$ and additional yaw moment $M_z$ are equally distributed on each wheel. The individual tire force can be obtained as follows:

$$\begin{array}{l}{F}_{xfl}=\frac{F_x}{4}-\frac{R}{4}\cdot \frac{M_{\mathrm{z}}}{(B_f+B_r+a\sin \delta )}\\ F_{xfr}=\frac{F_x}{4}+\frac{R}{4}\cdot \frac{M_{\mathrm{z}}}{(B_f+B_r+a\sin \delta )}\\ F_{xrl}=\frac{F_X}{4}-\frac{R}{4}\cdot \frac{M_{\mathrm{z}}}{(B_f+B_r+a\sin \delta )}\\ F_{xrr}=\frac{F_x}{4}+\frac{R}{4}\cdot \frac{M_{\mathrm{z}}}{(B_f+B_r+a\sin \delta )}\end{array}$$

(56)

Then, the motor torque $T_{dij}$ required by each wheel can be expressed as

$$T_{dij}=F_{xij}\cdot R$$

(57)

4.2. Optimal Distribution Torque (ODT)

The disadvantage of ADT is that the yaw moment is uniformly distributed on the four wheels according to the average rule. However, ODT is an integrated torque distribution method, and it takes into account the three factors of tire utilization rate, motor output torque, and road adhesion conditions. Thus, it can better ensure the stability of the driving process under extreme road conditions.

In optimal control allocation, the basic problem is defining the objective function and selecting the appropriate optimization target. In this paper, the concept of “tire utilization rate” is first introduced, which refers to the sum of four tires’ ratio of road adhesion on a single wheel to the maximum adhesion it can obtain. When this value is close to 1, it means that the working ability of the tire is close to its limit and cannot provide greater adhesion. Without control at this point, the vehicle is very likely to lose stability. The mathematical definition of a single tire attachment utilization rate is given as follows:

$$\chi =\Sigma \frac{F_{xi}^2+F_{yi}^2}{{\left(\mu F_{zi}\right)}^2},i=fl,fr,rl,rr$$

(58)

where $\mu$ is the road adhesion coefficient.

Masato Abe proposed that the minimum sum of squares of all tire utilization rates should be taken as the objective function to distribute the force on each tire to ensure that all wheels can maintain stability and have a certain stability margin [41]. The formula with the minimum sum of squares of the utilization rate of four tires as the objective function is shown as

$$\min J={\displaystyle \sum _{i=fl,fr,rl,rr}\frac{F_{xi}^2+F_{yi}^2}{{\left(\mu F_{zi}\right)}^2}}$$

(59)

For 4WID electric vehicles, the calculation is complicated when considering the optimization under the combined action of longitudinal force, transverse force, and vertical force, and the online real-time optimization cannot be realized. Vehicles with high real-time requirements will have great safety risks. Therefore, to increase the timeliness of vehicle control, the objective function of this paper is simplified to consider only the control longitudinal force:

$$\min J={\displaystyle \sum _{i=fl,fr,rl,rr}\frac{F_{xi}^2}{{\left(\mu F_{zi}\right)}^2}}$$

(60)

We propose the objective function by combining the tire utilization rate with the standard deviation coefficient. Adding the standard deviation coefficient to the optimization goal ensures that the mean and variance of the target value reach the optimal value simultaneously. The optimization objective function is expressed as

$$\min J={\displaystyle \sum _{i=fl,fr,rl,rr}\frac{F_{xi}^2}{{\left(\mu F_{zi}\right)}^2}}+\frac{\sqrt{\frac{1}{4}{\displaystyle \sum _{i=fl,fr,rl,rr}^{}{\left(\frac{F_{xi}^2}{{\left(\mu F_{zi}\right)}^2}-\frac{1}{4}{\displaystyle \sum _{i=fl,fr,rl,rr}^{}\frac{F_{xi}^2}{{\left(\mu F_{zi}\right)}^2}}\right)}^2}}}{\frac{1}{4}{\displaystyle \sum _{i=fl,fr,rl,rr}^{}\frac{F_{xi}^2}{{\left(\mu F_{zi}\right)}^2}}}$$

(61)

The Lagrange multiplier method in

Table 1. Lagrange multiplier method.

Lagrange Multiplier Method
$\min$ $f(x)$ $s.t$ $h(x)=0$ $g(x)\le 0$  By introducing slack variables $\Omega$, inequality constraints are transformed into equality constraints
$\min L(x,\lambda ,\Omega )=f(x)-{\lambda }_{1}h(x)-{\lambda }_2(g(x)+{\Omega }^2)$  where ${\lambda }_1,{\lambda }_2$ are the Lagrange computing factors, ${\Omega }^2$ is introduced to ensure that it is non-negative, so that the constraint function is less than or equal to zero is satisfied.
The derivative of $L(x,\lambda ,\Omega )$ with respect to each variable is 0 $\frac{\partial L}{\partial x}=0,\frac{\partial L}{\partial {\lambda }_1}=h(x)=0,\frac{\partial L}{\partial {\lambda }_2}=g(x)+{\Omega }^2=0,\frac{\partial L}{\partial \Omega }=-2\Omega$, By solving the equations above, we can obtain $x$.

is used to solve the optimization problem. The equality constraints are as follows:

$${\psi }_1(F_{xfl},F_{xfr},F_{xrl},F_{xrr})=(F_{xfl}\cos \delta +F_{xfr}\cos \delta +F_{xrl}+F_{xrr})-F_x=0$$

(62)

$${\psi }_2(F_{xfl},F_{xfr},F_{xrl},F_{xrr})=[\frac{B_f}{2}(-F_{x_{fl}}\cos \delta +F_{x_{fr}}\cos \delta +F_{x_{rr}}-F_{x_{rl}})+L_f\left(F_{xfl}\sin \delta +F_{x_{fr}}\sin \delta \right)]-M_z=0$$

(63)

First, an auxiliary solution function is constructed as

$$g\left(F_{xfl},F_{xfr},F_{xrl},F_{xrr},{\lambda }_1,{\lambda }_2\right)=J+{\lambda }_1{\psi }_1(F_{xfl},F_{xfr},F_{xrl},F_{xrr})+{\lambda }_2{\psi }_2(F_{xfl},F_{xfr},F_{xrl},F_{xrr})$$

(64)

where ${\lambda }_1,{\lambda }_2$ are the computing factors.

At the same time, the longitudinal force $F_{xij}$ should meet the constraint conditions of motor torque and ground longitudinal adhesion, expressed as

$$\left\{\begin{array}{l}{\psi }_3(F_{xfl},F_{xfr},F_{xrl},F_{xrr})=\left|F_{xi}\right|-\left(T_{\max }/R\right)\le 0,i=fl,fr,rl,rr\\ {\psi }_4(F_{xfl},F_{xfr},F_{xrl},F_{xrr})=F_{xi}^2+F_{yi}^2-{\left(\mu F_{zi}\right)}^2\le 0,i=fl,fr,rl,rr\end{array}\right.$$

(65)

where $T_{\max }$ is the maximum output torque of the motor.

By introducing slack variables ${\Omega }_i(i=1,2)$, inequality constraints are transformed into equality constraints, which can be solved quickly by the Lagrange multiplier method. The new Lagrange function is written as

$$g\left(F_{xfl},F_{xfr},F_{xrl},F_{xrr},{\lambda }_1,{\lambda }_2,{\lambda }_3,{\lambda }_4,{\Omega }_1,{\Omega }_2\right)=J-{\lambda }_1{\psi }_1(F_{xfl},F_{xfr},F_{xrl},F_{xrr})-{\lambda }_2{\psi }_2(F_{xfl},F_{xfr},F_{xrl},F_{xrr})$$

$$-{\lambda }_3({\psi }_3(F_{xfl},F_{xfr},F_{xrl},F_{xrr})+{\Omega }_1{}^2)-{\lambda }_4({\psi }_4(F_{xfl},F_{xfr},F_{xrl},F_{xrr})+{\Omega }_2{}^2)$$

(66)

where ${\lambda }_3,{\lambda }_4$ are the computing factors.

The partial derivative of the Lagrange function in Equation (67) is expressed as

$$\{\begin{array}{l}\frac{\partial g}{\partial F_{xfl}}=\frac{\partial J}{\partial F_{xfl}}-{\lambda }_1\frac{\partial {\psi }_1}{\partial F_{xfl}}-{\lambda }_2\frac{\partial {\psi }_2}{\partial F_{xfl}}-{\lambda }_3\frac{\partial {\psi }_3}{\partial F_{xfl}}-{\lambda }_4\frac{\partial {\psi }_4}{\partial F_{xfl}}=0\\ \frac{\partial g}{\partial F_{xfr}}=\frac{\partial J}{\partial F_{xfr}}-{\lambda }_1\frac{\partial {\psi }_1}{\partial F_{xfr}}-{\lambda }_2\frac{\partial {\psi }_2}{\partial F_{xfr}}-{\lambda }_3\frac{\partial {\psi }_3}{\partial F_{xfr}}-{\lambda }_4\frac{\partial {\psi }_4}{\partial F_{xfr}}=0\\ \frac{\partial g}{\partial F_{xrl}}=\frac{\partial J}{\partial F_{xrl}}-{\lambda }_1\frac{\partial {\psi }_1}{\partial F_{xrl}}-{\lambda }_2\frac{\partial {\psi }_2}{\partial F_{xrl}}-{\lambda }_3\frac{\partial {\psi }_3}{\partial F_{xrl}}-{\lambda }_4\frac{\partial {\psi }_4}{\partial F_{xrl}}=0\\ \frac{\partial g}{\partial F_{xrr}}=\frac{\partial J}{\partial F_{xrr}}-{\lambda }_1\frac{\partial {\psi }_1}{\partial F_{xrr}}-{\lambda }_2\frac{\partial {\psi }_2}{\partial F_{xrr}}-{\lambda }_3\frac{\partial {\psi }_3}{\partial F_{xrr}}-{\lambda }_4\frac{\partial {\psi }_4}{\partial F_{xrr}}=0\\ \frac{\partial g}{\partial {\lambda }_1}={\psi }_1(F_{xfl},F_{xfr},F_{xrl},F_{xrr})=0\\ \frac{\partial g}{\partial {\lambda }_2}={\psi }_2(F_{xfl},F_{xfr},F_{xrl},F_{xrr})=0\\ \frac{\partial g}{\partial {\lambda }_3}=({\psi }_3(F_{xfl},F_{xfr},F_{xrl},F_{xrr})+{\Omega }_1{}^2)=0\\ \frac{\partial g}{\partial {\lambda }_4}=({\psi }_4(F_{xfl},F_{xfr},F_{xrl},F_{xrr})+{\Omega }_2{}^2)=0\\ \frac{\partial g}{\partial {\Omega }_1}=-2{\Omega }_1,\frac{\partial g}{\partial {\Omega }_2}=-2{\Omega }_2\end{array}$$

(67)

The above equation can solve the longitudinal forces. By (57), the required motor torque can be obtained. The Lagrange multiplier method is given in Table 1. The results are solved using the MATLAB software to write programs.

5. Simulation

The simulation was carried out based on MATLAB/Simulink software. The three control strategies are compared under the conditions of a single lane change and J-turn to verify the effectiveness of the proposed control algorithm. It includes the CSMRLADT strategy composed of the conventional sliding mode reaching law and ADT, the NSMRLADT strategy composed of the new sliding mode reaching law and ADT, and the NSMRLODT strategy composed of the new sliding mode reaching law and ODT. The vehicle parameters are shown in

Table 2. Vehicle parameters.

Description	Symbol	Value
Total mass of the car vehicle	$m$	1480 (kg)
Distance from front axle to the center of centroid	$L_f$	1.2 (m)
Distance from rear axle to the center of centroid	$L_r$	1.4 (m)
Yaw moment of inertia of vehicle	$I_z$	1523 (kg$\cdot {\mathrm{m}}^2$)
Moment of inertia of the wheel	$J$	2.1 (kg$\cdot {\mathrm{m}}^2$)
Distance from the center of centroid to ground	$h_g$	0.5 (m)
Front track width	$B_f$	1.6 (m)
Rear track width	$B_r$	1.6 (m)
Rolling radius of the tire	$R$	0.354 (m)
Front wheel cornering stiffness	$C_f$	35796 (N/rad)
Rear wheel cornering stiffness	$C_r$	35400 (N/rad)
Peak torque of the motor tire	$T_{\max }$	400 (N$\cdot \mathrm{m}$)
Rolling resistance coefficient	$f$	0.018

.

5.1. Performance Analysis under Single Lane Change (SLC) Maneuver

The front wheel corner input is given a sine-wave input (3–7 s) with an amplitude of 0.08 rad, as shown in Figure 5. Road friction coefficients are set as 0.5 to simulate wet roads. The starting speed is set to 108 km/h.

The two stability index curves under the condition of SLC are shown in Figure 6. The peak value under control is lower than without control, the changes in the curve are relatively gentle, and the changing trend is more suitable for the SLC maneuver. It shows that the proposed control strategy can improve the vehicle’s stable driving and trajectory-keeping ability. In contrast, NSMRLODT has relatively lower peak values and smaller fluctuations, showing a suitable performance, which indicates that the NSMRLODT method has an adequate control effect.

The simulation comparison results between NSMRLADT and NSMRLODT are shown in Figure 7. The vehicle track is given in Figure 7a. NSMRLADT and NSMRLODT can better track the desired trajectory path. Vehicles with ODT control have better lateral stability compared to ADT control. The $\dot{\beta }-\beta$ phase diagram is shown in Figure 7b. The sideslip angle of the two methods is in the controllable range and presents a closed-loop convergence structure. Still, the overall control effect of NSMRLODT is superior to NSMRLADT. The ADT and ODT are shown in Figure 7c,d, respectively. The torque of the same side wheel under NSMRLADT control is the same, while the NSMRLODT control method can distribute the torque required to change among four respective in-wheel motors. It is possible to prevent the wheel with a small vertical load from quickly reaching the longitudinal force saturation. Although the torque distribution trend of the wheel on the same side under the NSMRLODT control method is the same, the amplitude is different. The front wheel torque is greater than the rear wheel torque, and this is because the front wheel is a steering wheel. As shown in Figure 7e, the sum of the load ratios of each wheel of NSMRLODT is smaller than that of NSMRLADT, indicating that the torque of each wheel is reasonably distributed under NSMRLODT. The maximum stability margin is kept as much as possible under the precondition of meeting the current yaw moment control. Leaving sufficient longitudinal force at the disposal of the car can make the vehicle more secure and reliable.

5.2. Performance Analysis under J-turn Maneuver

To further test the effectiveness of the proposed control strategy, we used the front steering angle for analog input, as shown in Figure 8.

The simulation results are shown in Figure 9, assuming driving at an initial speed of 80 km/h on a slippery road with an adhesion coefficient of 0.3.

The uncontrolled vehicle under the extreme road conditions of J-turn is unstable compared to the sinusoidal road conditions. In particular, the sideslip angle is worse than that of the SLC road condition. Compared with the uncontrolled vehicle’s yaw rate and sideslip curves, the CSMRLADT and NSMRLADT have converged the final tracking results. Still, the poor road conditions make the path tracking accuracy and stability poor. In contrast, the amplitude deviation of the yaw rate and sideslip angle is large, and there is chattering. The NSMRLODT effectively solves the chattering problem. The curves of yaw rate and sideslip angle are smoother and have suitable performance, which shows that the proposed control strategy can significantly improve yaw stability.

The simulation comparison results between NSMRLADT and NSMRLODT are shown in Figure 10. Note that the NSMRLODT introduced here is still effective in improving stability in extreme J-turn maneuvers. Based on the analysis of the above results, NSMRLODT is more suitable for the ideal state than NSMLRLADT.

6. Conclusions

The innovation of this paper is that new algorithms are proposed for the upper and lower controllers. In SMC, a more realistic 7DOF model is adopted, which contains two state variables, yaw rate and sideslip angle. NSMRL aims to reduce chattering and make sliding mode surface state variables converge faster. The lower controller adopts the utilization rate of four tires and its standard deviation coefficient as the objective function and prevents the stability margin gap between tires from being too large. The simulation results compare three hierarchical control strategies, CSMRLADT, NSMRLADT, and NSMRLODT, under single lane change and J-turn conditions. The comparison results show the NSMRLODT algorithm provides better trajectory tracking, and the torque of each wheel is reasonably distributed. The results show that the DYC control strategy designed in this paper can improve the system’s dynamic performance and tracking accuracy. At the same time, the stability and security of the system are considerably improved.