**1. Introduction**

With the continuous increase in car ownership, the modern intelligent transportation system has increasingly regarded vehicle safety as its key factor. Between 2002 and 2012, the lack of proper obstacle avoidance contributed to the deaths of millions of people in traffic accidents around the world, and the economic cost of these accidents amounted to \$0.48 trillion [1]. More than 90% of these accidents are caused by human factors [2]. Therefore, in order to avoid vehicle collisions and minimize the impact of accidents, a method of adding an increasing proportion of active safety systems is proposed [3].

The obstacle avoidance trajectory tracking control of intelligent vehicles is one of the key functions. According to the trajectory planned by the upper controller and the realtime state information, real-time vehicle control variables, such as front wheel angle and driving force/braking force, are generated [4]. Compared with traditional vehicles, the drive motors are directly installed in the drive wheels for the distributed drive electric vehicles. So they have outstanding advantages, such as short drive transmission chain, high transmission efficiency and compact structure. It can be more beneficial to realize the tracking control of the obstacle avoidance path by optimizing the distribution of the driving torque of each wheel [5]. Although distributed hub motors can lead to greater unsprung mass, a small amount of unsprung mass increase is negligible, given their advantages in dynamic control.

After the continuous development of path tracking control, different tracking control algorithms have been proposed and verified, such as sliding mode control, fuzzy control, and model predictive control. Reference [6] proposed a distributed drive unmanned vehicle path tracking and stability coordinated control strategy based on the layered control theory. In order to reduce the heading deviation and lateral deviation during the path-tracking

**Citation:** Wu, H.; Zhang, H.; Feng, Y. MPC-Based Obstacle Avoidance Path Tracking Control for Distributed Drive Electric Vehicles. *World Electr. Veh. J.* **2022**, *13*, 221. https://doi.org/ 10.3390/wevj13120221

Academic Editor: Yong Li

Received: 6 October 2022 Accepted: 16 November 2022 Published: 22 November 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

process, the sliding mode control method was used. Moreover, it can ensure the stability of the vehicle. On the basis of the traditional dynamic model and preview model, Zhang Bingli [7] designed a new trajectory-tracking controller based on neural networks and fuzzy control theory. Wang Yao [8] proposed a model prediction-fuzzy (MPC-Fuzzy) joint control strategy, which alleviated the conflict between the vehicle tracking accuracy and the controller's computational pressure under a single model prediction control and better solved the distributed driving vehicle path, tracking accuracy and stability issues in tracking. Hu [9] designed a fuzzy linear quadratic regulator (LQR) with preview PID angle compensation and verified that it can maintain good accuracy and stability at different vehicle speeds. However, it also had the problem of low tracking accuracy when the vehicle was running at a low speed. Kou [10] used a state expansion MPC and corner compensation fuzzy controller, designed a dual feedback MPC controller based on state expansion and verified that it had better path tracking performance at medium and low vehicle speeds. Wu [11] proposed a front-wheel steering controller and optimized the controller parameters so that the vehicle had better path-tracking performance at low and medium speeds and ensured that the vehicle had certain steering accuracy and driving stability when driving at high speed on ice and snow roads. In order to improve the stability of distributed drive electric vehicles in extreme working conditions, Li [12] proposed a nonlinear model predictive control (NMPC) based torque coordination control strategy. In order to coordinate the trajectory tracking accuracy and stability of distributed drive electric vehicles and improve the adaptive ability of the control algorithm to different working conditions, Zhuang [13] proposed a trajectory tracking control strategy based on Takagi–Sugeno fuzzy model predictive control (T-S FMPC).

The total required torque is calculated through the yaw moment obtained by the upper controller. The torque distributed to each drive wheel is obtained through different algorithms. There are many distribution methods for torque distribution. A fixed proportional distribution method was proposed in the literature [14,15]. Although this method was simple, it cannot implement different proportional distributions according to different road conditions. There was no optimization in the energy utilization efficiency of the motor, and the best distribution result could not be achieved. Ono E [16] proposed to realize the torque distribution to the four wheels according to the usage of the tire force of four tires. This method made the wear degree of the four tires more balanced, and the life of the tires and the motor were improved, but the performance of vehicle handling was greatly reduced. Mokhiamar O [17] also proposed to analyze the torque distribution problem by establishing a multi-objective optimization function, but this method may show different effects in different models and different operating conditions, so the consistency and robustness were poor. Reference [18] considered the dynamic characteristics of the vehicle when it is unstable and proposed a control allocation scheme based on feedback linearization. This scheme was relatively dependent on the accuracy of the model and had poor stability. Reference [19] distributed torque according to the vertical load distribution ratio of each wheel. Since the friction ellipse of the tire increases with the increase of the vertical load, the torque distribution according to the vertical load ratio of a single tire can avoid the saturation of the tire force and can simply and effectively improve the driving stability of the vehicle.

Based on the above analysis, the whole vehicle control process designed in this paper is shown in Figure 1. In Section 1, using the obstacle information obtained by the environment awareness layer and the vehicle pose information detected by the vehicle sensor, the obstacle avoidance path is planned. Anti-collision and anti-rollover constraints adding to a six-degree polynomial are used for path planning. The front wheel angle δf and additional yaw moment Δ Mz are output by the MPC controller in Section 3.2. By processing the information of the path planning layer, the torque distribution controller distributes the torque through the additional yaw moment and the wheel vertical force ratio in Section 3.3. The MPC controller and the torque distribution controller together form a path-tracking controller. During the obstacle avoidance process, the vehicle avoids obstacles through

the front wheel angle and four-wheel torque output by the controller. In Section 4, the simulation platform is established, and the effect of the controllers is verified.

**Figure 1.** Vehicle control flow chart. Where, *Xc*0 and *Yc*0 are the lateral distance and longitudinal distance from the obstacle; *vc* is the vehicle speed of the obstacle; *Xd* and *Yd* are the expected lateral and longitudinal distance of the vehicle; *ϕ*d is the expected yaw angle of the vehicle; Δ*Mz* is the additional yaw moment; *Tfl, Trl*, *Trr* and *Tfr* are the torques of the left front wheel, the left rear wheel, the right rear wheel and the right front wheel, respectively; *X* and *Y* are the transverse and vertical coordinates under the inertial coordinate system; *ϕ* is the yaw angle under the body coordinate system; *r* is the turning radius of the vehicle; *vx* and *vy* are the lateral and speed longitudinal of the vehicle; *ax* is the lateral acceleration of the vehicle.

#### **2. Path Planning**

In this paper, the six-degree polynomial path planning method [20] for a given time is used. In the sixth-degree polynomial path planning method, the vehicle's gentle arrival from the starting point to the target point is considered, the collision between the vehicle and the obstacle in the process of running is considered, and the smooth running according to the planned path is also considered. With the addition of anti-collision and anti-rollover conditions, the safety of obstacle avoidance path planning is improved.

Based on the longitudinal and lateral displacements at a given time, the initial expression for the sixth-degree polynomial is:

$$\begin{cases} X(t) = c\_0 + c\_1t + c\_2t^2 + c\_3t^3 + c\_4t^4 + c\_5t^5 + c\_6t^6\\\ Y(t) = b\_0 + b\_1t + b\_2t^2 + b\_3t^3 + b\_4t^4 + b\_5t^5 + b\_6t^6 \end{cases} \tag{1}$$

Adding anti-collision conditions [20]: set the coordinate center of the obstacle as (*Xob*, *Yob*), and the radius of the circumscribed circle of the vehicle and the obstacle as *Rcar* and *Rob*, this is shown in Figure 2,the constraint condition of anti-collision can be expressed as:

$$\left(R\_{\rm car} + R\_{\rm ob}\right)^2 < \left[X(t) - X\_{\rm obj}\right]^2 + \left[Y(t) - Y\_{\rm obj}\right]^2\tag{2}$$

**Figure 2.** Schematic diagram of collision avoidance.

Adding the anti-rollover condition constraint [20].

According to the simplified model diagram of vehicle rollover in Figure 3, the moment balance equation can be expressed as:

$$m\_s a\_y h + m\_s g \Delta y = (F\_{zI} - F\_{zr}) \frac{s}{2} \tag{3}$$

where, *ay* is the lateral acceleration, *ms* is the sprung mass, *h* is the distance from the center of mass to the ground, *Fzl* represents the sum of forces on the left front and rear wheels in the vertical direction, *Fzr* represents the sum of forces on the right front and rear wheels, *s* represents the track width, *g* is the gravitational acceleration, Δ*y* = *h* · sin *φ*.

**Figure 3.** Simplified model of automobile rollover.

The lateral transfer rate of the vehicle is:

$$\text{LTR} = \frac{F\_{\text{xl}} - F\_{\text{xx}}}{F\_{\text{xl}} + F\_{\text{xx}}} \tag{4}$$

Combining Equations (3) and (4), LTR can be expressed as:

$$LTR = \frac{2R}{s} \left[ \frac{a\_y}{\mathcal{g}} + \sin \phi \right] \tag{5}$$

where, *φ* is the vehicle roll angle.

The larger the absolute value of LTR indicates, the greater the rollover risk of the car, so the size of *ay* should be fully taken into account in path planning. Given the maximum roll risk LTRmax,

$$a\_{\mathcal{Y}} = \ddot{\mathcal{Y}}(t) \le \left(\frac{\mathbf{s} \cdot \text{LTR}\_{\text{max}}}{2\text{h}} - \sin\phi\right)\mathbf{g} \tag{6}$$

$$a\_{y\text{max}} = \left(\frac{\text{s} \cdot \text{LTR}\_{\text{max}}}{2h} - \sin\phi\right) \text{g} \tag{7}$$

Constructing the evaluation index function of anti-rollover:

$$J\_{anti-r} = \int\_{t\_0}^{t\_d} - \left(a\_{y\text{max}} - \left|a\_y\right|\right) dt \tag{8}$$

The polynomial optimization path planning problem is expressed as:

$$\begin{cases} \min l\_{anti-r} \\ \text{s.t.} (R\_{car} + R\_{ob})^2 < [X(t) - X\_{ob}]^2 + [Y(t) - Y\_{ob}]^2 \\ a\_y = \ddot{Y}(t) \le \left(\frac{s \cdot \text{LR}\_{\text{max}}}{2h} - \sin \phi\right) \text{g} \\ b\_3 t\_d^{-3} + b\_4 t\_d^{-4} + b\_5 t\_d^{-5} + b\_6 t\_d^{-6} = H \\ 3b\_3 t\_d^2 + 4b\_4 t\_d^3 + 5b\_5 t\_d^4 + 6b\_6 t\_d^5 = 0 \\ 6b\_3 t\_d + 12b\_4 t\_d^2 + 20b\_5 t\_d^3 + 30b\_6 t\_d^4 = 0 \end{cases} \tag{9}$$

where, *t*0 and *td* are the initial time and the end time, respectively, *ay* is the lateral acceleration, *Xob* and *Yob* are the coordinate centers of the obstacles, the coordinate center of

the obstacle:(*Xob*, *Yob*) = (100, 1), *Rcar* and *Rob* are the radiuses of the circumscribed circles between the vehicle and the obstacle, ∅ is the body roll angle, and *w* is the wheelbase, *h* is the distance from the center of mass to the ground, *g* is the acceleration of gravity, *H* is the self-conjugate matrix, *Janti-r* is the anti-rollover evaluation index function, and *LTR* is the roll risk coefficient.

In Equation (9), the coefficients *b*3, *b*4, *b*5 and *b*6 can be obtained by the optimization algorithm. The yaw angle can be expressed as:

$$\varphi(t) = \arctan \frac{\mathrm{d}Y(t)}{\mathrm{d}X(t)} = \arctan \frac{\dot{Y}(t)}{\dot{X}(t)}\tag{10}$$

#### **3. Path Tracking Control**

#### *3.1. Establishment of Vehicle Dynamics Model*

The three-degree-of-freedom vehicle dynamics model is shown in Figure 4, including lateral motion, yaw motion and longitudinal motion [21,22]. The aerodynamics, suspension system and steering system are ignored. Considering that the sideslip angle of the center of mass is small when the lateral acceleration is not high, the control accuracy will not be affected if it is ignored, so the sideslip angle of the center of mass is ignored [23]. According to the characteristics of the distributed drive vehicle, the additional yaw moment is considered.

**Figure 4.** Three-degree-of-freedom vehicle dynamics model.

In Figure 4, the coordinate axis *oxyz* is the fixed coordinate system of the vehicle, *xoz* is in the left-right symmetrical plane of the vehicle, the coordinate origin of the center of mass of the vehicle is *o*, the x-axis is the coordinate axis along the longitudinal direction of the vehicle, and the y-axis is the coordinate perpendicular to the longitudinal direction of the vehicle axis. The *z*-axis is the coordinate axis perpendicular to the *xoy* plane, and the direction satisfies the right-hand rule. *XOY* is the geodetic coordinate system. Considering that *Fxf* and *Fyf* are the force exerted on a single tire, it needs to be multiplied by 2. According to Newton's second law:

$$\begin{cases} m\ddot{\mathbf{x}} = m\dot{\mathbf{y}}\dot{\boldsymbol{\varphi}} + 2F\_{xf} + 2F\_{xr} \\ m\ddot{\mathbf{y}} = -m\dot{\mathbf{x}}\dot{\boldsymbol{\varphi}} + 2F\_{yf} + 2F\_{yr} \\ I\_z \ddot{\boldsymbol{\varphi}} = 2I\_f F\_{yf} - 2I\_r F\_{yr} \end{cases} \tag{11}$$

where *lf* and *lr* are the distances from the center of mass to the front and rear axles, respectively. *m* is the overall mass of the vehicle, and *Iz* is the central moment of inertia of the vehicle around the *z* axis.

The conversion relationship between the resultant force of the tire in the x-direction and y-direction, the longitudinal force and the lateral force is expressed as follows:

$$\begin{cases} F\_{xf} = F\_{lf}\cos\delta\_f - F\_{cf}\sin\delta\_f\\ F\_{xr} = F\_{lr}\cos\delta\_f + F\_{cr}\sin\delta\_r\\ F\_{yf} = F\_{lf}\sin\delta\_f + F\_{cf}\cos\delta\_f\\ F\_{yr} = F\_{cr}\cos\delta\_r - F\_{lr}\sin\delta\_f \end{cases} \tag{12}$$

where *Fxf* and *Fyf* are the resultant forces on the front wheels in the x and y directions, respectively. *Fxr* and *Fyr* are the resultant forces on the front wheels in the x and y directions, respectively, and *Flf* and *Flr* are the longitudinal forces on the front and rear wheels, respectively. *Fcf* and *Fcr* are the lateral forces on the front and rear wheels, respectively. *δf* is the front wheel angle, *δr* is the rear wheel angle.

The lateral and longitudinal forces of vehicle tires can be expressed as complex functions of parameters such as wheel slip angle, slip rate, road friction coefficient and vertical load:

$$F\_{\mathbb{C}} = f\_{\mathbb{C}}(\mathfrak{a}, \mathfrak{s}, \mathfrak{a}, F\_{\mathbb{Z}}) \tag{13}$$

$$F\_l = f\_l(\mathfrak{a}, \mathfrak{s}, \mathfrak{\mu}, F\_{\mathfrak{z}}) \tag{14}$$

where, *Fz* is the vertical load on the tire, *μ* is the road friction coefficient, *s* is the slip ratio, and *α* is the wheel side angle.

In order to simplify the model, the tire force is usually represented by a linear function under the condition of small tire longitudinal slip rate, lateral acceleration and slip angle. The force is expressed as follows:

⎧

⎨

⎩

$$\begin{cases} \quad F\_l = \mathbb{C}\_l s \\ \quad F\_\mathfrak{c} = \mathbb{C}\_\mathfrak{c} \mathfrak{a} \end{cases} \tag{15}$$

where *Cl* is the longitudinal stiffness of the vehicle tire, and *Cc* is the cornering stiffness of the vehicle tire.

The following approximate relationship is satisfied:

$$\begin{cases} \cos \theta \approx 1\\ \sin \theta \approx \theta\\ \tan \theta \approx \theta \end{cases} \tag{16}$$

where *θ* can be expressed as front wheel side angle and rear wheel side slip angle, etc.

The wheel side slip angle *α* can be obtained according to the geometric relationship:

$$\alpha = \tan^{-1} \frac{v\_c}{v\_l} \tag{17}$$

where *vc* and *vl* represent the lateral and longitudinal speeds of the tire, respectively, which can be represented by the speeds *vx* and *vy* in the direction of the coordinate system:

$$\begin{cases} \upsilon\_l = \upsilon\_y \sin \delta + \upsilon\_x \cos \delta\\ \upsilon\_c = \upsilon\_y \cos \delta - \upsilon\_x \sin \delta \end{cases} \tag{18}$$

where *δ* is the tire deflection angle. The tire speed is generally calculated by the speed conversion of the vehicle, and the conversion relationship is shown in Formula (19):

$$\begin{cases} \ v\_{yf} = \dot{y} + l\_f \dot{\phi} \\ \ v\_{yr} = \dot{y} - l\_r \dot{\phi} \end{cases} \tag{19}$$

$$\begin{cases} \ v\_{xf} = \dot{\mathbf{x}} \\ \ v\_{x\mathbf{r}} = \dot{\mathbf{x}} \end{cases} \tag{20}$$

Through the above formulas, the wheel side slip angle can be obtained as:

$$\begin{cases} \; \alpha\_f = \frac{\dot{y} + l\_f \dot{q}}{\dot{x}} - \delta\_f\\ \; \alpha\_r = \frac{\dot{y} - \dot{l}\_r \dot{q}}{\dot{x}} \end{cases} \tag{21}$$

From Equation (20), the tire lateral force of the front and rear wheels can be obtained as:

$$\begin{cases} \ \ \ \ \ \ \ \mathbf{F}\_{cf} = \mathbf{C}\_{cf} \left( \delta\_f - \frac{\dot{y} + l\_f \dot{\varphi}}{\dot{x}} \right) \\ \ \ \ \ \ \ \ \ \mathbf{F}\_{cr} = \mathbf{C}\_{cr} \frac{l\_r \dot{\varphi} - \dot{y}}{\dot{x}} \end{cases} \tag{22}$$

According to the above derivation formula, the following nonlinear vehicle dynamics model can be obtained as follows:

$$\begin{cases} m\ddot{y} = -m\dot{x}\dot{\varphi} + 2\left[\mathbb{C}\_{cf}\left(\delta\_{f} - \frac{\dot{y} + l\_{f}\dot{\varphi}}{\dot{x}}\right) + \mathbb{C}\_{cr}\frac{l\_{f}\dot{\psi} - \dot{y}}{\dot{x}}\right] \\\ m\ddot{x} = m\dot{y}\dot{\varphi} + 2\left[\mathbb{C}\_{lf}s\_{f} + \mathbb{C}\_{cf}\left(\delta\_{f} - \frac{\dot{y} + l\_{f}\dot{\varphi}}{\dot{x}}\right)\delta\_{f} + \mathbb{C}\_{lf}s\_{r}\right] \\\ l\_{z}\ddot{\varphi} = 2\left[l\_{f}\mathbb{C}\_{cf}\left(\delta\_{f} - \frac{\dot{y} + l\_{f}\dot{\varphi}}{\dot{x}}\right) - l\_{f}\mathbb{C}\_{cr}\frac{l\dot{\varphi} - \dot{y}}{\dot{x}}\right] + \Delta M\_{z} \\\ \dot{Y} = \dot{x}\sin\varphi + \dot{y}\cos\varphi \\\ \dot{X} = \dot{x}\cos\varphi - \dot{y}\sin\varphi \end{cases} \tag{23}$$

In this paper, due to the use of a linear time-varying model predictive control algorithm, the displacement and velocity in the longitudinal direction are not considered as control variables. Let .*ϕ* = r, .*x* = vx, .*y* = vy, Equation (22) can be simplified as:

$$\begin{cases} \dot{\upsilon}\_{y} = -\frac{\mathsf{C}\_{cf} + \mathsf{C}\_{cr}}{m\upsilon\_{x}} \upsilon\_{y} - \left(\frac{l\_{f}\mathsf{C}\_{cf} - l\_{r}\mathsf{C}\_{cr}}{m\upsilon\_{x}} + \upsilon\_{x}\right)r + \frac{\mathsf{C}\_{cf}}{m}\delta\_{\mathsf{f}}\\ \dot{r} = -\frac{l\_{f}\mathsf{C}\_{cf} - l\_{r}\mathsf{C}\_{cr}}{l\_{z}\upsilon\_{x}} \upsilon\_{y} - \frac{l\_{f}^{2}\mathsf{C}\_{cf} + l\_{r}^{2}\mathsf{C}\_{cr}}{l\_{z}\upsilon\_{x}}r + \frac{l\_{f}\mathsf{C}\_{cf}}{l\_{z}}\delta\_{f} + \frac{\Delta M\_{z}}{l\_{z}}\\ \dot{\varrho} = r\\ \dot{Y} = \upsilon\_{x}\varrho + \upsilon\_{y} \end{cases} \tag{24}$$

According to the above dynamic model, the following state space equation is established:

$$\begin{cases} \dot{\xi} = A\xi + Bu\\ y = C\xi \end{cases} \tag{25}$$

where the state quantity *ξ* = *vy r ϕ y*T, the control quantity *u* = *δf* Δ*Mz*T, **A**, **B**, and **C** are coefficient matrices, respectively:

$$A = \begin{bmatrix} -\frac{\mathbf{C}\_{cf} + \mathbf{C}\_{cr}}{m\mathbf{v}\_x} & -\left(\frac{l\_f \mathbf{C}\_{cf} - l\_r \mathbf{C}\_{cr}}{m\mathbf{v}\_x} + \mathbf{v}\_x\right) & \mathbf{0} & \mathbf{0} \\ -\frac{l\_f \mathbf{C}\_{cf} - l\_r \mathbf{C}\_{cr}}{l\_z \mathbf{v}\_x} & -\frac{l\_f^2 \mathbf{C}\_{cf} + l\_r^2 \mathbf{C}\_{cr}}{l\_z \mathbf{v}\_x} & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \mathbf{1} & \mathbf{0} & \mathbf{0} \\ \mathbf{1} & \mathbf{0} & \mathbf{v}\_x & \mathbf{0} \end{bmatrix} \tag{26}$$

$$\mathbf{B} = \begin{bmatrix} \frac{\mathbb{C}\_{cf}}{m} & 0\\ \frac{l\_f \mathbb{C}\_{cf}}{l\_z} & \frac{1}{l\_z} \\ 0 & 0 \\ 0 & 0 \end{bmatrix} \tag{27}$$

$$\mathcal{C} = \begin{bmatrix} 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{bmatrix} \tag{28}$$

#### *3.2. MPC Controller Design*

The linear time-varying model is shown in Figure 5. In the model, the control value is initial velocity *v0* = 18 m/s; The actual measured value is *vy*, *r*, *ϕ*, *y*; The control input value is *<sup>α</sup>f*, Δ*Mz*.

**Figure 5.** Flow chart of linear time-varying model control.

Since the linear time-varying model [13] is used, it is necessary to perform linear time-varying processing and discretization processing on Equation (24).

$$\text{Let } \mathbf{A}\_{\mathbf{k}} = \mathbf{I} + \mathbf{A}\_{\mathbf{T}\prime} \mathbf{B}\_{\mathbf{k}} = \mathbf{B} \mathbf{T}\prime \mathbf{C}\_{\mathbf{k}} = \mathbf{C} \dots$$

where T is the sampling time, I is the identity matrix,

$$\mathbf{A}\_{\mathbf{k}} = \begin{bmatrix} 1 - \frac{\mathbb{C}\_{cf} + \mathbb{C}\_{cr}}{m\upsilon\_{\mathbf{x}}}T & -\left(\frac{l\_f \mathbb{C}\_{cf} - l\_r \mathbb{C}\_{cr}}{m\upsilon\_{\mathbf{x}}} + \upsilon\_{\mathbf{x}}\right)T & 0 & 0\\ -\frac{l\_f \mathbb{C}\_{cf} - l\_r \mathbb{C}\_{cr}}{l\_z \upsilon\_{\mathbf{x}}}T & 1 - \frac{l\_f^2 \mathbb{C}\_{cf} + l\_r^2 \mathbb{C}\_{cr}}{l\_z \upsilon\_{\mathbf{x}}}T & 0 & 0\\ 0 & T & 1 & 0\\ T & 0 & \upsilon\_{\mathbf{x}}T & 1 \end{bmatrix} \tag{29}$$

$$B\_k = \begin{bmatrix} \frac{C\_{\varepsilon f}}{m} T & 0\\ \frac{I\_f \bar{C}\_{\varepsilon f}}{I\_z} T & \frac{T}{I\_z} \\ 0 & 0 \\ 0 & 0 \end{bmatrix} \tag{30}$$

$$\mathbf{C\_k} = \begin{bmatrix} 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{bmatrix} \tag{31}$$

So the following discrete state space equation is obtained:

$$\begin{cases} \dot{\mathfrak{F}}(k+1) = A\mathfrak{F}(k) + \mathfrak{B}u(k) \\ \mathfrak{y}(k) = \mathfrak{C}\_{k}\mathfrak{F}(k) \end{cases} \tag{32}$$

In order to avoid the phenomenon of a sudden change in the control quantity of the system, it is necessary to restrict the control increment in each sampling period, so the above formula of this paper is converted and deformed:

Note: *ζ* =, substitute into the above formula to obtain a new state space equation:

$$\begin{cases} \mathcal{J}(k+1) = \overline{A}\_k \mathcal{J}(k) + \overline{B}\_k \Delta u(k) \\ y(k) = \overline{\mathcal{C}}\_k \mathcal{J}(k) \end{cases} \tag{33}$$

In the formula,

$$
\tilde{A}\_k = \begin{bmatrix} A\_k & B\_k \\ 0\_{3 \times 4} & I\_2 \end{bmatrix} \tag{34}
$$

$$
\bar{B}\_k = \begin{bmatrix} B\_k \\ I\_2 \end{bmatrix} \tag{35}
$$

$$
\overline{\mathbf{C}}\_{k} = \begin{bmatrix} 0 & 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 \end{bmatrix} \tag{36}
$$

The optimization objective function of the path tracking and torque distribution controller is expressed as:

$$J = \sum\_{i=1}^{N\_P} \left\| y(k+i) - y\_{ref}(k+i) \right\|^2 Q + \sum\_{i=1}^{N\_\varepsilon - 1} \left\| \Delta u(k+i) \right\|^2 R \tag{37}$$

where *yref* is obtained by the path planning layer; **Q** and **R** are the trajectory tracking error weight matrix and the control amount increment weight matrix, respectively; *Np* and *Nc* are the prediction time domain and the control time domain, respectively. Some parameter settings in model predictive control are shown in Table 1.

**Table 1.** Part of the parameter settings of the controller.


#### *3.3. Torque Distribution Controller Design*

Torque distribution strategy has a grea<sup>t</sup> impact on vehicle lateral stability. Firstly, the additional yaw moment Δ*Mz* output by the MPC controller is used to calculate the total demand torque *Tall*. Then using the proportion of the vertical load, calculates the torque of the four wheels. On the premise of ensuring lateral stability, the flexible torque distribution is realized.

The additional yaw moment Δ*Mz* and the four-wheel torque should satisfy the following relationship:

$$
\Delta M\_z = \frac{\mathbf{b}}{2R} \left[ (T\_{Fr} - T\_{Fl}) \cos \delta\_f + (T\_{Rr} - T\_{Rl}) \right] \tag{38}
$$

where *R* is the wheel radius, *b* is the wheel track, and *TFl*, *TFr*, *TRl* and *TRr* are the torques of the left front wheel, the right front wheel, the left rear wheel and the right rear wheel, respectively.

$$F\_z = F\_{z,FI} + F\_{z,Fr} + F\_{z,RI} + F\_{z,Rr} \tag{39}$$

while:

$$\begin{array}{c} k\_1 = \frac{F\_{z,Fl}}{F\_z} \\ k\_2 = \frac{F\_{z,Fr}}{F\_z} \\ k\_3 = \frac{F\_{z,Fl}}{F\_z} \\ k\_4 = \frac{F\_{z,RF}}{F\_z} \end{array} \tag{40}$$

where *Fz,Fl*, *Fz,Fr*, *Fz,Rl*, *Fz,Rr* are the vertical loads of the left front wheel, the right front wheel, the left rear wheel and the right rear wheel, respectively. Force, *k*1, *k*2, *k*3, *k*4 are the proportional coefficients of the four-wheel vertical force and the total vertical force.

⎧⎪⎪⎪⎪⎨⎪⎪⎪⎪⎩

According to Equations (37)–(39), the total demand torque *Tall* and the distributed four-wheel torques *TFl*, *TFr*, *TRl*, *TRr* can be obtained:

$$T\_{all} = \frac{\Delta M\_z \cdot 2R}{d \cdot \left[ (k\_2 - k\_1) \cos \delta\_f + (k\_4 - k\_3) \right]} \tag{41}$$

$$\begin{cases} \begin{array}{c} T\_{Fl} = k\_1 T\_{all} \\ T\_{Fr} = k\_2 T\_{all} \\ T\_{Rl} = k\_3 T\_{all} \\ T\_{Rr} = k\_4 T\_{all} \end{array} \end{cases} \tag{42}$$

The solution of the final path tracking and torque distribution controller can be expressed as:

$$\begin{array}{ll}\min Z = \sum\_{i=1}^{N\_p} \left\lVert y(k+i) - y\_{ref}(k+i) \right\rVert^2 + \sum\_{i=1}^{N\_c-1} \left\lVert \Delta u(k+i) \right\rVert^2 \\ \text{s.t.} \; \zeta(k+1) = \overline{A}\_k \overline{\zeta}(k) + \overline{B}\_k \Delta u(k) \\ \Delta u\_{\min} \le \Delta u(k+i) \le \Delta u\_{\max} \\ u\_{\min} \le u(k+i) \le u\_{\max} \\ y\_{\min} \le y(k+i) \le y\_{\max} \\ y\_{\min} \le g(k+i) \le g\_{\max} \end{array} \tag{43}$$

The control quantity at time *k* + 1 can be expressed by the control quantity at time *k* plus the control increment at time *k*:

$$
\mathfrak{u}(k+1) = \mathfrak{u}(k) + \Delta\mathfrak{u}(k)\tag{44}
$$

#### **4. Simulation Analysis and Verification**

⎧⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎩

The C-class vehicle model is adopted in Carsim. PD18 in-wheel motor is adopted in the motor model. According to the characteristics of the distributed drive vehicle, the vehicle parameters are shown in Table 2.



#### *4.1. Verification of Co-Simulation Platform*

The initial speed of the vehicle is 60 km/h, the road adhesion coefficient *μ* = 0.8 and the maximum value of LTR is set at 0.1. The simulation results are shown in Figures 6 and 7.

It can be seen from Figure 5 that the controller can accurately track the reference path. The tracking error of the yaw angle is small, and the lateral acceleration appears slightly jittered between 0 and 1 s, but the value of the lateral acceleration is maintained throughout the process. Within a reasonable range, the front wheel corners shake slightly at 0 s to 1 s and 6 s to 7 s, but the changes in the corners are relatively gentle, and the corner values do not exceed the specified value. It can be concluded that the control effect of the controller is good and can meet the needs of vehicle obstacle avoidance.

**Figure 6.** Co-simulation results. (**a**) Path tracking. (**b**) Yaw angle tracking. (**c**) Lateral acceleration. (**d**) Front-wheel turning angle.

**Figure 7.** Changes in LTR value.

It can be seen from Figure 6 that the value of LTR is always in a small range during the entire obstacle avoidance process, which verifies the effectiveness of the set obstacle avoidance path and meets the requirements of anti-rollover in the obstacle avoidance process.

#### *4.2. Simulation of Tracking Control Performance*

The stability and robustness of the designed path tracking and torque distribution controller are verified by the simulation of different road adhesion coefficient conditions. Select high adhesion pavement coefficient *μ* = 0.8, low adhesion pavement coefficient *μ* = 0.2, *μ* = 0.1. The specific simulation results are shown in Figures 7 and 8.

It can be seen from Figure 8a–d that the tracking control effect of the controller is excellent when the high adhesion coefficient road surface *μ* = 0.8. The effect is slightly decreased when the low adhesion coefficient road surface *μ* = 0.2, and the tracking control error is slightly increased. The changes in lateral acceleration and front wheel angle are relatively gentle under both conditions. When *μ* = 0.1, the vehicle deviates greatly from the planned path and the planned yaw angle. The lateral acceleration and front wheel angle also produce a lot of jitter. This is because the road adhesion coefficient is too low to the adhesion limit. Under this condition, the controller cannot perform tracking control.

**Figure 8.** Simulation comparison under different road adhesion coefficient conditions. (**a**) Path tracking. (**b**) Yaw tracking. (**c**) Front-wheel angle. (**d**) Lateral acceleration.

It can be seen from Figure 9 that when *μ* = 0.8 and *μ* = 0.2, the vertical force changes of the wheels are relatively gentle, which meets the requirements of stable driving of the vehicle in the process of obstacle avoidance. When *μ* = 0.1, the tires are not sufficient to provide the vehicle with steady longitudinal and lateral force. At 2.5 s and 4.7 s, there is a large jump, which can no longer meet the stability requirements of the vehicle during driving.

**Figure 9.** Comparison of vertical forces on wheels. (**a**) vertical force of left wheel; (**b**) vertical force of right wheel.

Figure 10a–d shows that when *μ* = 0.8, the torque distribution of each wheel meets the requirements of this paper. When *μ* = 0.2, the torque of the left front wheel begins to increase after 1 s and reaches a peak at 1.5~3.5 s, the maximum value is 107 N·m. The torque starts to decrease continuously at 3.5~5 s and 5~5.5 s. It reaches the reverse maximum value of 116 N·<sup>m</sup> and then changes to 0 N·m. The torque changes of the remaining wheels are roughly the same as the torque changes of the left front wheel, so this conclusion description will not be repeated. It can be seen that under this working condition, the controller has a good torque distribution effect. When *μ* = 0.1, the torque of each wheel drops rapidly to about −100 N·<sup>m</sup> in 2 s, and the controller cannot continue to distribute the torque normally.

**Figure 10.** Comparison of torque distribution among wheels. (**a**) Torque comparison of left front wheel. (**b**) Torque comparison of left rear wheel. (**c**) Torque comparison of right front wheel. (**d**) Torque comparison of right rear wheel.

From the analysis of Figures 9 and 10, it can be seen that the controller designed in this paper can still control the vehicle to drive according to the correct planned path under the condition of low road adhesion coefficient *μ* = 0.2 and carry out reasonable torque distribution, which verifies the path. The tracking controller has good robustness. When *μ* = 0.1, the vehicle has already slipped, and the path tracking control can not be realized.
