On the Lateral Stability System of Four-Wheel Driven Electric Vehicles Based on Phase Plane Method

Abstract

To improve the handling and stability of four-wheel independent drive electric vehicles (FWID EVs), this paper introduces a hierarchical architecture lateral stability control system. The upper-level controller is responsible for generating the additional yaw moment required by the vehicle. This includes a control strategy based on feedforward control and a Linear Quadratic Regulator (LQR) for handling assistance control, an LQR-based stability control, a PID controller-based speed-following control, and a stability assessment method. The lower-level controller uses Quadratic Programming (QP) to optimally distribute the additional yaw moment to the four wheels. A “normalized” method was proposed to determine vehicle stability. After comparing it with the existing double-line method, diamond method, and curved boundary method through the open-loop Sine with Dwell test and the closed-loop Double Lane Change (DLC)test simulation, the results demonstrate that this method is more sensitive and accurate in determining vehicle stability, significantly enhancing vehicle handling and stability.

1. Introduction

According to a report released by Worldmetrics, as of 2024, there are over 1 billion vehicles operating globally [1]. IEA statistics indicate that global electric vehicle (EV) sales in 2023 approached 14 million units, bringing the cumulative number of EVs on the road to approximately 40 million, which accounts for around 4% of the total vehicle fleet worldwide [2]. It can be said that cars have become an indispensable means of transportation in people’s daily lives. However, with the increase in the number of vehicles, traffic accidents have also risen. According to statistics from the World Health Organization, approximately 1.19 million people die each year due to traffic accidents, with about 30% of these accidents resulting from vehicle instability caused by speeding [3]. To prevent traffic accidents caused by vehicle instability at high speeds, research on vehicle stability control systems has become one of the key research areas in the automotive industry [4,5,6].

The vehicle stability control system improves vehicle handling and stability by continuously monitoring the difference between the vehicle’s real-time operating state and the ideal reference model and applying certain control strategies. One strategy, researched by Li, X. et al. [7], is based on an adaptive sliding mode controller; the sliding surface’s weight is automatically adjusted according to the vehicle’s operating state, which effectively improves vehicle handling and stability. To enhance computational efficiency and reduce calculation time, Roberto, Z. et al. [8] employed a constrained parametric Model Predictive Controller (MPC), applying exponential parameterization to the control vector within the prediction range to reduce the complexity of solving optimization problems. This allows for the fast computation of additional yaw moment to maintain the vehicle’s lateral stability. Zhai, L. et al. [9] used a fuzzy Proportional-Integral-Derivative (PID) controller to calculate the yaw moment required to maintain vehicle stability, effectively improving the vehicle’s stability.

The stability criterion of a vehicle dictates the timing of intervention and withdrawal of the control system, which is an extremely crucial link in stability control. Therefore, research on the stability criterion of a vehicle is of great significance for improving vehicle stability and ensuring the safety of drivers and passengers [6,10,11]. Scholars have conducted extensive research on stability criterion of vehicles. In Reference [12], a phase plane of vehicle sideslip angle $\beta$ and $\dot{\beta }$ was constructed based on different vehicle operating conditions. The diamond method was used to divide the stable regions of the phase plane under various conditions. To avoid excessive data volume, a Support Vector Regression (SVR) model was established to achieve automatic regression of dynamic stability regions. In [13], based on real-time vehicle dynamics and road conditions, the phase plane of sideslip angle $\beta$ and yaw rate $\gamma$ was divided into three control zones: stable, critical, and unstable, enabling timely and precise vehicle stability control. This method adapts to rapid changes in vehicle status under extreme driving conditions and enhances overall vehicle stability. Reference [14] utilized the real-time estimated lateral tire stiffness of the front and rear wheels to obtain the understeer coefficient, which reflects the vehicle’s stability characteristics, to determine the vehicle’s stability. Reference [15] based on a two-degree-of-freedom vehicle model, created a $\beta -\dot{\beta }$ phase plane and applied the “normalization” method, comprehensively considering the effects of road surface adhesion coefficient $\mu$, vehicle speed $v_x$, and front-wheel steering angle $\delta$ on the size and location of the stable region in the phase plane. The stability determination was evaluated by using the distance between $\beta$ and $\gamma$ values from the boundary of their value ranges.

Due to the independent control of the four wheels, FWID EVs offer high flexibility and control precision. This paper conducts research on lateral stability control based on FWID EVs. The aforementioned extensive research has provided valuable references for this study. The hierarchical control framework proposed in References [7,8,9,14] laid the foundation for the controller structure design in this research. Although the different stability criterion methods proposed in References [12,13,16,17] enable a vehicle stability assessment, they ignore the influence of the front-wheel steering angle $\delta$ on the stability center and the size of the stability region in the phase plane when defining the stability boundaries. The “normalization” method in Reference [15] accounts for the influence of the front-wheel steering angle $\delta$, but it does not consider the impact of vertical load changes on the four wheels, caused by vehicle body roll during cornering, on the phase plane. Therefore, based on References [15,18], the main contributions of this paper are as follows: a four-wheel vehicle dynamics model and a Magic Formula Tire model are established to construct the $\beta -\dot{\beta }$ phase plane, which considers the changes in vertical tire load due to vehicle body roll; a “normalization” method is proposed for vehicle stability assessment, and a hierarchical lateral stability control system is designed. When the vehicle is in a stable state, the control system outputs additional yaw moment to improve vehicle handling. When the vehicle is in an unstable state, the control system outputs additional yaw moment to enhance vehicle stability.

The rest of this article is organized as follows: Section 2 introduces the four-wheel vehicle dynamics model that reflects the dynamic characteristics of the vehicle, as well as the two-degree-of-freedom (2-DOF) vehicle dynamics model used as a reference for designing the controller. Section 3 presents the $\beta -\dot{\beta }$ phase plane method and the stability criterion method. Section 4 focuses on the design of the controller, which is divided into an upper-level and lower-level controller. The upper-level controller includes the handling assistance controller, the stability controller, the speed-following controller, and the stability criterion, while the lower-level controller consists of the control allocation. Section 5 introduces the testing experiments and simulation results. Section 6 concludes the paper.

2. Vehicle Dynamics Model

A precise vehicle dynamics model can accurately represent the dynamic characteristics of a vehicle. In this paper, a four-wheel two-degree-of-freedom vehicle dynamics model and a nonlinear Magic Formula Tire model, known for its superior performance, are utilized to characterize the vehicle’s dynamic behavior. When designing the vehicle controller, a linear 2-DOF vehicle dynamics model, which requires less computational effort and can represent vehicle stability and handling performance, is used as the reference model.

2.1. Four-Wheel Vehicle Dynamics Model

In this paper, a four-wheel vehicle dynamics model is adopted to consider lateral load transfer between left-side and right-side wheels during steering. Neglecting the vehicle’s pitch motion, as shown in Figure 1 refined by Zhao et al. [4]. $v$ represents the velocity at the vehicle’s center of mass; $v_x$ and $v_y$ represent the vehicle’s longitudinal and lateral speeds, respectively; $F_{xij}$ and $F_{yij}$ represent the vehicle’s longitudinal and lateral forces acting on the wheels, respectively; ${\alpha }_{ij}$ is the slip angle at each wheel, where $i\in \left\{f,r\right\}$ represents the front or the rear tire and $j\in \left\{l,r\right\}$ represents the left or the right tire; $L$ is the vehicle’s wheelbase, with $l_f$ and $l_r$ indicating the distances from the vehicle’s center of mass to the front and rear axles, respectively; δ is the steering angle of the front axle; $\mathsf{\gamma }$ and $\beta$ represent the yaw rate and side slip angle of the vehicle, respectively; $t_f$ and $t_r$ are the front and rear track widths, respectively. A four-wheel vehicle dynamics model is established, as described in (1), and includes both the vehicle’s lateral and yaw motions.

$$\left\{\begin{array}{c}m{v}_x\left(\dot{\beta }+r \right)=\left(F_{yfl}+F_{yfr}\right)\cos \delta +F_{yrl}+F_{yrr}\hfill \\ I_z\dot{r}=l_f\left(F_{yfl}+F_{yfr}\right)\cos \delta -l_r\left(F_{yrl}+F_{yrr}\right)+{\displaystyle \frac{t_f}{2}}\left(F_{yfl}-F_{yfr}\right)\sin \delta +M_z\hfill \end{array}\right.$$

(1)

where $m$ represents the vehicle mass; $I_z$ is the yaw moment of inertia; $M_z$ represents the additional yaw moment.

Assuming the longitudinal vehicle speed remains constant during operation, and given that the lateral velocity is much smaller than the longitudinal velocity, the vehicle’s sideslip angle $\beta$ can be approximately expressed by (2). The slip angles of the front and rear wheels can be expressed by (3). Considering load transfer during vehicle cornering, the vertical forces $F_{zij}$ acting on each wheel can be expressed by (4). The lateral forces $F_{yij}$ can be described using the Magic Formula Tire model, as shown in (5). For detailed information, refer to [18].

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

(2)

$$\left\{\begin{array}{c}{\alpha }_{fl}={\alpha }_{fr}=\delta -\beta -{\displaystyle \frac{l_{f}r}{v_x}}\\ {\alpha }_{rl}={\alpha }_{rr}=\beta -{\displaystyle \frac{l_{r}r}{v_x}} \end{array}\right.$$

(3)

$$\left\{\begin{array}{c}{F}_{zfl}={\displaystyle \frac{mg{l}_r}{2L}}-{\displaystyle \frac{m_s{a}_{y}h}{t_f}}\hfill \\ \begin{array}{c}{F}_{zfr}={\displaystyle \frac{mg{l}_r}{2L}}+{\displaystyle \frac{m_s{a}_{y}h}{t_f}}\hfill \\ \begin{array}{c}{F}_{zrl}={\displaystyle \frac{mg{l}_f}{2L}}-{\displaystyle \frac{m_s{a}_{y}h}{t_r}}\\ F_{zrr}={\displaystyle \frac{mg{l}_f}{2L}}+{\displaystyle \frac{m_s{a}_{y}h}{t_r}}\end{array}\hfill \end{array}\hfill \end{array}\right.$$

(4)

$$\left\{\begin{array}{c}{F}_{yfj}=f\left({\alpha }_{fj},\mu ,F_{zfj}\right) \hfill \\ F_{yrj}=f\left({\alpha }_{rj},\mu ,F_{zrj}\right)\hfill \end{array}\right.$$

(5)

where $m_s$ represents the sprung mass of the vehicle; $h$ is the height of the vehicle’s center of gravity; $a_y$ is the vehicle’s lateral acceleration and $\mu$ is the road adhesion coefficient.

2.2. Linear 2-DOF Vehicle Dynamics Model

By combining the left and right wheels of Figure 1, a linear 2-DOF vehicle dynamics model is used, as shown in Figure 2. In order to simplify calculations, this model uses the linear tire model to calculate the lateral force as shown in (6), $C_{\alpha f}$ and $C_{\alpha r}$ represent the cornering stiffness of the front and rear wheels, respectively. Using (1), (3), and (6) the state-space representation of this linear 2-DOF vehicle dynamics model can be expressed in (7).

$$\left\{\begin{array}{c}{F}_{yf}=C_{\alpha f}\ast {\alpha }_f\hfill \\ F_{yr}=C_{\alpha r}\ast {\alpha }_r\hfill \end{array}\right.$$

(6)

$$\dot{\mathit{x}}=\mathit{A}\mathit{x}+\mathit{B}u+\mathit{G}w$$

(7)

where $\mathit{x}=\left[\begin{array}{c}\beta \\ r\end{array}\right]$, $\mathit{A}=\left[\begin{array}{cc}{a}_{11}& a_{12}\\ a_{21}& a_{22}\end{array}\right]$, $\mathit{B}=\left[\begin{array}{c}{b}_1\\ b_2\end{array}\right]$, $a_{11}=-{\displaystyle \frac{\left(C_{\alpha f}+C_{\alpha r}\right)}{m{v}_x}}$, $a_{12}$ = ${\displaystyle \frac{l_r{C}_{\alpha r}-l_f{C}_{\alpha f}}{m_{vx}^2}}-1$, $a_{21}={\displaystyle \frac{l_r{C}_{\alpha r}-l_f{C}_{\alpha f}}{I_z}}$, $a_{22}=-{\displaystyle \frac{\left(l_f^2{C}_{\alpha f}+l_r^2{C}_{\alpha r}\right)}{I_z{v}_x}}$, $b_1=0$, $b_2={\displaystyle \frac{1}{I_z}}$, $\mathit{G}=\left[\begin{array}{c}{\displaystyle \frac{C_{\mathsf{\alpha }f}}{m{v}_x}}\\ {\displaystyle \frac{l_f{C}_{\mathsf{\alpha }f}}{I_z}}\end{array}\right]$ $u=M_z$, $w=\delta$

When $\dot{\beta }=0$ and $\dot{\gamma }=0$, meaning the vehicle is in a stable state, the ideal vehicle sideslip angle ${\beta }_d$ and yaw rate ${\gamma }_d$ that serve as references for the controller can be obtained.

3. Vehicle Stability Criterion

This paper uses the phase plane method to determine vehicle stability. To improve the accuracy of the phase plane, a four-wheel vehicle dynamics model, which accounts for the changes in tire vertical forces caused by vehicle roll, along with the nonlinear Magic Formula Tire model, is employed to establish the $\beta -\dot{\beta }$ phase plane.

3.1. The Establishment of the Phase Plane

By combining (1), (3), (4), and (5), the vehicle state–space expression (8) is derived. When $M_z=0$, $v_x=80 \mathrm{km}/\mathrm{h}$ and $\mu =0.85$, with $\delta$ set to 0 $deg$, 1.5 $deg$, and 3 $deg, respectively,$ the differential solution of (8) yields the $\beta -\dot{\beta }$ phase plane diagram shown in Figure 3. It can be seen that $\delta$ has a certain effect on the position of the stable center and the stable region in the phase plane. As $\delta$ increases, the stable center shifts a certain distance to the left along the horizontal axis, the stable region on the left shrinks, and the stable region on the right expands. When $\delta$ increases beyond a certain value, the stable center disappears, and the system becomes unstable.

$$\left[\begin{array}{c}\dot{\beta }\\ \dot{\gamma }\end{array}\right]=\left[\begin{array}{c}{\displaystyle \frac{\left(F_{yfl}+F_{yfr}\right)·\cos \delta +\left(F_{yrl}+F_{yrr}\right)}{m{v}_x}}-\gamma \\ {\displaystyle \frac{\left(F_{yfl}+F_{yfr}\right)·l_f·\cos \delta -l_r·\left(F_{yrl}+F_{yrr}\right)+\frac{t_f}{2}\left(F_{yfl}-F_{yfr}\right)\sin \delta }{I_z}}+{\displaystyle \frac{M_z}{I_z}}\end{array}\right]$$

(8)

3.2. Phase Plane-Based Stability Region Estimation

Existing methods for determining vehicle stability based on the phase plane primarily involve dividing the stable region. When the point ($\beta$, $\dot{\beta }$) lies within the stable region, the vehicle is considered stable; otherwise, it is considered unstable. Since the front wheel steering angle $\delta$ changes within a small range when the vehicle is traveling at high speeds, its influence is typically ignored. These methods mainly focus on the effects of $v_x$ and $\mu$, which have a greater impact on the phase plane and stable region. Relatively little in-depth research has been conducted on the effect of the $\delta$ on the phase plane; however, the influence of the $\delta$ on the phase plane cannot be ignored [11].

3.2.1. The Typical Phase Plane Stable Region Determination Method

As shown in Figure 4, Case A, Case B, and Case C correspond to the double-line method, the diamond method, and the curved boundary method, respectively. The double-line method is shown in Figure 4a and can be specifically referenced in [16,17]; the diamond method is shown in Figure 4b and can be specifically referenced in [12]; and the curved boundary method fully considers the maximum allowable $\gamma$ under the constraint of $\mu$ (12) and the two extreme values of $\beta$ (${\beta }_{min}$, ${\beta }_{max}$) when $\delta =0$, obtained from the first equation of (8), to constrain the phase plane, resulting in the magenta curve in Figure 4c. The area outside the thick solid line in Figure 4 represents the instability zone, the area enclosed by the black dashed line represents the stable zone, and the area between the black dashed line and the thick solid line is the critical zone.

By input $\dot{\beta }=0$ into (8), $\gamma$ can be expressed as (9):

$$\gamma ={\displaystyle \frac{\left(F_{yfl}+F_{yfr}\right)·\cos \delta +\left(F_{yrl}+F_{yrr}\right)}{m{v}_x}}$$

(9)

While the vehicle is turning, ignoring the velocity changes along the longitudinal axis, combined with $\dot{\beta }=0$, then $\gamma$ can be presented as (10).

$$\gamma ={\displaystyle \frac{a_y}{v_x}}$$

(10)

Restricted by the road adhesion coefficient $\mu$, the lateral acceleration is given as shown in (11):

$$a_y\le 0.85 \mu g$$

(11)

Considering the vehicle turning left and right, the extreme values of the vehicle’s yaw rate $\gamma$ can be obtained as shown in (12).

$$\left\{\begin{array}{c}{\gamma }_{min}=-{\displaystyle \frac{0.85 \mu g}{v_x}}\\ {\gamma }_{max}={\displaystyle \frac{0.85 \mu g}{v_x}} \end{array}\right.$$

(12)

When the point ($\beta$, $\dot{\beta }$) is located in the stable region, the stability control weight coefficient $W=0$; when it is in the unstable region, $W=1$; and when it is in the critical region, to achieve a smooth transition between the stable and unstable regions, the expression for the stability control weight coefficient $W$ is given as shown in (13). $s_1$ is the distance from the point ($\beta$, $\dot{\beta }$) to the black dashed line, and $s_0$ is the distance between the corresponding solid boundary and the dashed boundary.

$$W={\displaystyle \frac{s_1}{s_0}}$$

(13)

3.2.2. Effect of $\delta$ on the Phase Plane

As discussed in Section 3.1, the influence of $\delta$ on the $\beta -\dot{\beta }$ phase plane cannot be ignored. When $\delta$ takes values of −2 $deg$, 0 $deg$, and 2 $deg$, respectively, and the Case C method is applied, the local stable phase plane diagram is shown in Figure 5. It is evident that the stable region and stable center in the phase plane vary significantly with different $\delta$. Ignoring the influence of the front wheel steering angle $\delta$ may lead to inaccurate vehicle stability determinations or even misjudgments.

3.2.3. The “Normalization” Method Design

To improve the accuracy of the vehicle stability determination, this paper adopts the “normalization” method as Case D for the vehicle stability assessment. This method comprehensively considers the effects of $\mu$, $v_x$, and $\delta$ on the size and position of the stable region in the phase plane, avoiding the misjudgments caused by using a fixed partitioning method to determine whether the real-time changing midpoint ($\dot{\beta }$, $\beta$) in the phase plane is stable. This approach transforms the stability determination problem from assessing whether the point ($\dot{\beta }$, $\beta$) lies within a specific region in the phase plane to evaluating the distance between the values of $\beta$ and $\gamma$ and the boundaries of their value ranges as a metric for stability.

Under different driving conditions, the value range of $\gamma$ can be determined by (12), and the value range of $\beta$ can be calculated using (5), (9), and (12). When $\delta$ reaches a certain value, causing the phase plane to have no stable center, it is defined that ${\beta }_{min}$= ${\beta }_{max}=0$. Figure 6 shows the relationship between ${\beta }_{min}$, ${\beta }_{max}$, and $v_x$, $\delta$ on a road surface with an adhesion coefficient $\mu =0.85$. A lookup table is established for subsequent use.

Through (14), the evaluation indices $I_{\beta }$ and $I_{\gamma }$, which represent the distances of the coordinates ($\dot{\beta }$, $\beta$) on the phase plane trajectory to their respective boundaries for $\gamma$ and $\beta$, are obtained. When the evaluation index falls within the range [0, 0.8], it indicates the stable region, as shown by the light-yellow area in Figure 7. When the evaluation index is in the range [0.8, 1], it indicates the critical region, as shown by the light red area in Figure 7. If the evaluation index exceeds 1, it indicates the unstable region. Substituting into (15), the maximum value of $I_{\beta }$ and $I_{\gamma }$ is used as the input coefficient $u$ for determining the stability control weight coefficient $W$ in subsequent discussions. To ensure a smooth transition of the stability control weight coefficient $W$, a smooth step function described by (16) is applied, as shown in Figure 8, determining the relationship between the input coefficient u and the vehicle stability control weight coefficient $W$.

$$\left\{\begin{array}{c}{I}_{\beta }=1-sign\left[\left({\beta }_{max}-\beta \right)\left(\beta -{\beta }_{min}\right)\right]\times {\displaystyle \frac{min\left(\left|\left({\beta }_{max}-\beta \right)\right|,\left|\left(\beta -{\beta }_{min}\right)\right|\right)}{0.5\times \left({\beta }_{max}-{\beta }_{min}\right)}}\\ I_{\gamma }=1-sign\left[\left({\gamma }_{max}-\gamma \right)\left(\gamma -{\gamma }_{min}\right)\right]\times {\displaystyle \frac{min\left(\left|\left({\gamma }_{max}-\gamma \right)\right|,\left|\left(\gamma -{\gamma }_{min}\right)\right|\right)}{0.5\times \left({\gamma }_{max}-{\gamma }_{min}\right)}} \end{array}\right.$$

(14)

$$u=\max \left(I_{\beta },I_{\gamma }\right)$$

(15)

$$W=\left\{\begin{array}{cc} 0,\hfill & \hfill 0\le u<0.8\\ \begin{array}{c}{\displaystyle \frac{1}{2}}\left(1-cos\left(\pi {\displaystyle \frac{u-0.8}{1-0.8}}\right)\right),\end{array}& \hfill 0.8\le u<1\\ 1,\hfill & \hfill u > 1\hfill \end{array}\right.$$

(16)

4. Controller Design

Figure 9 shows a block diagram of the overall system. The controller consists of an upper-level controller and a lower-level controller. The upper-level controller includes four components: the handling assistance controller, the stability controller, the speed-following controller, and the stability criterion method. The lower-level controller mainly handles the control allocation. The upper-level controller, based on the driver’s inputs, uses the stability criterion method to calculate the proportional weight of the additional yaw moment generated by the handling assistance control and stability control, producing the necessary yaw moment to enhance vehicle handling or maintain vehicle stability. The lower-level controller distributes the additional yaw moment and the torque required to maintain vehicle speed to the individual wheels.

4.1. Handling Assistance Controller

To improve the handling performance of the vehicle during stable driving, this paper proposes a handling assistance controller. The block diagram is shown in Figure 10. Based on the driver’s input, such as front wheel steering angle $\delta$ and target vehicle speed $v_t$, the system generates additional yaw moment $M_{z,ff}$ applied to the vehicle through Feed-Forward Control, making the vehicle sideslip angle $\beta$ approach the ideal reference sideslip angle ${\beta }_{d,h}$. Considering that the vehicle may deviate from the ideal value due to factors such as road conditions or wind during actual operation, the handling assistance controller adopts a Linear Quadratic Regulator (LQR). It calculates the required additional yaw moment $M_{z,fb,hand}$ to improve vehicle handling by using the difference between the ideal reference vehicle yaw rate $r_{d,h}$, the reference sideslip angle ${\beta }_{d,h}$, and the actual yaw rate $\gamma$ and actual sideslip angle $\beta$ as input values.

Taking the Laplace transform of (7), then obtain the following (17):

$$\left\{\begin{array}{c}\beta \left(s\right)={\displaystyle \frac{\left(g_{1}s-g_1{a}_{22}+g_2{a}_{12}\right){\delta }_f\left(s\right)+b_2{a}_{12}{M}_z\left(s\right)}{\left(s-a_{11}\right)\left(s-a_{22}\right)-a_{12}{a}_{21}}} \\ \gamma \left(s\right)={\displaystyle \frac{\left(g_{2}s+g_1{a}_{21}-g_2{a}_{11}\right){\delta }_f\left(s\right)+b_2\left(s-a_{11}\right)M_z\left(s\right)}{\left(s-a_{11}\right)\left(s-a_{22}\right)-a_{12}{a}_{21}}} \end{array}\right.$$

(17)

To make the vehicle sideslip angle $\beta$ approach the ideal value, a feed-forward control is designed. The relationship between the additional control yaw moment $M_{z,ff}$ and $\delta$ is given by (18):

$$M_{z,ff}=G_{ff}\delta \left(s\right)$$

(18)

Substituting (18) into the first equation of (17) we can obtain the (19):

$${\displaystyle \frac{{\beta }_0}{{\delta }_0}}={\displaystyle \frac{-g_1{a}_{22}+g_2{a}_{12}+b_2{a}_{12}{G}_{ff}}{a_{11}{a}_{22}-a_{12}{a}_{21}}}$$

(19)

When the vehicle is in a steady state, ${\beta }_0=0$, then $G_{ff}$ can be expressed as:

$$G_{ff}={\displaystyle \frac{g_1{a}_{22}-g_2{a}_{12}}{b_2{a}_{12}}}$$

(20)

Thus, $M_{z,ff}$ is showed as (21):

$$M_{z,ff}={\displaystyle \frac{g_1{a}_{22}-g_2{a}_{12}}{b_2{a}_{12}}}\delta$$

(21)

In the feedback control of the handling assistance control system, a Linear Quadratic Regulator (LQR) is used. When $u=0$ in (7), the ideal reference values ${\beta }_{d,h}$ and $r_{d,h}$ in the handling assistance control system can be expressed as:

$$\left[\begin{array}{c}{\dot{\beta }}_{d,h}\\ {\dot{\gamma }}_{d,h}\end{array}\right]=\left[\begin{array}{cc}{a}_{11}& \hfill a_{12}\\ a_{21}& \hfill a_{22}\end{array}\right]\left[\begin{array}{c}{\beta }_{d,h}\\ r_{d,h}\end{array}\right]+\left[\begin{array}{c}{g}_1\\ g_2\end{array}\right]\delta$$

(22)

Subtracting (22) from (7), then we arrive at the new Equation:

$$\dot{\mathit{x}}=\mathit{A}\mathit{x}+\mathit{B}u$$

(23)

where $\mathit{x}=\left[\begin{array}{c}{\beta }_{d,h}-\beta \\ r_{d,h}-r\end{array}\right]$, $\mathit{A}=\left[\begin{array}{cc}{a}_{11}& a_{12}\\ a_{21}& a_{22}\end{array}\right]$, $\mathit{B}=\left[\begin{array}{c}{b}_1\\ b_2\end{array}\right]$ and $u=M_z$.

The objective function of the controller is expressed as (24)

$$J={\displaystyle \frac{1}{2}}{{\displaystyle \int }}_0^{\infty }\left[x^T{Q}_{hand}x+u^{T}Ru\right]dt$$

(24)

where $Q_{hand}=\left[\begin{array}{cc}{Q}_1& 0\\ 0& Q_2\end{array}\right]$ and $R$ are the weighting matrices for the state deviations and the input error. Since the handling assistance control system primarily follows the yaw rate, $Q_2$ is much greater than $Q_1$. By solving the Algebraic Riccati Equation (25), the gain matrix ${\mathit{K}}_h$ for the handling assistance control is obtained. Substituting into (27), the required additional yaw moment $M_{z,fb,hand}$ is calculated.

$$A^{T}S+SA-S\mathit{B}{R}^{-1}{\mathit{B}}^{T}S+Q=0$$

(25)

$${\mathit{K}}_h=R^{-1}{\mathit{B}}^{T}S$$

(26)

$$M_{z,fb,hand}={\mathit{K}}_{h}x$$

(27)

In summary, combining (21) and (27), the additional yaw moment $M_{z,hand}$ in the handling assistance control system is calculated as:

$$M_{z,hand}=M_{z,ff}+M_{z,fb,hand}$$

(28)

4.2. Stability Controller

The block diagram of the stability controller from this paper is shown in Figure 11. The reference model adopts the handling assistance control reference model. To prevent the loss of control, the target value for the reference vehicle sideslip angle ${\beta }_{d,s}$ is set to 0, while maintaining a certain level of steering capability, with the yaw rate $r_{d,s}=r_{d,h}$. The LQR controller is used to calculate the additional yaw moment $M_{z,stab}$. To emphasize limiting the vehicle sideslip angle, $Q_{stab}$ must satisfy the condition in (29) and $Q_3\gg Q_4$. Solving the Algebraic Riccati Equation yields the gain matrix ${\mathit{K}}_s$, and finally, the additional yaw moment from the stability controller $M_{z,stab}$ is obtained.

$$Q_{stab}=\left[\begin{array}{cc}{Q}_3& 0\\ 0& Q_4\end{array}\right]$$

(29)

$$M_{z,stab}={\mathit{K}}_{s}x$$

(30)

4.3. Speed-Following Controller

This speed-following controller uses PID control to obtain the required speed-following torque $T_{vx}$, achieving control over the target vehicle speed. The specific flowchart is shown in Figure 12.

4.4. Control Allocation

Combining Section 3.2, Section 4.1 and Section 4.2, the additional yaw moment $M_{z,total}$, which improves vehicle handling and enhances stability, is expressed in (31):

$$M_{z,total}=\left(1-W\right)·M_{z,hand}+W·M_{z,stab}$$

(31)

In order to effectively distribute the additional yaw moment generated by the upper-level controller to the four wheels, the relationship shown in (32) is derived based on the analysis of forces acting on the wheels. To prevent tire forces from reaching saturation and causing vehicle instability, this paper uses Control Allocation and applies Quadratic Programming (QP) to solve for the optimal torque values distributed to each wheel.

$$\left\{\begin{array}{c}{M}_{z,total}={\displaystyle \frac{T_{fl}}{R}}\left(-{\displaystyle \frac{t_f}{2}}cos\delta +l_f\sin \delta \right)+{\displaystyle \frac{T_{fr}}{R}}\left({\displaystyle \frac{t_f}{2}}cos\delta +l_{f}sin\delta \right)-{\displaystyle \frac{T_{rl}{t}_r}{2\mathrm{R}}}+{\displaystyle \frac{T_{rr}{t}_r}{2\mathrm{R}}}\\ T_{vx}=T_{fl}cos\delta +T_{fr}cos\delta +T_{rl}+T_{rr} \end{array}\right.$$

(32)

where $T_{ij}$ is the torque applied to the corresponding wheel of the vehicle; $R$ is the radius of the wheel.

From the tire adhesion limit circle, we obtain the (33). Since (33) is a nonlinear inequality, the computational complexity is high, and the possibility of no solution may arise during the solving process. In this paper, the tire adhesion limit circle is simplified into a linear octagon, as shown in Figure 13. The octagon is inscribed within the adhesion limit circle, converting the nonlinear inequality in (33) into a linear inequality constraint problem, as shown in (34).

$$\sqrt{{\left(F_{xij}\right)}^2+{\left(F_{yij}\right)}^2}\le {\mu }_{ij}\ast F_{zij}$$

(33)

From Figure 13, we can obtain $R_F=\mu ·F_{zij}$, and then the constraint in (33) can be transformed into (34):

$$\left\{\begin{array}{c}-\sqrt{2}\times \cos \left({22.5}^{°}\right)\times \mu ·F_{zij}-F_{yij}\le F_{xij}\le \sqrt{2}\times \cos \left({22.5}^{°}\right)\times \mu ·F_{zij}-F_{yij}\\ -\cos \left({22.5}^{°}\right)\times \mu ·F_{zij}\le F_{xij}\le \cos \left({22.5}^{°}\right)\times \mu ·F_{zij}\end{array}\right.$$

(34)

The objective of torque distribution is to minimize the tire-road adhesion utilization of the four wheels. Assuming the adhesion coefficient for all four wheels is $\mu$, the objective function can be expressed as shown in (35):

$$J=min\sum {\displaystyle \frac{1}{{\left(\mu \ast F_{zij}\right)}^2}}\ast \left[{\left(F_{xij}\right)}^2+{\left(F_{yij}\right)}^2\right]$$

(35)

Due to the coupling effects between the tire’s longitudinal force, lateral force, and vertical force, to reduce computational complexity, this paper considers only the optimization of the tire’s longitudinal force. Therefore, (35) is transformed into (36):

$$J=min\sum {\displaystyle \frac{1}{{\left(\mu \ast F_{zij}\right)}^2}}\ast \left[{\left(F_{xij}\right)}^2\right]$$

(36)

From (32), (34), and (36), we can derive the quadratic objective function of the optimization problem as (37):

$$\mathit{min}{\displaystyle \frac{1}{2}}{T}_x^T H T_x, subject to\left\{\begin{array}{c}{\mathit{A}}_{\mathit{e}\mathit{q}}·{\mathit{T}}_x={\mathit{b}}_{\mathit{e}\mathit{q}}\\ \mathit{A}·{\mathit{T}}_x\le \mathit{b}\end{array}\right.$$

(37)

where ${\mathit{A}}_{\mathit{e}\mathit{q}}$ is the control matrix that describes the relationship between the command vector ${\mathit{b}}_{\mathit{e}\mathit{q}}$ and the input vector ${\mathit{T}}_x$ and Matrix $\mathit{A}$ describes the linear inequality constraint relationship between the input vector ${\mathit{T}}_x$ and the command vector $\mathit{b}$.

The parameters for the optimization problem are:

$${\mathit{A}}_{\mathit{e}\mathit{q}}=\left[\begin{array}{cccc}{\displaystyle \frac{1}{R}}\left(-{\displaystyle \frac{t_f}{2}}cos\delta +l_f\sin \delta \right)\hfill & \hfill {\displaystyle \frac{1}{R}}\left({\displaystyle \frac{t_f}{2}}cos\delta +l_{f}sin\delta \right)\hfill & \hfill -{\displaystyle \frac{t_r}{2\mathrm{R}}}& \hfill {\displaystyle \frac{t_r}{2\mathrm{R}}}\\ cos\delta & cos\delta & 1& 1\end{array}\right], {\mathit{b}}_{\mathit{e}\mathit{q}}=\left[\begin{array}{c}{M}_{z,total}\\ T_{vx}\end{array}\right]$$

$${\mathit{T}}_x={\left[\begin{array}{cccc}{T}_{fl}& T_{fr}& T_{rl}& T_{rr}\end{array}\right]}^T, \mathit{A}=\left[\begin{array}{c}\mathit{A}1\\ \mathit{A}2\end{array}\right], \mathit{b}=\left[\begin{array}{c}\mathit{b}1\\ \mathit{b}2\end{array}\right]$$

$$\mathit{A}1={\displaystyle \frac{1}{R}}\times \left[\begin{array}{c}\begin{array}{cccc}1& 0& 0& 0\\ -1& 0& 0& 0\end{array}\\ \begin{array}{cccc}0& 1& 0& 0\\ 0& -1& 0& 0\end{array}\\ \begin{array}{cccc}1& 0& 0& 0\\ -1& 0& 0& 0\end{array}\\ \begin{array}{cccc}1& 0& 0& 0\\ -1& 0& 0& 0\end{array}\end{array}\right], \mathit{A}2={\displaystyle \frac{1}{R}}\times \left[\begin{array}{c}\begin{array}{cccc}1& 0& 0& 0\\ -1& 0& 0& 0\end{array}\\ \begin{array}{cccc}0& 1& 0& 0\\ 0& -1& 0& 0\end{array}\\ \begin{array}{cccc}1& 0& 0& 0\\ -1& 0& 0& 0\end{array}\\ \begin{array}{cccc}1& 0& 0& 0\\ -1& 0& 0& 0\end{array}\end{array}\right], \mathit{b}1={\displaystyle \frac{1}{R}}\times \left[\begin{array}{c}\cos \left(22.5\right)\mu ·F_{zfl}\\ \cos \left(22.5\right)\mu ·F_{zfl}\\ \cos \left(22.5\right)\mu ·F_{zfr}\\ \cos \left(22.5\right)\mu ·F_{zfr}\\ \cos \left(22.5\right)\mu ·F_{zrl}\\ \cos \left(22.5\right)\mu ·F_{zrl}\\ \cos \left(22.5\right)\mu ·F_{zrr}\\ \cos \left(22.5\right)\mu ·F_{zrr}\end{array}\right]$$

$$\mathit{b}2={\displaystyle \frac{1}{R}}\times \left[\begin{array}{c}\sqrt{2}\cos \left(22.5\right)\mu ·F_{zfl}-F_{yfl}\\ \sqrt{2}\cos \left(22.5\right)\mu ·F_{zfl}-F_{yfl}\\ \sqrt{2}\cos \left(22.5\right)\mu ·F_{zfr}-F_{yfr}\\ \sqrt{2}\cos \left(22.5\right)\mu ·F_{zfr}-F_{yfr}\\ \sqrt{2}\cos \left(22.5\right)\mu ·F_{zrl}-F_{yrl}\\ \sqrt{2}\cos \left(22.5\right)\mu ·F_{zrl}-F_{yrl}\\ \sqrt{2}\cos \left(22.5\right)\mu ·F_{zrr}-F_{yrr}\\ \sqrt{2}\cos \left(22.5\right)\mu ·F_{zrr}-F_{yrr}\end{array}\right], \mathit{H}={\displaystyle \frac{1}{{\left(\mu R\right)}^2}}\left[\begin{array}{cccc}{\displaystyle \frac{1}{{\left(F_{zfl}\right)}^2}}& 0& 0& 0\\ 0& {\displaystyle \frac{1}{{\left(F_{zlr}\right)}^2}}& 0& 0\\ 0& 0& {\displaystyle \frac{1}{{\left(F_{zrl}\right)}^2}}& 0\\ 0& 0& 0& {\displaystyle \frac{1}{{\left(F_{zrr}\right)}^2}}\end{array}\right]$$

5. Results

Based on different vehicle stability criterion methods, including the double-line method, the diamond method, the curved boundary method, and the “normalization” method, four test groups were designed using MATLAB and Simulink (version 2023b) to build the controllers: Case A, Case B, Case C, and Case D. These controllers were tested in a Carsim–Simulink co-environment. The evaluation of the controllers was conducted through open-loop control using the Sine with Dwell steering test and closed-loop control using the DLC test.

The ramp input steering test can intuitively demonstrate the linearity of the vehicle’s response to steering wheel input. With slight modifications to the ‘Slowly Increasing Steer Procedure’ from Euro NCAP, the test was conducted under the following conditions: road surface adhesion coefficient μ = 0.85, vehicle speed $v_x=80 \mathrm{km}/\mathrm{h}$, and steering wheel angle ${\delta }_{hand}$ input at 13.5 $\mathrm{deg}/\mathrm{s}$. The results are shown in Figure 14. When $a_y=0.3 \mathrm{g}$, all four test vehicles were in a stable state, so only the handling assistance control was activated, with the corresponding ${\delta }_{hand}$ for point $\mathrm{A}$ being 22.92 $\mathrm{deg}$, which prepares for the subsequent tests.

5.1. Sine with Dwell Steering Test

Sine with Dwell steering is an open-loop testing method used to assess the dynamic performance of a vehicle under extreme conditions. According to FMVSS 126, which serves as the standard for ESC, there are three criteria for passing this test. First, the vehicle’s lateral displacement must exceed 1.83 m at 1.07 s after the sine wave begins. Second, at 1 s and 1.75 s after the sine wave input ends, the vehicle’s yaw rate ${\gamma }_{T_1}$ and ${\gamma }_{T_2}$ must be less than 35% and 25% of the peak value, respectively. The final criterion is that the amplitude of the input sine wave at the end of the test must be greater than 270 $\mathrm{deg}$. From the ramp input test, the value of A was determined to be 22.92 $\mathrm{deg}$, and here the sine wave input is conducted with an amplitude of 12 A, equal to 275 $\mathrm{deg}$, as shown in Figure 15.

The results are shown in Figure 16 and

Table 1. Sine with Dwell steering test results.

$Case$	${\mathit{Y}}_{1.07}\left(\mathit{m}\right)$	${\mathit{r}}_{\mathit{p}\mathit{e}\mathit{a}\mathit{k}}\left(\mathit{d}\mathit{e}\mathit{g}/\mathit{s}\right)$	${\mathit{\gamma }}_{{\mathit{T}}_1}\left(\mathit{d}\mathit{e}\mathit{g}/\mathit{s}\right)$	${\mathit{P}}_{{\mathit{\gamma }}_{\mathit{T}1}}$	${\mathit{\gamma }}_{{\mathit{T}}_2}\left(\mathit{d}\mathit{e}\mathit{g}/\mathit{s}\right)$	${\mathit{P}}_{{\mathit{\gamma }}_{\mathit{T}2}}$	Result
A	$3.41$	51.54	6.35	12.32%	0.24	0.47%	Pass
B	$3.42$	50.90	6.08	11.94%	$0.21$	0.41%	Pass
C	$3.31$	45.92	0.41	0.89%	0.00	0%	Pass
D	$3.23$	42.31	0.07	0.16%	0.00	0%	Pass

. At 1.07 s after the sine wave began, the lateral displacements of the four test groups were 3.41 m, 3.42 m, 3.31 m, and 3.23 m, respectively, with all exceeding 1.83 m. At 1 s after the sine wave ended, the yaw rates ${\gamma }_{T_1}$ for the four test groups were 12.32%, 11.94%, 0.89%, and 0.16% of their respective peak values—all below 35% of the peak value. At 1.75 s after the sine wave ended, all four test vehicles had returned to a stable state, with yaw rates below 25% of the peak value, meeting the requirements of the sine with dwell steering test.

As seen in Figure 17, the phase trajectory of Case D is smaller compared to the other Cases, indicating that Case D has a superior control effect. Specifically, from Figure 18 and the data in

Table 2. Control variable table for the four test cases.

$Case$	${\mathit{\beta }}_{\mathit{m}\mathit{a}\mathit{x}}\left(\mathit{d}\mathit{e}\mathit{g}\right)$	${\mathit{M}}_{\mathit{z}\mathit{m}\mathit{a}\mathit{x}}\left(\mathit{N}\mathit{m}\right)$	${\mathit{T}}_{\mathit{m}\mathit{a}\mathit{x}}\left(\mathit{N}\mathit{m}\right)$	${\mathit{\lambda }}_{\mathit{m}\mathit{a}\mathit{x}}$
A	$14.78$	1808.84	322.80	14.94%
B	$14.62$	1631.65	296.56	14.76%
C	$8.74$	1612.86	292.45	11.06%
D	$7.19$	1536.19	275.24	9.74%

, it can be observed that in terms of the maximum vehicle sideslip angle ${\beta }_{max}$, Case D reduced it from 14.78 $\mathrm{deg}$ (in Case A, B, and C) to 7.19 $\mathrm{deg}$, a decrease of 51.35%. The smaller sideslip angle provides the driver with a better driving experience. In terms of the additional torque required to maintain stable vehicle driving $M_{zmax}$, it decreased from 1808.84 $\mathrm{Nm}$ to 1536.19 $\mathrm{Nm}$—a reduction of 15.07%. The maximum torque distributed to the four wheels $T_{max}$ decreased from 322.80 $\mathrm{Nm}$ to 275.24 $\mathrm{Nm}$—a reduction of 14.73%. Lastly, the vehicle’s maximum slip ratio ${\lambda }_{max}$ decreased from 14.94% to 9.74%, improving the tire’s longitudinal adhesion capability and significantly enhancing vehicle stability.

Case D demonstrates better vehicle stability control compared to Cases A, B, and C. The reason, as illustrated in Figure 17 and Figure 18a, is that Case A’s ‘double-line method’ has an open stable region, where a misjudgment of stability occurs when the vehicle is actually in an unstable state due to a large $\beta$, leading to a delay in the activation or deactivation of the stability control system. Although Case B’s ‘diamond method’ has a closed region, it also fails to account for the $\gamma$ and $\beta$ limits during actual vehicle operation, resulting in an overly large stable region, which similarly delays the intervention or deactivation of the stability control system. While Case C provides better stability control than Cases A and B, it does not consider the influence of a large front wheel steering angle $\delta$ (which can reach up to 16.26 $\mathrm{deg}$ in this test) on the position of the stability center and the size of the stable region. This leads to a delay in recognizing vehicle instability at 1.13 s compared to Case D by 0.32 s, and at 1.81 s, when the steering wheel rapidly crosses 0 $\mathrm{deg}$, it misjudges the vehicle as stable, requiring a larger additional torque $M_{zmax}$ to maintain vehicle stability in subsequent moments. Case D comprehensively considers the effects of road adhesion coefficient $\mu$, vehicle speed $v_x$, and front wheel angle $\delta$ on the size and position of the stability region in the phase plane. This approach avoids potential misjudgments that could arise when assessing the stability of dynamic points ($\beta$, $\dot{\beta }$) within the phase plane by using a fixed region partitioning method. Furthermore, a smooth step function is employed to ensure a seamless transition for the stability control weighting factor $W$.

5.2. Double Lane Change Test

The DLC test is a closed-loop test method. According to the international standard (ISO 3888-1 [19]), the test route is shown in Figure 19. Using the driver model provided by Carsim, four test cases were conducted under high vehicle speed conditions ($v_x=80 km/h$) and a low road adhesion coefficient ($\mu = 0.3$).

The test results are shown in Figure 20. From the driving trajectory in Figure 20a, it can be seen that all four test vehicles follow the predetermined target route steadily, with Case D switching lanes more promptly and closely following the target path. From the lateral acceleration $a_y$ dynamic graph in Figure 20b, it is evident that all test vehicles reach 2.82 m/s². Case D exhibits more sensitive acceleration compared to the other three cases, with its acceleration increasing from −2.82 m/s² to 2.82 m/s² 0.1 s earlier at 4.55 s, indicating better handling performance for Case D. From Figure 20, Figure 20, and

Table 3. Control variable table for the four test cases in DLC test with $v_x=80 km/h$ and $\mu = 0.3$.

Case	${\mathit{\beta }}_{\mathit{m}\mathit{a}\mathit{x}}\left(\mathit{d}\mathit{e}\mathit{g}\right)$	${\mathit{r}}_{\mathit{m}\mathit{a}\mathit{x}}\left(\mathit{d}\mathit{e}\mathit{g}/\mathit{s}\right)$	${\mathit{M}}_{\mathit{z}\mathit{m}\mathit{a}\mathit{x}}\left(\mathit{N}\mathit{m}\right)$	${\mathit{T}}_{\mathit{m}\mathit{a}\mathit{x}}\left(\mathit{N}\mathit{m}\right)$	${\mathit{\lambda }}_{\mathit{m}\mathit{a}\mathit{x}}$
A	7.86	26.34	1236.92	208.02	6.16%
B	7.46	25.19	1236.92	208.22	6.08%
C	5.70	22.51	1029.65	181.44	3.96%
D	2.17	12.92	693.73	138.40	1.30%

, it is clear that Case D has a significant advantage over the other three cases. The maximum sideslip angle ${\beta }_{max}$ and maximum yaw rate ${\gamma }_{max}$ are reduced to 2.17 $\mathrm{deg}$ and 12.92 $\mathrm{deg}∕\mathrm{s}$, representing decreases of 72.39% and 50.95%, respectively, compared to the other cases. In terms of the yaw moment $M_{zmax}$ required to maintain stable driving, Case D reduced $M_{zmax}$ from 1236.92 $\mathrm{Nm}$ to 693.73 $\mathrm{Nm}$—a decrease of 43.91%. The maximum torque $T_{max}$ distributed to the wheels was reduced from 208.22 $\mathrm{Nm}$ to 138.40 $\mathrm{Nm}$, a reduction of 33.53%, while the maximum wheel slip rate decreased from 6.16% to 1.30% as the distributed torque decreased, significantly improving the vehicle’s stability during the rapid lane change on the low-adhesion surface.

Figure 21a and Figure 22 clearly explain why Case D performed the best among the four test cases. In Cases A and B, the stability boundaries were inaccurately determined during vehicle stability assessment. Additionally, the introduction of $\mu$ from Reference [20] caused a smooth transition between stability control and handling assistance control, but when $\mu$ is small, the transition zone (critical region) becomes too large, resulting in delays and inaccuracies in the switching and intervention of the handling assistance and stability control systems. Although Case C considered the limitations of $\gamma$ and $\beta$, it did not account for the influence of δ on the position of the stability center and the size of the stable region. This led to a 0.26 s delay in identifying vehicle instability compared to Case D at 1.06 s. At 3.13 s, after completing a single lane change, a misjudgment of vehicle stability occurred, requiring a larger additional torque $M_{zmax}$ to maintain vehicle stability.

6. Conclusions

The following research results were achieved in this study:

A hierarchical lateral stability control system was designed. The upper-level controller automatically selects the handling assistance controller or stability controller based on driver input and vehicle operating conditions through a stability assessment method. It calculates the yaw moment needed to maintain vehicle following of the ideal reference model and transmits it to the lower-level controller. Additionally, to prevent sudden acceleration or deceleration during driving, a speed-following controller uses a PID controller to calculate the torque required to maintain speed and transmits it to the lower-level controller. The lower-level controller then solves the optimal distribution of the additional yaw moment to the four wheels.

To accurately assess vehicle stability, a $\beta -\dot{\beta }$ phase plane was established using a four-wheel vehicle dynamics model and a nonlinear Magic Formula Tire model, which takes into account the effects of vehicle body roll. A “normalization” method was proposed to assess vehicle stability. Based on different stability assessment methods, four sets of experimental groups were designed using the double-line method, diamond method, curve boundary method, and “normalization” method. An Open-loop Sine with Dwell steering test and a closed-loop DLC test were conducted. The results show that the proposed “normalization” method is more sensitive and accurate in assessing vehicle stability. Moreover, the results indicate that the “normalization” method significantly reduces the additional torque required for distribution to the wheels when improving vehicle handling and enhancing stability.

The proposed control system requires the accurate acquisition of relevant parameters reflecting vehicle motion, such as the vehicle’s sideslip angle $\beta$ and speed $v_x$. Accurately obtaining these parameters is a prerequisite for the proper functioning of the control system. In this test, the parameters were directly obtained through vehicle simulation software. However, it is difficult to accurately measure these parameters using traditional sensors on real vehicles. After obtaining the kinematic parameters via measurement or the state-estimation method, experiments should be conducted in the following research to test the proposed control system. Additionally, during vehicle acceleration/deceleration and cornering, the vertical load on each of the four wheels varies due to load transfer (longitudinal transfer during acceleration/deceleration and lateral transfer during cornering). This affects the forces on the wheels, subsequently impacting the vehicle’s dynamic characteristics. Therefore, it is essential in future research to continuously and accurately estimate the vertical load on each wheel.