Trajectory Tracking Coordinated Control of 4WID-4WIS Electric Vehicle Considering Energy Consumption Economy Based on Pose Sensors

Abstract

In order to improve the stability and economy of 4WID-4WIS (four-wheel independent drive—four-wheel independent steering) electric vehicles in trajectory tracking, this paper proposes a trajectory tracking coordinated control strategy considering energy consumption economy. First, a hierarchical chassis coordinated control architecture is designed, which includes target planning layer, and coordinated control layer. Then, the trajectory tracking control is decoupled based on the decentralized control structure. Expert PID and Model Predictive Control (MPC) are employed to realize longitudinal velocity tracking and lateral path tracking, respectively, which calculate generalized forces and moments. In addition, with the objective of optimal overall efficiency, the optimal torque distribution for each wheel is achieved using the Mutant Particle Swarm Optimization (MPSO) algorithm. Additionally, the modified Ackermann theory is used to distribute wheel angles. Finally, the control strategy is simulated and verified using Simulink. Comparing the control results of the average distribution strategy and the wheel load distribution strategy, it can be concluded that the proposed coordinated control not only provides good trajectory tracking but also greatly improves the overall efficiency of the motor operating points, which enhances the energy economy and realizes the multi-objective coordinated control of the chassis.

1. Introduction

4WID-4WIS EV is a novel electric vehicle that reduces mechanical components, such as differentials and half shafts. It also integrates four-wheel steering technology based on a distributed four-wheel drive system. This allows four wheels to be driven and steered independently, increasing the controllable degree of freedom. In a word, 4WID-4WIS EV has unique advantages in vehicle dynamics control [1,2].

To fully exploit the control potential of 4WID-4WIS EV and improve the overall performance of the vehicle, coordinated multi-objective control of the chassis has become a current research focus [3,4,5]. In particular, stability and economy are two important performance characteristics. It is difficult to improve one performance by controlling only one performance, and in extreme cases, the effect deteriorates. Therefore, there is a need for coordinated control of the two performances to make the vehicle stable and economical at the same time, which has been researched and studied by several scholars [6]. In [7], a 4WID EV torque coordination control strategy is developed, which considers both stability and economy and uses MPC to calculate the generalized force. The demand torque of each wheel is determined by an online solution using the control distribution error, tire utilization rate, and power consumption of the drive system as objective functions. In [8], a multi-objective online optimization of energy management strategy for 4WID EV is proposed. It considers the efficiency of the drive system, tire slip energy consumption, wheel torque fluctuation, yaw rate tracking error, etc. The weights of each element are dynamically adjusted using the fuzzy control method, and the effectiveness of the strategy is verified by simulations. In [9], a cooperative game-based actuator fault-tolerant control strategy for a 4WID EV is proposed. This strategy minimizes tire energy dissipation to ensure economic efficiency by designing a two-dimensional game controller that simultaneously satisfies the generalized forces required for vehicle stability.

Due to the great potential of intelligent driving technology to enhance vehicle safety [10,11], improve traffic efficiency, and reduce energy consumption [12], the research on trajectory tracking control of 4WID-4WIS electric vehicles has received increasing attention in the automotive industry, and numerous control methods have been developed [13,14,15]. In [16], LTV-MPC (Linear Time-Varying Model Predictive Control) based on DYC (Direct Yaw Control) is used to realize velocity tracking and trajectory tracking, which improves stability during trajectory tracking. In [17], a robust path-tracking controller for 4WID-4WIS agricultural robotic vehicles is designed by combining the backstepping SMC (Sliding Mode Control) theory to improve the robustness of trajectory tracking. In [18], a four-wheel steering controller and a velocity tracking controller for 4WID-4WIS EV are designed based on the SMC algorithm to improve the accuracy of trajectory tracking. In [19], the feedback linear quadratic regulator (LQR) controller is used to realize emergency avoidance under 4WS high-speed conditions.

There are certain shortcomings in the aforementioned research which can be summarized as follows. Most studies on chassis coordinated control is related to traditional 4WID EV, while there is relatively little research on 4WID-4WIS EV. In addition, most existing research establishes multi-objective optimization problems and solves them online to obtain the control variables. This not only makes the objective function too complex and difficult to solve but also has a negative impact on the real-time performance of the controller. Last but not least, many studies on the trajectory tracking control of 4WID-4WIS EV mostly only consider the trajectory tracking performance. Some studies also consider stability, but there is almost no research that simultaneously considers trajectory tracking performance, stability and economy.

To improve the stability and economy of 4WID-4WIS EV in trajectory tracking, the contributions of this paper are summarized as follows:

• A hierarchical architecture for chassis coordinated control is designed including target planning layer, and coordinated control layer.
• The MPC path tracking controller is constructed, which takes into account stability constraints, such as yaw rate and tire slip angle.
• The MPSO algorithm is also used to create a mapping of the distribution coefficients by optimizing the torque distribution between the front and rear wheels offline, thereby reducing the computational cost.

The paper is organized as follows: Section 2 proposes a hierarchical chassis coordination control. Section 3 highlights the simulation validation and compares the simulation results of different torque distribution strategies. Section 4 sums up the conclusions.

2. Trajectory Tracking Coordinated Control

2.1. Vehicle Model

4WID-4WIS EV uses in-wheel motors and steering motors to replace the traditional transmission and steering mechanism, giving it more control freedom and thus enabling more advanced control methods to improve the overall performance of the vehicle. Considering the strong coupling between the vehicle subsystems and the strong nonlinearity of the tires and motor actuators, it is of great significance to establish a reasonable and accurate vehicle model to study the control strategy of 4WID-4WIS EV.

2.1.1. Vehicle Dynamic Model

The vehicle is a complex multi-body dynamic system, and the vehicle dynamic model is the basis for dynamic analysis, active control, function realization, and parameter optimization. According to the research objective and concern, many scholars have proposed various linear and nonlinear vehicle models with different degrees of complexity [20,21]. This work focuses on the longitudinal and lateral coupling control of the 4WID-4WIS EV. The coupling effects between longitudinal motion and lateral motion of the vehicle include three categories: kinematic coupling, tire force coupling, and load transfer coupling. Kinematic coupling refers to the longitudinal motion being influenced by the longitudinal component of the steering wheel lateral deflection force due to the presence of the wheel steering angle. Additionally, the lateral motion is also influenced by the longitudinal velocity. Tire force coupling is the interaction between tire lateral and longitudinal forces, the combined force of which is constrained by the friction ellipse. The vertical load redistribution is caused by longitudinal acceleration or lateral acceleration, which in turn influences the longitudinal and lateral dynamics. Therefore, a vehicle dynamic model with seven degrees of freedom is adopted, i.e., longitudinal motion, lateral motion, yaw motion, and the rotation of four wheels.

The vehicle is assumed to be left-right symmetric about the center plane and driven on a flat horizontal road, ignoring the vertical motion of the body. At the same time, it is assumed that the suspension system is a rigid structure and the body pitch and roll motion is neglected. After conducting the above modifications, the simplified dynamic model is developed, as shown in Figure 1.

The longitudinal motion equation is:

$$m\left({\dot{v}}_x-v_y{\omega }_r\right)=F_{xfl}+F_{xfr}+F_{xrl}+F_{xrr}$$

(1)

The lateral motion equation is:

$$m\left({\dot{v}}_y+v_x{\omega }_r\right)=F_{yfl}+F_{yfr}+F_{yrl}+F_{yrr}$$

(2)

The equation of yaw motion is:

$$I_z{\dot{\omega }}_r=\left(F_{xfr}+F_{xrr}-F_{xfl}-F_{xrl}\right)\frac{d}{2}+\left(F_{yfl}+F_{yfr}\right)l_f-\left(F_{yrl}+F_{yrr}\right)l_r$$

(3)

where,

$$F_{xij}=F_{txij}\cos {\delta }_{ij}-F_{tyij}\sin {\delta }_{ij}$$

(4)

$$F_{yij}=F_{txij}\sin {\delta }_{ij}+F_{tyij}\cos {\delta }_{ij}$$

(5)

$$i\in \left\{f,r\right\}, j\in \left\{l,r\right\}$$

(6)

The equation of wheel motion is:

$$I_w{\dot{\omega }}_{ij}=T_{dij}-F_{txij}{R}_{eff}-T_{bij}$$

(7)

where, $m$ is the vehicle mass, $v_x$ is the longitudinal velocity, $v_y$ is the lateral velocity, ${\omega }_r$ is the yaw rate, $F_{xij}$ and $F_{yij}$ are the longitudinal force and lateral force of each wheel, respectively, $I_z$ is the yaw inertia, $l_f$ is the distance from the center of mass to the front axle, $l_r$ is the distance from the center of mass to the rear axle, $d$ is the wheelbase of the vehicle, $F_{txij}$ and $F_{tyij}$ are the longitudinal and lateral forces of each wheel in the tire coordinate system, respectively, ${\delta }_{ij}$ is the wheel angle, $I_w$ is the wheel moment of inertia, ${\omega }_{ij}$ is the speed of each wheel, $T_{dij}$ is the driving torque of each wheel, $T_{bij}$ is the braking torque of each wheel, and $R_{eff}$ is the effective rolling radius of the wheel.

2.1.2. Tire Model

As the only component that connects the vehicle to the ground, any state of vehicle motion depends on the interaction forces between the tire and the road surface. The tire has a critical impact on the vehicle’s performance, so it is necessary to perform accurate dynamic modeling of tires to better describe the effects of tire forces on vehicle dynamics. The tire model describes the relationship between the tire motion parameters and tire forces, which can be mainly divided into three types: theoretical model, empirical model, and semi-empirical model [22,23]. The theoretical model is formed by studying the deformation mechanism of tires based on the physical essence of tire mechanics, but its complex structure is not conducive to simulation research. The empirical model is based on the experimental tire data, but as a result, it also lacks theoretical support and is poorly scalable. The semi-empirical model combines the advantages of both models with theoretical support and experimental data to obtain key parameters that ensure good accuracy and scalability.

In this study, the most commonly used semi-empirical tire model MF (Magic Formula) is selected. The combination formula of trigonometric functions is used to describe the relationship between tire force, slip rate, and slip angle. This model is suitable for the operating condition with combined longitudinal and lateral forces. The general expression of the MF tire model is as follows:

$$\left\{\begin{array}{l}y=D\sin \left\{C\arctan \left[Bx-E\left(Bx-\arctan \left(Bx\right)\right)\right]\right\}\hfill \\ Y\left(X\right)=y\left(x\right)+S_V\hfill \\ x=X+S_H\hfill \end{array}\right.$$

(8)

where, $B$ is the stiffness factor, $C$ is the shape factor, $D$ is the peak factor, $E$ is the curvature factor, $Y$ is the output variable, that is, longitudinal force, lateral force or aligning torque, $X$ is the input variable, that is, slip rate or slip angle, $S_V$ is the vertical offset, and $S_H$ is the horizontal offset.

As input variables for the MF tire model, slip rate and side slip angle are important parameters for calculating the tire effect, which must be calculated from the vehicle state as shown in Equations (9) and (10), respectively.

$$s_{ij}=\left\{\begin{array}{l}\begin{array}{ll}\frac{{\omega }_{ij}{R}_{eff}-u_{ij}}{{\omega }_{ij}{R}_{eff}}& {\omega }_{ij}{R}_{eff}\ge u_{ij}\left(Drive\right)\end{array}\\ \begin{array}{ll}\frac{{\omega }_{ij}{R}_{eff}-u_{ij}}{u_{ij}}& {\omega }_{ij}{R}_{eff}<u_{ij}\left(Brake\right)\end{array}\end{array}\right.$$

(9)

where, $u_{ij}$ is the speed of the wheel center.

$$\left\{\begin{array}{c}{\alpha }_{fl}={\tan }^{-1}\left(\frac{v_y+{\omega }_r{l}_f}{v_x-{\omega }_r\cdot d/2}\right)-{\delta }_{fl}\\ {\alpha }_{fr}={\tan }^{-1}\left(\frac{v_y+{\omega }_r{l}_f}{v_x+{\omega }_r\cdot d/2}\right)-{\delta }_{fr}\\ {\alpha }_{rl}={\tan }^{-1}\left(\frac{v_y-{\omega }_r{l}_r}{v_x-{\omega }_r\cdot d/2}\right)-{\delta }_{rl}\\ {\alpha }_{rr}={\tan }^{-1}\left(\frac{v_y-{\omega }_r{l}_r}{v_x+{\omega }_r\cdot d/2}\right)-{\delta }_{rr}\end{array}\right.$$

(10)

2.1.3. Motor Model

The 4WID-4WIS EV drive system consists of four in-wheel motors, and independent four-wheel steering is achieved by four steering motors. Therefore, to investigate the longitudinal and lateral coordinated control strategy, the characteristics of the motors are first analyzed, and suitable motor models are established.

• In-wheel motor model

The research focuses on the coordinated control of the vehicle with multiple objectives. For the in-wheel motor model, the charging and discharging efficiency of the motor is only for the steady-state characteristics. The transient characteristics of the motor are not considered. Therefore, the charging and discharging efficiency of the motor can be expressed as in Equation (11).

$${\eta }_{em}=\eta \left(n_{em},T_{em}\right)$$

(11)

where, $n_{em}$ is the motor speed, and $T_{em}$ is the motor torque.

Since the driving state of the vehicle is mainly studied in this paper, the motor efficiency characteristics are shown in Figure 2. The external characteristic curve and the efficiency map are used to characterize the in-wheel motor. The external characteristic curve can be used to determine the maximum torque as a function of wheel speed in real-time. The discharge efficiency of the motor can be determined using a two-dimensional look-up table of motor speed and torque.

2.

Steering Motor Model

Since the steering motor in this study is not focused on the efficiency characteristics, it is simplified from the perspective of the motor actuation effect as a wheel angle tracking model. The delay caused by the steering mechanism is considered as shown in Equation (12).

$${\delta }_{ij\_out}={\delta }_{ij\_req}\cdot \frac{1}{1+\tau s}$$

(12)

where, ${\delta }_{ij\_out}$ is the actual output angle, ${\delta }_{ij\_req}$ is the required output angle, and $\tau$ is the response time constant.

2.2. Vehicle State Acquisition

Before performing chassis coordination control, it is necessary to observe the lateral and longitudinal states of the vehicle, including the vehicle’s position information and the body’s pose information. The position information of the vehicle consists of the coordinates of the vehicle’s lateral and longitudinal axes, velocity and acceleration information, and the pose information of the body includes yaw rate, sides lip angle, pitch angle, roll angle and so on.

The implementation of the control strategy is based on the acquisition of vehicle state information by pose sensors. There have been many studies on how to obtain the pose information of the vehicle. For example, RTK (Real Time Kinematic) and INS (Inertial Navigation System) are combined for positioning. Body combination sensors, including longitudinal acceleration sensors, lateral acceleration sensors and yaw rate sensors are used for pose calculation.

Considering that the signals from pose sensors contain uncertain noise interference, it is necessary to design a reasonable signal filter to obtain accurate pose signals. At the same time, the acquisition of road condition information, such as road adhesion coefficient and road slope, can be realized by designing corresponding state observers based on the information from pose sensors. Since the focus of this paper is on the trajectory tracking coordinated control, this part will not be elaborated in detail.

2.3. Chassis Control Architecture

Based on the vehicle states obtained by pose sensors, the trajectory tracking controller receives the expected path and velocity information sent by the decision planning layer and controls the longitudinal and lateral movement of the vehicle so that the vehicle follows the expected path. In addition, 4WID-4WIS EV has multiple degrees of freedom, and the driving and steering of each wheel are controllable, which provides a basis for realizing multi-objective coordinated control of the chassis. In summary, this research aims to achieve the expected trajectory tracking while improving driving stability and energy economy.

The chassis control architecture designed in this paper is shown in Figure 3, which has a layered structure and is divided into target planning layer and coordinated control layer.

Based on the decentralized control structure, the target planning layer receives the expected path and velocity information and decouples the longitudinal and lateral control targets. In the longitudinal direction, an expert PID control method is used to track the desired velocity and output the generalized longitudinal force. In the lateral direction, the MPC algorithm is used to perform the multi-objective real-time rolling optimization considering the stability constraints, such as yaw rate and tire slip angle. It outputs the generalized yaw moment and the generalized wheel angle.

The coordinated control layer distributes the torque and wheel angle of each wheel according to the generalized force, generalized yaw moment, and generalized wheel angle, which realizes the calculation of the control objective of the motor actuator. Considering the efficiency characteristics of the in-wheel motor, the torque distribution control uses the MPSO algorithm to optimize the torque of the front and rear wheels offline. Then, the two-dimensional diagram of the optimal torque distribution coefficient is generated and real-time table look-up is performed according to the vehicle state to calculate the torque for each wheel. The wheel angle distribution strategy considers the wheel slip angle and distributes the generalized wheel angle based on the modified Ackerman theory to obtain the wheel angles. Finally, the torque and angle of each wheel are output to in-wheel motors and steering motors to realize the closed-loop control of the strategy.

2.4. Target Planning

The target planning layer completes the resolution of the chassis control targets based on the expected trajectory. The expected vehicle trajectory includes information about the expected path and velocity, i.e., the expected coupling of the vehicle’s longitudinal and lateral motion states. Trajectory tracking strategies can be divided into decentralized control and centralized control according to different control structures. Decentralized control [24,25] refers to the decoupling of longitudinal and lateral vehicle motion, i.e., it decomposes the trajectory tracking problem into longitudinal velocity tracking and lateral path tracking problems and designs corresponding control objectives. Centralized control [26,27,28], on the other hand, means that the longitudinal and lateral motion in the trajectory tracking problem is considered uniformly, with a global controller computing the longitudinal and lateral control targets. The centralized trajectory tracking controller considers the system holistically and can better account for the coupling properties between longitudinal and lateral motion control. However, the higher dimensionality and complexity of the system model make it more challenging to design control laws than decentralized control, and it is also more computationally intensive and costly to implement physically. In contrast, the decentralized trajectory tracking strategy is relatively simple and easy to implement in engineering practice. In this work, the trajectory tracking strategy adopts a decentralized structure, decouples the trajectory tracking problem into a longitudinal velocity tracking problem and a lateral path tracking problem, and designs the tracking strategy to calculate the desired generalized force and moment.

2.4.1. Longitudinal Velocity Tracking

The uncertainty of vehicle system parameters, uncertain external disturbances, and nonlinear coupling between systems lead to difficulties in tuning PID control parameters and low robustness of dynamic and steady-state performance of vehicle velocity tracking [29]. In this work, the expert PID algorithm is used in conjunction with the expert system theory to track the longitudinal velocity.

The expert PID algorithm is based on various knowledge of the controlled object and the control laws. If the precise model of the controlled object is not available, the PID parameters are designed using expert experience. The deviations and increments of deviations are used to determine the current state of the controlled system. Rules are designed to adjust the output of the controller to achieve faster and smoother convergence of the controlled system.

Assume that $e(k)=v(k)-v_d(k)$ is the velocity error, and $\mathsf{\Delta }e(k)=e(k)-e(k-1)$ is the velocity error increment. At the same time, based on engineering practice experience, the maximum deviation value is set as $M_{max}$, the middle deviation value is set as $M_{mid}$ and the minimum deviation value is set as $M_{\min }$. According to the setting of deviation, deviation increment and extreme value of deviation, the setting rules are as follows:

(1)

When $\left|e\left(k\right)\right|>M_{max}$, it means that the velocity error is unacceptably large. At this time, the controller should be directly output at full load, that is, $u(k)=F_{\mathrm{xmax}}$.

(2)

When $e\left(k\right)\ast \Delta e\left(k\right)>0, \Delta e\left(k\right)=0$, it means that the velocity deviation is changing in the direction of increasing the absolute value of the deviation, or the deviation is a certain fixed value, then

$$u(k)=u(k-1)+k_1\left\{k_{i}e(k)+k_p\mathsf{\Delta }e(k)+k_d\Delta \mathsf{\Delta }e(k)\right\}$$

(13)

If $\left|e\left(k\right)\right|>M_{mid}$, the velocity deviation is also large, and it may be considered to increase the output gain $k_1$ of the controller.

(3)

When $e(k)\mathsf{\Delta }e(k)<0, \mathsf{\Delta }e(k)\mathsf{\Delta }e(k-1)>0,e(k)=0$, it means that the absolute value of the velocity deviation is changing in the direction of decreasing, or has reached the equilibrium state. Then, the controller output remains unchanged, that is, $u\left(k\right)=u\left(k-1\right)$.

(4)

When $e(k)\mathsf{\Delta }e(k)<0, \mathsf{\Delta }e(k)\mathsf{\Delta }e(k-1)<0$, it means that the velocity deviation is in the limit state, then,

$$u(k)=u(k-1)+k_2{k}_{i}e(k)$$

(14)

If $\left|e\left(k\right)\right|>M_{mid}$, it means that the velocity deviation is also large, and it may be considered to increase the output gain $k_2$ of the controller.

(5)

When $\left|e\left(k\right)\right|<M_{min}$, it means that the absolute value of the velocity deviation is very small. In order to reduce the static error of the system, PI control is implemented:

$$u(k)=u(k-1)+k_{i}e(k)+k_p\mathsf{\Delta }e(k)$$

(15)

To sum up, the generalized longitudinal force is expressed as:

$$F_{xd}(k)=u(k)$$

(16)

2.4.2. Lateral Path Tracking

Intelligent driving faces a complex and changing environment where conditions, such as critical safety constraints and actuator constraints, must be met [30,31,32]. In this paper, MPC is used to realize the lateral path tracking considering the driving stability constraints.

• Predictive model design

The predictive model is the basis of the MPC used to predict the future output of the controlled system. The vehicle is a complex coupled nonlinear dynamic model. More degrees of freedom can improve the modeling accuracy, but the complexity of the model will also increase, making it difficult to meet the requirements for fast model solutions and real-time control. Therefore, dynamic modeling must strike a balance between improving model accuracy and reducing model complexity. Considering the accuracy and real-time requirements for lateral control of vehicles, a nonlinear three degrees-of-freedom vehicle model with lateral and longitudinal coupling is established in this paper to realize the prediction of system output as shown in Figure 4.

Based on the established dynamic model, the force analysis is performed:

$$\left\{\begin{array}{l}\begin{array}{l}m\left({\dot{v}}_y+v_x{\omega }_r\right)=C_f\left(\frac{v_y+l_f{\omega }_r}{v_x}-{\delta }_f\right)+C_r\left(\frac{v_y-l_r{\omega }_r}{v_x}-{\delta }_r\right)-m{v}_x{\omega }_r\hfill \\ I_z{\dot{\omega }}_r=l_f{C}_f\left(\frac{v_y+l_f{\omega }_r}{v_x}-{\delta }_f\right)-l_r{C}_r\left(\frac{v_y-l_r{\omega }_r}{v_x}-{\delta }_r\right)+M_{zd}\hfill \end{array}\hfill \\ \dot{\phi }={\omega }_r\hfill \\ \dot{x}=v_x\cos \phi -v_y\sin \phi \hfill \\ \dot{y}=v_x\sin \phi +v_y\cos \phi \hfill \end{array}\right.$$

(17)

where, $\phi$ is the heading angle, $C_f$ and $C_r$ are the equivalent cornering stiffnesses of the front axle and the rear axle, respectively, ${\delta }_f$ is the front wheel angle, ${\delta }_r$ is the rear wheel angle, $M_{zd}$ is the generalized yaw moment, $x$ and $y$ are the longitudinal and lateral coordinates of the centroid in the geodetic coordinate system, respectively.

In the path tracking process, the desired heading angle and lateral coordinate can be obtained from the longitudinal coordinate of the current vehicle position. The state variable $X$, the control variable $u$ and the output variable $Y$ of the controller can be shown as:

$$\begin{array}{l}X={\left[v_y,\phi ,{\omega }_r,y,x\right]}^T\\ u={\left[{\delta }_f,{\delta }_r,M_{zd}\right]}^T\\ Y={\left[\phi ,y\right]}^T\end{array}$$

(18)

Assuming that the longitudinal velocity remains constant during path tracking, in order to meet the real-time requirements of high-speed control, this paper uses Taylor’s formula to linearize the nonlinear system. The Taylor expansion is performed at the reference point, and the high-order differential components are rounded off. Then, the discretization is carried out using the first-order difference quotient method, as shown in Equation (19).

$$X\left(k+1\right)=aX\left(k\right)+bu\left(k-1\right)+b\Delta u\left(k\right)+d\left(k\right)$$

(19)

where,

$$\begin{array}{l}a=I+T\frac{\partial f}{\partial X}, b=T\frac{\partial f}{\partial u}\\ d\left(k\right)=X_r\left(k+1\right)-a{X}_r\left(k\right)-b{u}_r\left(k\right)\end{array}$$

(20)

the augmented state variables are constructed:

$$\xi \left(k+1\right)=\left[\begin{array}{c}X\left(k+1\right)\\ u\left(k\right)\end{array}\right]=\left[\begin{array}{cc}a& b\\ 0_{N_u\times N_x}& I_{N_u\times N_u}\end{array}\right]\xi \left(k\right)+\left[\begin{array}{c}b\\ I_{N_u\times N_u}\end{array}\right]\Delta u\left(k\right)+\left[\begin{array}{c}d\left(k\right)\\ 0_{N_u\times 1}\end{array}\right]$$

(21)

where, $N_x$ is the number of state variables, and $N_u$ is the number of control variables. Equation (21) can be simplified as

$$\xi \left(k+1\right)=A\xi \left(k\right)+B\Delta u\left(k\right)+\overline{d}\left(k\right)$$

(22)

where,

$$A=\left[\begin{array}{cc}a& b\\ 0_{N_u\times N_x}& I_{N_u\times N_u}\end{array}\right], B=\left[\begin{array}{c}b\\ I_{N_u\times N_u}\end{array}\right], \overline{d}\left(k\right)=\left[\begin{array}{c}d\left(k\right)\\ 0_{N_u\times 1}\end{array}\right]$$

(23)

Therefore, the discretized output equation can be obtained as:

$$Y\left(k+1\right)=C\xi \left(k+1\right)$$

(24)

where,

$$C=\left[\begin{array}{cccccccc}0& 1& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 1& 0& 0& 0& 0\end{array}\right]$$

(25)

The prediction time-domain length is set as $p$, and the control time-domain length is set as $c$ ($p\ge c$). Assuming that the control increment outside the control time domain is zero, the predicted output equation of the system can eventually be obtained by iteration, as shown in Equation (26).

$$Y_p\left(k\right)=\Phi \xi \left(k\right)+\Theta \Delta U\left(k\right)+\Gamma D\left(k\right)$$

(26)

where,

$$Y_p\left(k\right)=\left[\begin{array}{l}Y\left(\left.k+1\right|k\right)\hfill \\ Y\left(\left.k+2\right|k\right)\hfill \\ Y\left(\left.k+3\right|k\right)\hfill \\ \begin{array}{c}⋮\\ Y\left(\left.k+c\right|k\right)\\ ⋮\end{array}\hfill \\ Y\left(\left.k+p\right|k\right)\hfill \end{array}\right], \Delta U\left(k\right)=\left[\begin{array}{c}\Delta u\left(k\right)\\ \Delta u\left(k+1\right)\\ \Delta u\left(k+2\right)\\ \begin{array}{c}⋮\\ \Delta u\left(k+c-1\right)\end{array}\end{array}\right], D\left(k\right)=\left[\begin{array}{c}\overline{d}\left(k\right)\\ \overline{d}\left(k+1\right)\\ \overline{d}\left(k+2\right)\\ ⋮\\ \overline{d}\left(k+c-1\right)\\ ⋮\\ \overline{d}\left(k+p-1\right)\end{array}\right]$$

(27)

$$\Phi =\left[\begin{array}{l}CA\hfill \\ C{A}^2\hfill \\ C{A}^3\hfill \\ \begin{array}{c}⋮\\ C{A}^c\\ ⋮\end{array}\hfill \\ C{A}^P\hfill \end{array}\right], \Theta =\left[\begin{array}{ccccc}CB& 0_{N_y\times N_u}& 0_{N_y\times N_u}& \cdots & 0_{N_y\times N_u}\\ CAB& CB& 0_{N_y\times N_u}& \cdots & 0_{N_y\times N_u}\\ C{A}^{2}B& CAB& CB& \cdots & 0_{N_y\times N_u}\\ ⋮& ⋮& ⋮& ⋮& ⋮\\ C{A}^{c-1}B& C{A}^{c-2}B& C{A}^{c-3}B& \cdots & CB\\ ⋮& ⋮& ⋮& ⋮& ⋮\\ C{A}^{p-1}B& C{A}^{p-2}B& C{A}^{p-3}B& \cdots & C{A}^{p-c}B\end{array}\right]$$

(28)

$$\Gamma =\left[\begin{array}{ccccccc}C& 0_{N_y\times \left(N_x+N_u\right)}& 0_{N_y\times \left(N_x+N_u\right)}& \cdots & 0_{N_y\times \left(N_x+N_u\right)}& \cdots & 0_{N_y\times \left(N_x+N_u\right)}\\ CA& C& 0_{N_y\times \left(N_x+N_u\right)}& \cdots & 0_{N_y\times \left(N_x+N_u\right)}& \cdots & 0_{N_y\times \left(N_x+N_u\right)}\\ C{A}^2& CA& C& \cdots & 0_{N_y\times \left(N_x+N_u\right)}& \cdots & 0_{N_y\times \left(N_x+N_u\right)}\\ ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮\\ C{A}^{c-1}& C{A}^{c-2}& C{A}^{c-3}& \cdots & C& \cdots & 0_{N_y\times \left(N_x+N_u\right)}\\ ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮\\ C{A}^{p-1}& C{A}^{p-2}& C{A}^{p-3}& \cdots & C{A}^{p-c}& \cdots & C\end{array}\right]$$

(29)

2.

Objective function design

The control objective of lateral path tracking is to ensure that the vehicle tracks the desired path smoothly, accurately, and quickly. Therefore, the objective function in this paper is designed as shown in Equation (30).

$$J={\left[Y_p\left(k\right)-Y_{ref}\right]}^T{Q}_Q\left[Y_p\left(k\right)-Y_{ref}\right]+\Delta U{\left(k\right)}^T{R}_R\Delta U\left(k\right)$$

(30)

where, $Q_Q=I_p\otimes Q$, $R_R=I_c\otimes R$, $Q$ is the weight coefficient matrix of the output of the control system, $R$ is the weight coefficient matrix of the control increment, and $\otimes$ represents the Kronecker product.

The first in the objective function represents the deviation of the heading angle and the deviation of the lateral displacement, which characterizes the accuracy of the vehicle in tracking the desired path. The second limits the control increment, which not only ensures the stability and continuity of the control but also takes into account the response capability of the actuator. Equation (30) can be simplified as,

$$J={\left[E+\Theta \Delta U\left(k\right)\right]}^{T}Q\left[E+\Theta \Delta U\left(k\right)\right]+\Delta U{\left(k\right)}^{T}R\Delta U\left(k\right)$$

(31)

where,

$$E=\Phi \xi \left(k\right)+\Gamma D\left(k\right)-Y_{ref}$$

(32)

Considering that E is a constant matrix at each sampling moment, it can be omitted. The final objective function in standard quadratic programming form is obtained, as shown in Equation (33).

$$J=\Delta U{\left(k\right)}^T\left({\Theta }^{T}Q\Theta +R\right)\Delta U\left(k\right)+2E^{T}Q\Theta \Delta U\left(k\right)$$

(33)

3.

Constraints design

A major advantage of MPC over other control methods is that it handles constraints better. MPC is a rolling optimization process that can naturally incorporate constraints into the optimization problem. For the path-following controller in this paper, the constraints are mainly considered in the form of the constraint on the control increment, the constraint on the control set, and the stability constraints.

In order to make the control process more stable and improve the stability and comfort of the vehicle when tracking the desired path, the control increment should be constrained. Due to the constraints of the actuator, there are also constraints on the control variables, such as front wheel angle, rear wheel angle, and generalized yaw moment. So, in summary, the following can be said:

$$U_{\min }\left(k\right)\le U_{ct}\left(k-1\right)+A_{ct}\ast \Delta U\left(k\right)\le U_{\max }\left(k\right)$$

(34)

where,

$$U_{ct}\left(k-1\right)=\left[\begin{array}{c}u\left(k-1\right)\\ u\left(k-1\right)\\ ⋮\\ u\left(k-1\right)\end{array}\right], A_{ct}=\left[\begin{array}{cccc}1& 0& \cdots & 0\\ 1& 1& \cdots & 0\\ ⋮& ⋮& ⋮& ⋮\\ 1& 1& 1& 1\end{array}\right]\otimes I_c$$

(35)

To ensure that the vehicle has good stability while tracking the desired path, stability constraints must be imposed on the MPC optimization problem. The stability constraints generally start with the two parameters yaw rate and tire slip angle. Therefore, the constraint output variable is defined as $Y_b={\left[{\omega }_r,{\alpha }_r\right]}^T$. The expression for the constraint output can be obtained, as shown in Equation (36).

$$Y_b\left(k+1\right)=C_b\xi \left(k+1\right)$$

(36)

where,

$$C_b=\left[\begin{array}{cccccccc}0& 0& 1& 0& 0& 0& 0& 0\\ \frac{1}{v_x}& 0& -\frac{l_r}{v_x}& 0& 0& 0& -1& 0\end{array}\right]$$

(37)

The predicted output equation for the stability constraint is shown in Equation (38).

$$Y_{p,b}\left(k\right)={\Phi }_b\xi \left(k\right)+{\Theta }_b\Delta U\left(k\right)+{\Gamma }_{b}D\left(k\right)$$

(38)

where,

$$Y_{p,b}\left(k\right)=\left[\begin{array}{l}{Y}_b\left(\left.k+1\right|k\right)\hfill \\ Y_b\left(\left.k+2\right|k\right)\hfill \\ Y_b\left(\left.k+3\right|k\right)\hfill \\ \begin{array}{c}⋮\\ Y_b\left(\left.k+c\right|k\right)\\ ⋮\end{array}\hfill \\ Y_b\left(\left.k+p\right|k\right)\hfill \end{array}\right], {\Phi }_b=\left[\begin{array}{l}{C}_{b}A\hfill \\ C_b{A}^2\hfill \\ C_b{A}^3\hfill \\ \begin{array}{c}⋮\\ C_b{A}^c\\ ⋮\end{array}\hfill \\ C_b{A}^P\hfill \end{array}\right]$$

(39)

$${\Theta }_b=\left[\begin{array}{ccccc}{C}_{b}B& 0_{N_{y1}\times N_u}& 0_{N_{y1}\times N_u}& \cdots & 0_{N_{y1}\times N_u}\\ C_{b}AB& C_{b}B& 0_{N_{y1}\times N_u}& \cdots & 0_{N_{y1}\times N_u}\\ C_b{A}^{2}B& C_{b}AB& C_{b}B& \cdots & 0_{N_{y1}\times N_u}\\ ⋮& ⋮& ⋮& ⋮& ⋮\\ C_b{A}^{c-1}B& C_b{A}^{c-2}B& C_b{A}^{c-3}B& \cdots & C_{b}B\\ ⋮& ⋮& ⋮& ⋮& ⋮\\ C_b{A}^{p-1}B& C_b{A}^{p-2}B& C_b{A}^{p-3}B& \cdots & C_b{A}^{p-c}B\end{array}\right]$$

(40)

$${\Gamma }_b=\left[\begin{array}{ccccccc}{C}_b& 0_{N_{y1}\times \left(N_x+N_u\right)}& 0_{N_{y1}\times \left(N_x+N_u\right)}& \cdots & 0_{N_{y1}\times \left(N_x+N_u\right)}& \cdots & 0_{N_{y1}\times \left(N_x+N_u\right)}\\ C_{b}A& C_b& 0_{N_{y1}\times \left(N_x+N_u\right)}& \cdots & 0_{N_{y1}\times \left(N_x+N_u\right)}& \cdots & 0_{N_{y1}\times \left(N_x+N_u\right)}\\ C_b{A}^2& C_{b}A& C_b& \cdots & 0_{N_{y1}\times \left(N_x+N_u\right)}& \cdots & 0_{N_{y1}\times \left(N_x+N_u\right)}\\ ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮\\ C_b{A}^{c-1}& C_b{A}^{c-2}& C_b{A}^{c-3}& \cdots & C_b& \cdots & 0_{N_{y1}\times \left(N_x+N_u\right)}\\ ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮\\ C_b{A}^{p-1}& C_b{A}^{p-2}& C_b{A}^{p-3}& \cdots & C_b{A}^{p-c}& \cdots & C_b\end{array}\right]$$

(41)

Based on the road adhesion conditions and the tire slip angle limits, it can be shown in Equation (42).

$$Y_{b,\min }\le {\Phi }_b\xi \left(k\right)+{\Theta }_b\Delta U\left(k\right)+{\Gamma }_{b}D\left(k\right)\le Y_{b,\max }$$

(42)

In summary, the finite time domain optimization problem for MPC in this paper can be transformed into a standard quadratic programming problem:

$$\begin{array}{l}\underset{\Delta U\left(k\right)}{\min }\Delta U{\left(k\right)}^T\left({\Theta }^{T}Q\Theta +R\right)\Delta U\left(k\right)+2E^{T}Q\Theta \Delta U\left(k\right)\\ \begin{array}{ll}s.t.\hfill & \Delta U_{\min }\left(k\right)\le \Delta U\left(k\right)\le \Delta U_{\max }\left(k\right)\hfill \\ \hfill & U_{\min }\left(k\right)\le U_{ct}\left(k-1\right)+A_{ct}\ast \Delta U\left(k\right)\le U_{\max }\left(k\right)\hfill \\ \hfill & Y_{b,\min }\le {\Phi }_b\xi \left(k\right)+{\Theta }_b\Delta U\left(k\right)+{\Gamma }_{b}D\left(k\right)\le Y_{b,\max }\hfill \end{array}\end{array}$$

(43)

The optimal sequence of control increments is obtained by solving the standard quadratic programming problem, and the first element $\Delta u^{\ast }\left(k\right)$ is output. The control variable at the current moment can be described as Equation (44).

$$u\left(k\right)=u\left(k-1\right)+\Delta u^{\ast }\left(k\right)$$

(44)

2.5. Coordinated Control

The coordinated control layer receives the generalized force, generalized yaw moment, and generalized wheel angle from the upper layer. The wheel torque and the wheel angle are calculated and output to in-wheel motors and steering motors. Depending on the actuating subsystem, the coordinated control layer consists of two parts: torque distribution control and wheel angle distribution control.

2.5.1. Torque Distribution Control

The torque distribution control receives the generalized longitudinal force and the generalized yaw moment output from the upper level and distributes and controls the driving torque of the individual wheels. The current research in torque distribution mainly considers the stability index and economic index as the control objective and uses online real-time optimization to complete the torque coordinated distribution [33,34]. However, current research has the following limitations:

• The rule-based torque distribution control does not consider the operating efficiency of motor, resulting in unnecessary power loss.
• The real-time optimization of torque distribution places a large burden on the controller. The solution speed may be slow and this problem may even be unsolvable under certain working conditions.
• The economic evaluation generally uses the size of the control variables as the index, ignoring the efficiency characteristics of the motor. This is mainly due to the nonlinearity of the motor model, which makes it impossible to optimize the system efficiency in real-time.

To address these issues, this paper simplifies the real-time optimization problem of torque distribution in order to balance practical and optimization objectives. A combination of rule-based strategy and offline optimization strategy is used for the distribution strategy. A rule-based strategy is designed to distribute the torque between the left and right wheels. Considering the stability constraints, the optimization of the torque distribution coefficients between the front and rear wheels based on the MPSO algorithm leads to an optimal wheel torque, improving the efficiency while ensuring the response speed of the vehicle. First, the torque distribution between left and right wheels is carried out based on the generalized longitudinal force and generalized yaw moment requirements. Then, the MPSO algorithm is used to optimize the torque distribution offline between the front and rear wheels on the same side to obtain the optimal distribution coefficient of the front axle. Finally, the driving torque of each wheel is calculated.

• Left-right distribution

A rule-based strategy is used for the generalized longitudinal force to distribute it equally between the left and right sides of the vehicle. The corresponding demand torques are obtained as follows.

$$T_{dlF}=T_{drF}=\frac{F_{xd}}{2}{R}_{eff}$$

(45)

where, $T_{dlF}$ and $T_{drF}$ are the left-hand and right-hand demand torques, respectively, corresponding to the generalized longitudinal forces.

The same uniform distribution is used for the generalized yaw moment. The yaw moment generated by the wheels on both sides is equal and has an opposite direction, which can shorten the response time of the vehicle and speed up the response speed. It can be represented as Equation (46).

$$\begin{array}{l}{T}_{dlM}=-\frac{M_{zd}}{d}{R}_{eff}\\ T_{drM}=\frac{M_{zd}}{d}{R}_{eff}\end{array}$$

(46)

where, $T_{dlM}$ and $T_{drM}$ are the left-hand and right-hand demand torques, respectively corresponding to the generalized yaw moment.

The total demand torques are,

$$\begin{array}{l}{T}_{dl}=\frac{F_{xd}}{2}{R}_{eff}-\frac{M_{zd}}{d}{R}_{eff}\\ T_{dr}=\frac{F_{xd}}{2}{R}_{eff}+\frac{M_{zd}}{d}{R}_{eff}\end{array}$$

(47)

where, $T_{dl}$ and $T_{dr}$ are the left-hand and right-hand demand torques, respectively.

2.

Front-rear distribution

For the front and rear wheels on the same side, the total demand torque is fixed. The independent drive function allows the torque of the front and rear wheels to be optimally distributed so that in-wheel motors can operate in the high-efficiency range as far as possible. Considering the high time complexity of online optimization, this paper uses the MPSO algorithm to optimize the distribution coefficient of the front axle offline, and it is put into a tabular form for subsequent table lookup operations which ensures the real-time performance of the system.

(1)

The MPSO algorithm

PSO (Particle Swarm Optimization) is a type of swarm intelligent optimization algorithm [35]. In PSO, each solution of the optimization problem is abstracted as a particle, and all particles search for the optimal solution in the solution space. Each particle is assigned a fitness function to determine the fitness of the current location, and a speed property to determine the distance and direction of flight, after which the optimal solution is determined by iteration.

The traditional PSO algorithm has the advantage of fast convergence, but it can easily fall into local optimal solutions in some complex situations. Therefore, the MPSO algorithm is introduced to avoid the problem of premature convergence. The MPSO algorithm combines the traditional PSO algorithm with the idea of mutation in the genetic algorithm. The mutation occurs when the population location is updated, thereby jumping out of the local optimal solution, which is conducive to finding the global optimal solution and reduces the probability of premature convergence. The process of the MPSO is shown in Figure 5.

(2)

Distribution coefficient optimization

Since the four in-wheel motors are the same, the offline optimization problems on the left and right sides are essentially the same, so this paper takes the right-side wheel as an example. In order to carry out the optimization, the economy needs to be characterized and this paper uses the overall motor efficiency as the economic index. The expression of the overall motor efficiency is,

$$J=\frac{T_{fr}{n}_{fr}+T_{rr}{n}_{rr}}{\frac{T_{fr}{n}_{fr}}{{\eta }_{fr}}+\frac{T_{rr}{n}_{rr}}{{\eta }_{rr}}}$$

(48)

where, $T_{fr}$ and $T_{rr}$ are the output torques of the right front and right rear motors, respectively, $n_{fr}$ and $n_{rr}$ are the output speeds of the right front and right rear motors, respectively, and ${\eta }_{fr}$ and ${\eta }_{rr}$ are the output efficiencies of the right front and right rear motors, respectively.

To simplify the complexity of the optimization problem, the concept of the front axle distribution coefficient $\lambda$ is introduced, which is the ratio of the front wheel torque to the total demand torque.

During normal driving, there is little difference in speed between the front and rear wheels on the same side, so the optimization problem can be translated into Equation (49).

$$\min J=\frac{\lambda }{{\eta }_{fr}}+\frac{1-\lambda }{{\eta }_{rr}}$$

(49)

Considering the constraints on the external characteristics of the in-wheel motor, the following constraints are made:

$$\begin{array}{l}0\le \lambda T_{dr}\le T_{\max }\\ 0\le \left(1-\lambda \right)T_{dr}\le T_{\max }\end{array}$$

(50)

where, $T_{\max }$ is the peak torque of the in-wheel motor.

In addition, depending on the vehicle model, the center of mass is closer to the front axle, which means that the vertical load on the front axle is greater than that on the rear axle. This results in a greater tire adhesion limit for the front wheel, with a greater longitudinal force limit. Therefore, the front wheels should take on a larger portion of the required torque. The following restriction applies to the distribution coefficient of the front axle.

$$0.5\le \lambda \le 1$$

(51)

Factors affecting the front axle distribution coefficient include the velocity and the demand torque, wherein the velocity is represented by wheel speed. Based on the MPSO algorithm, the offline optimization of the front axle distribution coefficient is realized by programming in MATLAB. Finally, the mapping of the optimal distribution coefficient for the front axle is shown in Figure 6. The optimal distribution coefficient for the front axle is determined by the velocity and the demand torque, which realizes the torque distribution between the front and rear wheels on the same side. From the optimization results, when the demand torque is low, the front wheel drive is selected; when the demand torque is high, it tends to be four-wheel drive.

(3)

Torque calculation for each wheel

Based on the optimal coefficient for the distribution of front axle torque determined via offline optimization, the torque for each wheel can be calculated. The offline optimization of the front axle distribution coefficients takes into account the constraints on the peak torque of the in-wheel motors but still requires stability corrections according to the constraints on road adhesion, as shown in Equation (52).

$$\left\{\begin{array}{l}{T}_{fl}=\min \left(\left|{\lambda }_l{T}_{dl}\right|,\mu F_{zfl}\right)\cdot \mathrm{sgn}\left(T_{dl}\right)\hfill \\ T_{fr}=\min \left(\left|{\lambda }_r{T}_{dr}\right|,\mu F_{zfr}\right)\cdot \mathrm{sgn}\left(T_{dr}\right)\hfill \\ T_{rl}=\min \left(\left|\left(1-{\lambda }_l\right)T_{dl}\right|,\mu F_{zrl}\right)\cdot \mathrm{sgn}\left(T_{dl}\right)\hfill \\ T_{rr}=\min \left(\left|\left(1-{\lambda }_r\right)T_{dr}\right|,\mu F_{zrr}\right)\cdot \mathrm{sgn}\left(T_{dr}\right)\hfill \end{array}\right.$$

(52)

where, ${\lambda }_l$ and ${\lambda }_r$ are the optimal front axle distribution coefficients for the left and right side, respectively. $T_{ij}$ is the torque of each wheel, $\mu$ is the coefficient of road adhesion, and $F_{zij}$ is the vertical load of each wheel, which can be obtained using Equation (53).

$$\left\{\begin{array}{c}{F}_{zfl}=\frac{mg{l}_r}{2l}-\frac{m{a}_{x}h}{2l}-\frac{m{a}_{y}h{l}_r}{dl}\\ F_{zfr}=\frac{mg{l}_r}{2l}-\frac{m{a}_{x}h}{2l}+\frac{m{a}_{y}h{l}_r}{dl}\\ F_{zrl}=\frac{mg{l}_f}{2l}+\frac{m{a}_{x}h}{2l}-\frac{m{a}_{y}h{l}_f}{dl}\\ F_{zrr}=\frac{mg{l}_f}{2l}+\frac{m{a}_{x}h}{2l}+\frac{m{a}_{y}h{l}_f}{dl}\end{array}\right.$$

(53)

2.5.2. Wheel Angle Distribution Control

The generalized wheel angle calculated by the target planning layer is the corresponding angle at the center of the front and rear axles, which is referred to as the equivalent angle of the front and rear wheels. Traditional control methods often assume that the left and right wheels on the same axis have the same angle, which directly corresponds to the equivalent angle. Although this method can simplify the design process, it does not take into account the actual steering geometry relationship, which not only increases tire wear but also causes unnecessary energy consumption. For this reason, the wheel angle distribution method based on Ackermann steering geometry is applied in this paper to convert the equivalent angle into the angle of each wheel.

By distributing the wheel angles based on the Ackermann steering geometry, the tires can be put into a pure rolling condition as much as possible. However, the ideal Ackermann steering principle is limited by the fact that it ignores the cornering characteristics of the tires. To further improve vehicle stability and reduce excessive tire wear, the Ackermann principle must be modified to obtain wheel angles.

When the vehicle is steering, the lateral elastic deformation of the tires results in the slip angles of the individual wheels on the front and rear axles. As shown in Equation (54), the corresponding relationship that each wheel angle should satisfy can be called the modified Ackermann principle.

$$\left\{\begin{array}{l}\tan ({\delta }_{fl}-{\alpha }_{fl})=\frac{\tan {\delta }_f}{1-\frac{d}{2l}(\tan {\delta }_f-\tan {\delta }_r)}\\ \tan ({\delta }_{fr}-{\alpha }_{fr})=\frac{\tan {\delta }_f}{1+\frac{d}{2l}(\tan {\delta }_f-\tan {\delta }_r)}\\ \tan ({\delta }_{rl}-{\alpha }_{rl})=\frac{\tan {\delta }_r}{1-\frac{d}{2l}(\tan {\delta }_f-\tan {\delta }_r)}\\ \tan ({\delta }_{rr}-{\alpha }_{rr})=\frac{\tan {\delta }_r}{1+\frac{d}{2l}(\tan {\delta }_f-\tan {\delta }_r)}\end{array}\right.$$

(54)

3. Simulation and Results

This section is divided into subheadings. It should provide a concise and precise description of the experimental results, their interpretation, as well as the experimental conclusions that can be drawn.

3.1. Environment and Configuration

In this section, the simulation model is built based on the MATLAB R2022a software produced by MathWorks (Natick, MA, USA), as shown in Figure 7. The vehicle model is the seven degrees-of-freedom vehicle model presented in Section 2. The simulations were performed using single-lane change and slalom test. The results verify the effect of trajectory tracking and economic optimization, which fully proves the effectiveness of the proposed control strategy.

To validate the economic performances of the coordinated control, simulations comparing the average distribution strategy and the wheel load distribution strategy are performed. These two strategies are simple and efficient and are widely used in practice. Therefore, the advantages of the proposed distribution control strategy can be further emphasized.

3.1.1. Average Distribution Strategy

The generalized longitudinal force and the generalized yaw moment output from the upper layer are distributed equally to each wheel. The wheel torque in this strategy is as follows:

$$\left\{\begin{array}{c}{T}_{fl}=T_{rl}=\frac{F_{xd}}{4}{R}_{eff}-\frac{M_{zd}}{2d}{R}_{eff}\\ T_{fr}=T_{rr}=\frac{F_{xd}}{4}{R}_{eff}+\frac{M_{zd}}{2d}{R}_{eff}\end{array}\right.$$

(55)

3.1.2. Wheel Load Distribution Strategy

The generalized force and yaw moment are distributed equally between the left and right sides of the vehicle. Then, the distribution between the front and rear wheels is carried out on the same side in proportion to the vertical load on the tires. The wheel torque with this strategy is shown in Equation (56).

$$\left\{\begin{array}{c}{T}_{fl}=\frac{F_{zfl}}{F_{zfl}+F_{zrl}}\left(\frac{F_{xd}}{2}{R}_{eff}-\frac{M_{zd}}{d}{R}_{eff}\right)\\ T_{fr}=\frac{F_{zfr}}{F_{zfr}+F_{zrr}}\left(\frac{F_{xd}}{2}{R}_{eff}+\frac{M_{zd}}{d}{R}_{eff}\right)\\ T_{rl}=\frac{F_{zrl}}{F_{zfl}+F_{zrl}}\left(\frac{F_{xd}}{2}{R}_{eff}-\frac{M_{zd}}{d}{R}_{eff}\right)\\ T_{rr}=\frac{F_{zrr}}{F_{zfr}+F_{zrr}}\left(\frac{F_{xd}}{2}{R}_{eff}+\frac{M_{zd}}{d}{R}_{eff}\right)\end{array}\right.$$

(56)

3.2. Results and Analysis

3.2.1. Single-Lane Change

As one of the common conditions, the single-lane change condition is usually used to simulate a vehicle lane change scenario.

• Trajectory tracking effect

Three speeds of 40 km/h, 80 km/h, and 120 km/h are simulated to verify the effect of the controller. The longitudinal velocity tracking and lateral path tracking results are analyzed to verify the performance of the vehicle at low, medium, and high speeds. The simulation results are shown in Figure 8 and Figure 9.

For a single-lane change, the longitudinal velocity tracking algorithm can track the expected velocity very well according to Figure 8. The tracking error is less than 0.2 km/h, which satisfies the longitudinal tracking accuracy requirements. As the velocity increases, the tracking error gradually increases. This is mainly because the nonlinearity of the tire gradually increases and is associated with the strong coupling properties between the subsystems, resulting in a decrease in tracking accuracy.

The data in Figure 9c,d have been processed to obtain the tracking error at different velocities to facilitate quantitative analysis of the error, as shown in

Table 1. Comparison of path tracking errors under single-lane change.

	Performance Index	40 km/h	80 km/h	120 km/h
Lateral displacement tracking error (m)	Maximum	0.0115	0.0171	0.0234
	Average	0.0024	0.0036	0.0053
	Standard deviations	0.0040	0.0058	0.0075
Heading angle tracking error (rad)	Maximum	0.0012	0.0036	0.0042
	Average	0.0002	0.0009	0.0009
	Standard deviations	0.0004	0.0012	0.0012

.

From Figure 9 and Table 1, it can be seen that the lateral path tracking algorithm is very good at tracking the desired path for single-lane change. As the speed increases, the lateral displacement error gradually increases. The maximum lateral displacement error is 0.0115 m at 40 km/h, 0.0171 m at 80 km/h, and 0.0234 m at 120 km/h, which is mainly due to the gradual increase in the nonlinearity of the tire and the enhancement of the coupling effect of the subsystem. At the same time, although the lateral displacement error increases, it is still generally small and within acceptable limits. Similarly, the heading angle error also increases with velocity but remains at a low level overall. The maximum heading angle error at 120 km/h is only 0.0042 rad, which demonstrates the good path-tracking performance of the controller.

The trajectory tracking controller is designed to achieve good stability while providing accurate and reliable path tracking. To check the stability effect, the lateral acceleration of the vehicle, side slip angle, and yaw rate are selected as indicators, and the simulation results are analyzed.

From Figure 10, Figure 11, Figure 12 and Figure 13, it can be seen that the vehicle is in a stable state due to the stability constraints and no instability occurs even at a high speed of 120 km/h. The lateral acceleration is kept within 0.4 g, which ensures the accuracy of the linear tire model. The side slip angles at different speeds are all within 0.01 rad, and the yaw rate is all within 0.15 rad/s, indicating that the vehicle has good driving stability.

The simulation results shown in Figure 13 illustrate the advantages of 4WIS in terms of dynamic control. At low speeds, the front and rear wheels rotate in the opposite direction, reducing the steering radius and improving the vehicle’s maneuverability. At high speeds, the front and rear wheels turn in the same direction, thereby increasing the lateral stability margin and improving the vehicle’s driving stability.

2.

Economy optimization effect

The above three distribution strategies are simulated at different velocities to validate the economy optimization effect, using the overall efficiency of the motor and the battery SOC as performance indicators. It is assumed that Rule 1 is the average distribution strategy, Rule 2 is the wheel load distribution strategy, and Rule 3 is the torque distribution strategy proposed in this paper. At the same time, the control strategy does not include a braking energy recovery strategy. The initial value of SOC is set to 0.8. The simulation results are shown in Figure 14, Figure 15 and Figure 16 and

Table 2. Comparison of overall motor efficiency under single-lane change.

Velocity (km/h)	Strategy	Maximum	Average
40	Rule 1	0.7356	0.7303
	Rule 2	0.7347	0.7283
	Rule 3	0.8411	0.8375
80	Rule 1	0.8969	0.8877
	Rule 2	0.8970	0.8882
	Rule 3	0.9276	0.9241
120	Rule 1	0.8973	0.8791
	Rule 2	0.8972	0.8793
	Rule 3	0.9147	0.9107

.

For a single-lane change, the change in SOC is almost identical for the average and wheel load distribution strategies at different speeds. The maximum and average values of the overall motor efficiency are also almost identical, indicating that the economic effects of these two strategies are almost the same.

The change in battery SOC represents the electrical energy consumed by the in-wheel motors while driving. Comparing the change curves of SOC at different speeds, the SOC change in Rule 3 is smaller than those in Rule 1 and Rule 2, and the SOC in Rule 3 falls relatively slowly. This shows that the in-wheel motors consume less electric energy in the Rule 3 strategy, which proves the energy economy of the proposed strategy.

The average value of the overall motor efficiency of Rule 3 is 14.68% higher than that of Rule 1 at 40 km/h. Compared with the results of Rule 1, the average values at 80 km/h and 120 km/h have increased by 4.10% and 3.5%, respectively. This shows that the coordinated control strategy proposed in this paper can maintain the overall motor efficiency at a relatively high level and effectively improve the vehicle economy.

Comparing the effects of economy optimization at different speeds, it can be found that the economy optimization at low speeds is better, and the difference in SOC decline is significant. This is because the demand torque is relatively low at low speeds. In the average distribution strategy, the torque of each wheel is very low, so in-wheel motors operate in the low-efficiency region. However, the proposed strategy optimizes the torque distribution and tends to select the front wheel drive mode. This allows the front in-wheel motors to operate in a relatively efficient range, and the overall efficiency can be significantly improved. However, at high speeds, the in-wheel motors’ speeds are also high. As shown in Figure 2, the motor’s efficiency characteristic range gradually narrows as the speed increases when the motor’s speed is high. The distribution difference of the motor working points decreases for different strategies, resulting in a less significant improvement effect.

3.2.2. Slalom Test

• Trajectory tracking effect

Two speeds of 30 km/h and 60 km/h are simulated to verify the effect of the controller. The longitudinal velocity tracking and lateral path tracking results are analyzed to verify the performance of the vehicle at different speeds. The simulation results are shown in Figure 17 and Figure 18.

The data in Figure 18c,d have been processed to obtain the tracking error at different velocities to facilitate quantitative analysis of the error, as shown in

Table 3. Comparison of path tracking errors under slalom test.

	Performance Index	30 km/h	60 km/h
Lateral displacement tracking error (m)	Maximum	0.0412	0.0603
	Average	0.0158	0.0241
	Standard deviations	0.0147	0.0214
Heading angle tracking error (rad)	Maximum	0.0058	0.0129
	Average	0.0011	0.0033
	Standard deviations	0.0010	0.0029

.

From Figure 18 and Table 3, it can be seen that the lateral path tracking algorithm is very good at tracking the desired path in a slalom test. As the speed increases, the path-tracking error gradually increases. The average lateral displacement error is 0.0158 m and the average heading angle error is 0.0011 rad at 30 km/h, while the average lateral displacement error is 0.0241 m and the average heading angle error is 0.0033 rad at 60 km/h. Although the tracking error increases, it remains low overall, which shows the good tracking performance of the designed controller.

Figure 19, Figure 20, Figure 21 and Figure 22 show that the vehicle is in a stable state due to stability constraints. The side slip angle is maintained at a low level, which provides good tracking capability. The lateral acceleration is kept within 0.4 g, ensuring the accuracy of the linear tire model. The velocities under both simulation conditions are not particularly high, so the front and rear wheels turn in opposite directions, which improves the maneuverability of the vehicle while meeting the stability requirements.

2.

Economy optimization effect

The above three distribution strategies are simulated at different speeds to validate the economy optimization effect, using the overall efficiency of the motor and the battery SOC as performance indicators. The simulation results are shown in Figure 23 and Figure 24 and

Table 4. Comparison of overall motor efficiency under slalom test.

Velocity (km/h)	Strategy	Maximum	Average
30	Rule 1	0.6404	0.6151
	Rule 2	0.6409	0.6121
	Rule 3	0.7857	0.7484
60	Rule 1	0.8973	0.7995
	Rule 2	0.8971	0.7987
	Rule 3	0.9345	0.8829

.

In the slalom test, the change in SOC is almost identical for the average and wheel load distribution strategies at different speeds. The maximum and average values of the overall motor efficiency are also almost identical, indicating that the economic effects of these two strategies are almost the same.

Comparing the change curves of SOC at different speeds, the change in SOC in Rule 3 is smaller than in Rule 1 and Rule 2, and the SOC in Rule 3 drops relatively slowly. This shows that the in-wheel motors consume electric energy under the Rule 3 strategy, which proves the energy economy of the proposed strategy.

The coordinated control proposed in this paper is significantly better than the other two strategies, as shown in Table 4. Compared with the results of Rule 1, the average value of the overall motor efficiency of Rule 3 has increased by 21.67% at 30 km/h. At 60 km/h, it has increased by 10.43%, which proves that the coordinated control strategy can effectively improve the vehicle economy.

4. Conclusions

In this paper, a 4WID-4WIS EV trajectory tracking coordinated control strategy considering energy consumption economy is proposed to improve vehicle stability and economy during trajectory tracking. A hierarchical chassis coordinated control architecture for 4WID-4WIS EV is designed. In the target planning layer, the longitudinal velocity tracking and lateral path tracking are achieved by using expert PID and MPC, respectively, considering stability constraints, such as yaw rate and tire slip angle. In the coordinated control layer, the optimal torque distribution for each wheel is performed to achieve the optimal overall motor efficiency based on the MPSO algorithm. Then, the angle distribution of each wheel is performed in combination with the modified Ackermann theory. From the simulation results, it can be concluded that the proposed coordinated control strategy not only achieves good trajectory tracking but also ensures driving stability and improves the energy consumption economy.

In future work, the coupling characteristics between multiple subsystems will be considered. The influence of vertical dynamics will be incorporated into the control strategy to achieve unified control of multi-dimensional lateral, longitudinal, and vertical dynamics. At the same time, the motor model will be further improved, and research on redundant control of actuator failures will be conducted. The test with a real vehicle will be conducted to further verify the effectiveness of the control strategy.