**1. Introduction**

The rapid development of control-by-wire technology provides a strong technical guarantee for the realization of electronic, intelligent, and electrified vehicles, and makes x-by-wire vehicles a current research hotspot [1,2]. The x-by-wire electric vehicle is a highly redundant system, whose steer angle and driving torque of the four wheels can be precisely controlled; thus, compared with traditional vehicles, it can theoretically achieve better control effects [3]. The x-by-wire electric vehicle has a variety of driving modes that can be switched, with better trajectory tracking performance and a larger turning curvature limit, which can easily realize in situ steering, differential steering, oblique driving, and other working conditions, and can meet the requirements of higher precision trajectory tracking and critical working conditions [4–6]. The high degree of controllability and flexibility and high execution ability of x-by-wire electric vehicles provide a platform for the research into driverless vehicles, which has grea<sup>t</sup> development potential.

The problem of trajectory tracking control is one of the three key technologies related to intelligent driving vehicles [7]. As far as the trajectory tracking control of traditional front-wheel steered vehicles is concerned, the relevant theories and research methods are relatively mature, the relevant simulation experiments have been quite abundant, and

**Citation:** Wang, Z.; Li, Y.; Kaku, C.; Zheng, H. Trajectory Tracking Control of Intelligent X-by-Wire Vehicles. *World Electr. Veh. J.* **2022**, *13*, 205. https://doi.org/10.3390/ wevj13110205

Academic Editor: Joeri Van Mierlo

Received: 21 September 2022 Accepted: 26 October 2022 Published: 1 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/).

many control strategies have been verified by real vehicles in specific environments. The control strategies include the following: pure tracking control, traditional proportionalintegral-derivative (PID) control, sliding mode control, dynamic feedback control, fuzzy control, and model predictive control (MPC) [8–10]. Al-Mayyahi et al. [11] proposed a PID-based fractional-order (FOPID) controller for the trajectory tracking problem, using two FOPID controllers to calculate the control input torque of the vehicle and using a particle swarm optimization algorithm to control the parameters in the FOPID controller. Kapania et al. [12] added the vehicle's side-slip characteristics to the feedforward controller to improve the dynamic characteristics of the tire within the friction limit while keeping the vehicle within the handling limit, minimizing lateral tracking deviation. Wu [13] proposed a back stepping sliding mode controller to reduce the chattering problem of the sliding mode control. Using the fuzzy control algorithm for the overdrive system, the simulation experiment proves that the back stepping sliding mode control method has higher control precision and a smoother control process than the traditional sliding mode control. Mallem et al. [14] proposed a fast terminal sliding mode dynamic inverse control method based on PID, which makes the position and direction of the movement around the desired trajectory asymptotically stable, and used a terminal sliding mode control method to ensure finite time convergence of the trajectory tracking error to zero. Rokonuzzaman et al. [15] proposed to use the large amount of data provided by the vehicle to design MPC with a neural-network-based vehicle learning dynamic model to improve the tracking performance; the results under various road conditions show that the proposed method outperforms the MPC of the traditional vehicle model. Funke et al. [16] proposed a new control structure based on model predictive control and feedback control that integrates path tracking, vehicle stabilization, and collision avoidance, and coordinated these conflicting goals through the priority of collision avoidance. Experimental data show that the controller drives safely within the operational limits of the vehicle.

The x-by-wire electric vehicle has the advantages of independent four-wheel drive/brake and steering control, so it has been widely studied by universities and enterprises. However, the research on trajectory tracking control mainly focuses on traditional vehicles and there are few studies on trajectory tracking control algorithms for x-by-wire electric vehicles.

Based on proportional-integral (PI) control and adaptive model predictive control, Ahn et al. [17] designed an integrated autonomous driving system independent of each wheel for vehicles equipped with four-wheel independent motors to improve vehicle stability and path tracking performance. Li et al. [7] proposed a sliding mode drive controller based on PID control and sliding mode control for 4WIS/4WID vehicles. Hiraoka et al. [18] proposed an automatic controller for four-wheel steering vehicle path tracking based on sliding mode control theory. Compared with active front-wheel steering, the four-wheel steering controller has a more stable and accurate trajectory tracking capability. Zheng et al. [19] designed a trajectory tracking strategy based on a hierarchical control method. The path tracking layer adopts a nonlinear state feedback controller, and a neural network proportional-integral-derivative controller is designed to track the desired path and obtain the desired yaw rate. Chen et al. [20] designed a new adaptive linear quadratic optimal regulator (LQR) as a coordinated controller for 4WIS/4WID electric vehicle stability control. According to different vehicle speeds and road conditions, the phase plane method is used to calculate the center of mass slip angle and stability margin.

On the other hand, the current research on the trajectory tracking control of x-by-wire electric vehicles rarely considers constraints such as tire adhesion. As each wheel of the x-by-wire vehicle can be independently controlled, a reasonable control algorithm and control strategy can achieve independent longitudinal and lateral control of each wheel, improving the accuracy of trajectory tracking and the driving stability of the vehicle.

This paper takes the x-by-wire electric vehicle as the research carrier and designs the trajectory tracking control strategy based on the hierarchical control architecture, which includes the trajectory tracking layer, the tire force distribution layer, and the actuator control layer. The main contribution of this paper is as follows:


The rest of this article is structured as follows. Section 2 introduces a nonlinear x-bywire chassis vehicle model. Section 3 introduces the trajectory tracking controller based on a hierarchical architecture. Section 4 introduces the simulation under the DLC test. Section 5 presents the conclusion of the article.

#### **2. Vehicle Dynamic Model and Tire Model**

#### *2.1. Vehicle Dynamic Model*

The nonlinear dynamic model of the vehicle established in this paper is as follows:

$$\begin{cases} \dot{\upsilon}\_{y} = -\upsilon\_{x}r + \frac{1}{m}\sum F\_{y} \\ \dot{\upsilon}\_{x} = \upsilon\_{y}r + \frac{1}{m}\sum F\_{x} \\ \dot{r} = \frac{1}{T\_{\natural}}\sum M\_{z} \\ \dot{\varepsilon}\_{\theta} = r - \dot{s}\kappa\_{r} \\ \dot{\varepsilon}\_{l} = \upsilon\_{x}\sin\varepsilon\_{\theta} + \upsilon\_{y}\cos\varepsilon\_{\theta} \\ \dot{s} = \frac{\upsilon\_{x}\cos\varepsilon\_{\theta} - \upsilon\_{y}\sin\varepsilon\_{\theta}}{1 - l\kappa\_{r}} \end{cases} \tag{1}$$

where *m* is the mass of the vehicle. *Iz* is the moment of inertia around the center of mass. *s* is the distance of the desired path. *κr* is the curvature of the desired path. *r*, *vx*, and *vy* are the yaw rate, longitudinal speed, and lateral speed of the vehicle, respectively. *el* and *<sup>e</sup>ϕ* are the lateral deviation and the heading angle deviation between the center of mass of the vehicle and the reference waypoint, respectively.∑ *Fx*, ∑ *Fy*, and ∑ *Mz* are the lateral force, longitudinal force, and yaw moment received by the center of mass of the vehicle, respectively.

The above vehicle model is nonlinear. Considering that the heading deviation *<sup>e</sup>ϕ* is generally small, it is assumed that cos*<sup>e</sup>ϕ* ≈ 1, sin*<sup>e</sup>ϕ* ≈ 0. At the same time, considering that the road curvature and lateral deviation are generally small, the lateral speed of the vehicle is generally much smaller than the longitudinal speed. Based on this, the following can be obtained:

$$\begin{cases} \dot{s} \approx \upsilon\_{\text{x}} - \upsilon\_{y}\varepsilon\_{\text{\(\sigma\)}} \approx \upsilon\_{\text{x}}\\ \dot{\varepsilon}\_{l} \approx \upsilon\_{\text{x}}\varepsilon\_{\text{\(\sigma\)}} + \upsilon\_{y} \\ \dot{\upsilon}\_{\text{\(\sigma\)}} = \frac{1}{m}\sum F\_{\text{x}} \end{cases} \tag{2}$$

The state equation of the vehicle lateral model is as follows:

$$
\frac{d}{dt} \begin{bmatrix} v\_y \\ w \\ e\_\varphi \\ e\_l \end{bmatrix} = \begin{bmatrix} 0 & -v\_x & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 1 & 0 & v\_x & 0 \end{bmatrix} \begin{bmatrix} v\_y \\ w \\ e\_\varphi \\ e\_l \end{bmatrix} + \begin{bmatrix} \frac{1}{m} & 0 \\ 0 & \frac{1}{l\_z} \\ 0 & 0 \\ 0 & 0 \end{bmatrix} \begin{bmatrix} \sum F\_y \\ \sum M\_z \end{bmatrix} + \begin{bmatrix} 0 \\ 0 \\ -v\_x \\ 0 \end{bmatrix} \kappa\_r \tag{3}
$$

The state equation of the vehicle longitudinal model is as follows:

$$
\frac{d}{dt} \begin{bmatrix} s \\ v\_x \end{bmatrix} = \begin{bmatrix} 0 & 1 \\ 0 & 0 \end{bmatrix} \begin{bmatrix} s \\ v\_x \end{bmatrix} + \frac{1}{m} \sum F\_x \tag{4}
$$

#### *2.2. Tire Model*

At present, there are three types of tire models: the physical tire model, empirical– semi-empirical tire model, and finite element tire model. In this paper, the Pacejka tire model is selected, which is a semi-empirical tire model and can be expressed as follows [21]:

$$Y(\mathbf{x}) = D \sin \{ \text{Cortanh}[B\mathbf{x} - E(B\mathbf{x} - \text{arctan}(B\mathbf{x}))] \}$$

where *<sup>Y</sup>*(*x*) is the lateral tire force. *x* is the slip angle of tire. *B* is the stiffness factor. *C*, *D*, *E* are curve shape factor, peak factor, and curvature factor, respectively.

#### **3. Trajectory Tracking Controller Based on Hierarchical Architecture**

The hierarchical control algorithm structure has the advantages of clear algorithm design and being convenient for subsequent algorithm updating. For the x-by-wire chassis vehicle, which is a research object with high integration and high controllable degrees of freedom, the method of controlling each subsystem separately cannot give full play to its own performance advantages. The control strategy of the hierarchical structure can realize the coordinated control of each controller of the vehicle with an x-by-wire chassis. As shown in Figure 1, in the hierarchical trajectory tracking controller, it is divided into the trajectory tracking layer, tire force distribution layer, and actuator control layer. In the trajectory tracking layer, which is based on the MPC, the expected total force and moment are calculated by taking the longitudinal force ∑ *Fx*, lateral force ∑ *Fy*, and yaw moment ∑ *Mz* in the vehicle coordinate system as the control variables. The tire force distribution layer solves for the longitudinal force *Fl*,*i*,*<sup>j</sup>* and lateral force *Fc*,*i*,*<sup>j</sup>* in the vehicle coordinate system of each wheel. The actuator control layer obtains the tire steer angle *δij* and driving torque *Tij* based on the inverse tire force model and controls the controlled vehicle to track the trajectory.

**Figure 1.** Trajectory tracking controller architecture.

#### *3.1. Trajectory Tracking Layer*

#### 3.1.1. Prediction Model

In this paper, the trajectory tracking layer is established based on the MPC algorithm and the state equation established based on Equation (1) is used as the prediction model as follows:

$$
\frac{d}{dt} \begin{bmatrix} v\_{\mathcal{Y}} \\ r \\ e\_{\mathcal{P}} \\ e\_{\mathcal{I}} \\ s \\ v\_{\mathcal{X}} \end{bmatrix} = \begin{bmatrix} 0 & -v\_{\mathcal{X}} & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 & 0 \\ 1 & 0 & v\_{\mathcal{X}} & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 \\ 0 & 0 & 0 & 0 & 0 & 0 \end{bmatrix} \begin{bmatrix} v\_{\mathcal{Y}} \\ r \\ e\_{\mathcal{P}} \\ e\_{\mathcal{I}} \\ s \\ v\_{\mathcal{X}} \end{bmatrix} + \begin{bmatrix} \frac{1}{m} & 0 & 0 \\ 0 & \frac{1}{l\_{\mathcal{X}}} & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & \frac{1}{m} \end{bmatrix} \begin{bmatrix} \sum F\_{\mathcal{Y}} \\ \sum M\_{\mathcal{Z}} \\ \sum F\_{\mathcal{X}} \\ \sum F\_{\mathcal{X}} \end{bmatrix} + \begin{bmatrix} 0 \\ 0 \\ -v\_{\mathcal{X}} \\ 0 \\ 0 \\ 0 \end{bmatrix} \kappa\_{r} \tag{5}
$$

Equation (5) is rewritten as follows:

$$\begin{cases} \dot{\mathbf{x}}(t) = A\_c \mathbf{x}(t) + B\_{c\mu} \mathbf{u}(t) + B\_{cd} d(t) \\ y(t) = \mathbb{C}\_c \mathbf{x}(t) \end{cases} \tag{6}$$

*Ac* is the state matrix, *Bcu* is the control matrix, *Bcd* is the disturbance matrix, and *Cc* is the output matrix, respectively:

$$A\_c = \begin{bmatrix} 0 & -\upsilon\_x & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 & 0 \\ 1 & 0 & \upsilon\_x & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 \\ 0 & 0 & 0 & 0 & 0 & 0 \end{bmatrix}, B\_{cu} = \begin{bmatrix} \frac{1}{m} & 0 & 0 \\ 0 & \frac{1}{T\_x} & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & \frac{1}{m} \end{bmatrix}, C\_c = I\_{6 \times 6}$$

Discretization of Equation (5):

$$\begin{cases} \boldsymbol{\chi}(k+1) = A\_d \boldsymbol{\chi}(k) + B\_{du} \boldsymbol{\mu}(k) + B\_{dd} d(k) \\ \boldsymbol{\chi}(t) = \mathbb{C} \boldsymbol{\chi}(t) \end{cases} \tag{7}$$

The discrete sampling time is defined as *Ts*, then the discretized coefficient matrix is as follows:

$$\begin{array}{l} A\_{d} = \mathfrak{e}^{A\_{\mathfrak{c}}T\_{\mathfrak{s}}} \\ B\_{d\mathfrak{u}} = B\_{\mathfrak{c}\mathfrak{u}} \int\_{0}^{T\_{\mathfrak{s}}} \mathfrak{e}^{A\_{\mathfrak{c}}\tau} d\tau \\ B\_{d d} = B\_{\mathfrak{c}d} \int\_{0}^{T\_{\mathfrak{s}}} \mathfrak{e}^{A\_{\mathfrak{c}}\tau} d\tau \\ C\_{d} = C\_{\mathfrak{c}} \end{array} \tag{8}$$

In order to reduce the vibration of the control amount, the control amount is rewritten as an incremental type:

$$
\Delta u(k) = u(k) - u(k-1) \tag{9}
$$

The state equation is as follows:

$$\mathbf{x}(k+1) = A\_d \mathbf{x}(k) + B\_{du} (\mu(k-1) + \Delta u(k)) + B\_{dd} d(k) \tag{10}$$

The augmented state is as follows:

$$\mathfrak{F}(k) = \begin{bmatrix} \mathfrak{x}(k) \\ \mathfrak{u}(k-1) \end{bmatrix} \tag{11}$$

Then, the equation of state can be rewritten as follows:

$$\begin{cases} \xi(k+1) = \overline{A}\_d \overline{\xi}(k) + \overline{B}\_{du} \Delta u(k) + \overline{B}\_{dd} d(k) \\ y(k) = \overline{\mathbb{C}}\_d \overline{\xi}(k) \end{cases} \tag{12}$$

Among:

$$\begin{aligned} \;^c \overline{A}\_d = \begin{bmatrix} A\_d & B\_{du} \\ O & I \end{bmatrix}, \; \overline{B}\_{du} = \begin{bmatrix} B\_{du} \\ I \end{bmatrix}, \; \overline{B}\_{dd} = \begin{bmatrix} B\_{dd} \\ O \end{bmatrix}, \; \overline{C}\_d = \begin{bmatrix} C\_d & O \end{bmatrix}. \end{aligned}$$

After applying the control variables, the model prediction process is as follows:

$$
\zeta(k+1) = \overline{A}\_d \zeta(k) + \overline{B}\_{du} \Delta \mu(k) + \overline{B}\_{dd} d(k) \tag{13}
$$

$$\begin{array}{ll} \mathfrak{F}(k+2) &= \overline{A}\_{d}\mathfrak{F}(k+1) + \overline{B}\_{du}\Delta u(k+1) + \overline{B}\_{dd}d(k+1) \\ &= \overline{A}\_{d}^{2}\mathfrak{F}(k) + \overline{A}\_{d}\overline{B}\_{du}\Delta u(k) + \overline{B}\_{du}\Delta u(k+1) + \overline{A}\_{d}\overline{B}\_{dd}d(k) + \overline{B}\_{dd}d(k+1) \end{array} \tag{14}$$

$$\begin{aligned} \mathfrak{F}(k+N\_{\mathsf{C}}) &= \overline{A}\_{d}\mathfrak{F}(k+N\_{\mathsf{C}}-1) + \overline{B}\_{d\mathsf{u}}\Delta u(k+N\_{\mathsf{C}}-1) + \overline{B}\_{d\mathsf{d}}d(k+N\_{\mathsf{C}}-1) \\ &= \overline{A}\_{d}^{N\_{\mathsf{C}}}\mathfrak{F}(k) + \sum\_{i=0}^{N\_{\mathsf{C}}-1} \overline{A}\_{d}^{i}\overline{B}\_{d\mathsf{u}}\Delta u(k+N\_{\mathsf{c}}-1-i) + \sum\_{i=0}^{N\_{\mathsf{c}}-1} \overline{A}\_{d}^{i}\overline{B}\_{d\mathsf{d}}d(k+N\_{\mathsf{c}}-1-i) \end{aligned} \tag{15}$$

$$\begin{split} \tilde{\xi}\left(\mathbf{k} + \mathbf{N\_{p}}\right) &= \overline{A}\_{d}\tilde{\xi}\left(\mathbf{k} + \mathbf{N\_{p}} - 1\right) + \overline{B}\_{\text{dir}}\Delta\mu\left(\mathbf{k} + \mathbf{N\_{c}} - 1\right) + \overline{B}\_{\text{dir}}d\left(\mathbf{k} + \mathbf{N\_{p}} - 1\right) \\ &= \overline{A}\_{d}^{\mathbf{N\_{p}}}\tilde{\xi}\left(\mathbf{k}\right) + \sum\_{i=0}^{\mathbf{N\_{c}}-1} \overline{A}\_{d}^{i + \mathbf{N\_{p}} - \mathbf{N\_{c}}} \overline{B}\_{\text{dir}}\Delta\mu\left(\mathbf{k} + \mathbf{N\_{c}} - 1 - i\right) + \sum\_{i=0}^{\mathbf{N\_{p}}-1} \overline{A}\_{d}^{i} \overline{B}\_{\text{dd}}d\left(\mathbf{k} + \mathbf{N\_{p}} - 1 - i\right) \end{split} \tag{16}$$

Definition:

$$Y(k+1) = \begin{bmatrix} y(k+1) \\ y(k+2) \\ \vdots \\ y(k+N\_{\mathcal{C}}) \\ \vdots \\ y(k+N\_{\mathcal{P}}) \end{bmatrix} = \begin{bmatrix} \overline{\mathsf{C}}\_{d}\mathfrak{s}(k+1) \\ \overline{\mathsf{C}}\_{d}\mathfrak{s}(k+2) \\ \vdots \\ \overline{\mathsf{C}}\_{d}\mathfrak{s}(k+N\_{\mathcal{C}}) \\ \vdots \\ \overline{\mathsf{C}}\_{d}\mathfrak{s}(k+N\_{\mathcal{P}}) \end{bmatrix}, \\ X(k) = \begin{bmatrix} \mathfrak{s}(k+1) \\ \mathfrak{s}(k+2) \\ \vdots \\ \mathfrak{s}(k+N\_{\mathcal{C}}) \\ \vdots \\ \mathfrak{s}(k+N\_{\mathcal{P}}) \end{bmatrix}, \Delta D(k) = \begin{bmatrix} d(k) \\ d(k+1) \\ \vdots \\ d(k+N\_{\mathcal{C}}) \\ \vdots \\ d(k+N\_{\mathcal{P}}) \end{bmatrix}$$

where *Y*(*k* + 1)is the output vector, *X*(*k*)is the state vector, and Δ*D*(*k*)is the disturbance vector.

$$Y(k+1) = S\_{\mathbf{x}}\mathcal{Z}(k) + S\_{\mathbf{u}}\Delta l I(k) + S\_{d}D(k) \tag{17}$$

The state matrix *Sx*, control matrix *Su*, and disturbance matrix *Sd* are, in respective order, as follows:

*Sx* = ⎡⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣ *CdAd CdA*<sup>2</sup>*d* ... *CdANc d* ... *CdANp d* ⎤⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦*Su* = ⎡⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣ *CdBdu O* ··· *O CdAdBdu CdBdu* ··· *O* ... ... ... ... *CdANc*−<sup>1</sup> *d Bdu CdANc*−<sup>2</sup> *d Bdu* ··· *CdBdu* ... ... ... ... *CdANp*−<sup>1</sup> *d Bdu CdANp*−<sup>2</sup> *d Bdu* ··· *CdANp*−*Nc d Bdu*⎤⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦*Sd* = ⎡⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣ *CdBdd O* ··· *O CdAdBdd CdBdd* ··· *O* ... ... ... ... *CdANc*−<sup>1</sup> *d Bdd CdANc*−<sup>2</sup> *d Bdd* ··· *O* ... ... ... ... *CdANp*−<sup>1</sup> *d Bdd CdANp*−<sup>2</sup> *d Bdd* ··· *CdBdd*⎤⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦

3.1.2. Cost Function

In order to ensure the safe and accurate driving of the x-by-wire electric vehicle according to the predetermined trajectory, the cost function is as follows:

$$\begin{split} \mathcal{I} &= \sum\_{i=1}^{N\_{\rm tr}} \left[ \Gamma\_{\rm ri} \left( y(k+i) - y\_{\rm ref} \right) \right]^2 + \sum\_{i=0}^{N\_{\rm tr}-1} \left[ \Gamma\_{\rm ir} (\Delta u(k+i)) \right]^2 \\ &= \left( Y(k+i) - Y\_{\rm ref} \right)^T \Gamma\_{\rm r} \left( Y(k+i) - Y\_{\rm ref} \right) + \Delta \mathcal{U}(k)^T \Gamma\_{\rm n} \Delta \mathcal{U}(k) \\ &= \left( \mathcal{S}\_{\rm r} \mathcal{I}(k) + \mathcal{S}\_{\rm n} \Delta \mathcal{U}(k) + \mathcal{S}\_{\rm d} \mathcal{D}(k) - Y\_{\rm ref} \right)^T \Gamma\_{\rm r} \left( \mathcal{S}\_{\rm r} \mathcal{I}(k) + \mathcal{S}\_{\rm u} \Delta \mathcal{U}(k) + \mathcal{S}\_{\rm d} \mathcal{D}(k) - Y\_{\rm ref} \right) + \Delta \mathcal{U}(k)^T \Gamma\_{\rm u} \Delta \mathcal{U}(k) \end{split} \tag{18}$$

The error *E* is defined as follows:

$$E = S\_x \zeta(k) + S\_d D(k) - \chi\_{ref} \tag{19}$$

Equation (18) be rewritten as follows:

$$\begin{array}{rcl} J &=& \left( S\_{\mathsf{u}} \Delta \mathsf{U} I(k) + E(k) \right)^{T} \Gamma\_{\mathsf{x}} \left( S\_{\mathsf{u}} \Delta \mathsf{U} I(k) + E(k) \right) + \Delta \mathsf{U} I(k)^{T} \Gamma\_{\mathsf{u}} \Delta \mathsf{U} I(k) \\ &=& \Delta \mathsf{U}(k)^{T} \left( \Gamma\_{\mathsf{u}} + S\_{\mathsf{u}}{}^{T} \Gamma\_{\mathsf{x}} S\_{\mathsf{u}} \right) \Delta \mathsf{U} I(k) + 2 \Delta \mathsf{U} (k)^{T} S\_{\mathsf{u}}{}^{T} \Gamma\_{\mathsf{x}} E(k) + E(k)^{T} \Gamma\_{\mathsf{x}} E(k) \end{array} \tag{20}$$

The symmetric matrix is as follows:

$$\Gamma\_x = \operatorname{diag} \left( \Gamma\_{x1\nu} \cdot \cdot \cdot \Gamma\_{xN\_p} \right) \\ \Gamma\_u = \operatorname{diag} \left( \Gamma\_{u1\nu} \cdot \cdot \cdot \Gamma\_{uN\_c} \right)$$

To simplify further,

$$\begin{array}{lcl}H &= \Gamma\_{\boldsymbol{\mu}} + \mathcal{S}\_{\boldsymbol{\mu}}{}^{T}\Gamma\_{\boldsymbol{X}}\mathcal{S}\_{\boldsymbol{\mu}}\\f &= 2\mathcal{S}\_{\boldsymbol{\mu}}{}^{T}\Gamma\_{\boldsymbol{X}}E(k) + E(k)^{T}\Gamma\_{\boldsymbol{X}}E(k)\end{array} \tag{21}$$

Substituting Equation (21) into Equation (18), the cost function is as follows:

$$J = \Delta l I(k)^T H \Delta l I(k) + \Delta l I(k)^T f + E(k)^T \Gamma\_x E(k) \tag{22}$$

#### 3.1.3. Security Constraints

When the vehicle is driving, it is necessary to ensure the stability and safety of the vehicle. Vehicle stability characteristics can be captured by nonlinear tire state models. Bobier [22] uses the center of mass slip angle-yaw rate to handle the vehicle's stability constraint envelope.

When the additional yaw moment generated by the longitudinal force is not considered, in the steady state, the relationship between the yaw rate of the vehicle and the lateral force is as follows:

.

.

*r*

$$
\dot{r} = \frac{\sum F\_y}{m v\_x} \tag{23}
$$

As the road friction coefficient limit is ∑ *Fy* ≤ *μmg*, the limit of the yaw rate at which the vehicle is stable is as follows:

$$t \le \frac{\mu \text{g}}{v\_{\text{x}}} \tag{24}$$

The safe driving environment of the vehicle is the safe envelope area considering road boundaries, obstacles, and traffic vehicles. The obstacle avoidance and stability control framework proposed by Erlien et al. [23] uses the safe driving envelope to delineate a collision-free area for the vehicle to exercise and the environmental safety constraints are constrained by road boundaries and traffic vehicles. This paper is based on the road boundary of the tracked trajectory as a constraint, which can be expressed as follows:

$$H\_{\rm env} \mathfrak{F} \le \begin{bmatrix} \mathfrak{e}\_{\rm max} - \frac{B}{2} - l\_{\rm buffer} \\ -\mathfrak{e}\_{\rm min} + \frac{B}{2} + l\_{\rm buffer} \end{bmatrix} \tag{25}$$

*Henv* = 000 1 00 000 −<sup>100</sup> . *e*max and *e*min represent the upper and lower boundariesoftheroad,respectively.*lbuffer*istheboundarysafetymargin.

#### *3.2. Tire Force Distribution Layer*

As the four wheels of an x-by-wire electric vehicle can be independently controlled, the control variables obtained based on the trajectory tracking layer are the longitudinal force ∑ *Fx*, lateral force ∑ *Fy*, and yaw moment ∑ *Mz* at the center of mass. Therefore, the resultant force and torque need to be distributed to each wheel. The tire force distribution layer adopts an optimized method to distribute the force/moment to each wheel under the condition of satisfying the tire force constraint.

#### 3.2.1. Tire Force Constraints

During the driving of the vehicle, the resultant force generated by each tire must satisfy the constraints of the trajectory tracking layered force equation:

$$\begin{cases} \begin{aligned} \mathbf{F}\_{\mathbf{x},fl} + \mathbf{F}\_{\mathbf{x},fr} + \mathbf{F}\_{\mathbf{x},rl} + \mathbf{F}\_{\mathbf{x},rr} &= \sum \mathbf{F}\_{\mathbf{x}} \\ \mathbf{F}\_{\mathbf{y},fl} + \mathbf{F}\_{\mathbf{y},fr} + \mathbf{F}\_{\mathbf{y},rl} + \mathbf{F}\_{\mathbf{y},rr} &= \sum \mathbf{F}\_{\mathbf{y}} \\ \left(\mathbf{F}\_{\mathbf{y},fl} + \mathbf{F}\_{\mathbf{y},fr}\right)l\_{f} - \left(\mathbf{F}\_{\mathbf{y},rl} + \mathbf{F}\_{\mathbf{y},rr}\right)l\_{r} + \frac{\mathbf{B}}{2}\left(-\mathbf{F}\_{\mathbf{x},fl} + \mathbf{F}\_{\mathbf{x},fr} - \mathbf{F}\_{\mathbf{x},rl} + \mathbf{F}\_{\mathbf{x},rr}\right) &= \sum \mathbf{M}\_{z} \end{aligned} \end{cases} \tag{26}$$

When distributing the tire force, in order to ensure the stability of the vehicle, it is necessary to consider the adhesion capacity of the wheel and tire. Assuming that the road friction coefficient is *μ*, the maximum resultant force that the tire can generate is *F* ≤ *μFz*, where *Fz* is the vertical load.

$$F\_x^2 + F\_y^2 \le \left(\mu F\_z\right)^2\tag{27}$$

As shown in Figure 2, the longitudinal force and lateral force of the tire must satisfy the friction circle constraint. When the wheel is in the extreme condition, the tire force will reach the limit of the friction circle and the wheel will slip, resulting in the loss of vehicle stability. The constraint of tire adhesion is nonlinear and the tire force problem can be defined as an optimization problem with quadratic constraints, but it is often time-consuming. In order to improve the real-time performance of the algorithm, the tire force constraint is simplified on the premise of ensuring the driving stability of the vehicle. Therefore, this paper adopts the circumscribed regular octagon of the friction circle to approximately describe the tire friction constraint and converts the quadratic constraint into a linear constraint.

**Figure 2.** Tire friction constraint.

As shown in Figure 2, *R* is the radius of friction circle and *Re* is the radius of the simplified octagon.

·

$$R\_{\varepsilon} = R \cdot \sec 22.5^{\circ} \approx 1.08R$$

$$\frac{S\_{\text{octagon}}}{S\_{circle}} = \frac{8}{\pi \left(1 + \sqrt{2}\right)} \approx 1.0548$$

According to the geometric calculation, the radius of the regular octagon is 1.08 times the radius of the friction circle and the area of the octagon is about 5.5% larger than that of the friction circle. This linearization treatment has little effect on tire force and ensures the driving stability of an x-by-wire chassis vehicle. The regular octagon is replaced by the friction circle and the octagon linear inequality constraint of the tire friction circle is

shown in Equation (28). *Fx*, *Fy* must be inside the regular octagon *A*1*A*2*A*3*A*4*A*5*A*6*A*7*A*8 and the constraints of the attached ellipse can be expressed as follows:

$$\begin{cases} -\mu F\_z \le F\_x \le \mu F\_z \\ -\mu F\_z \le F\_y \le \mu F\_z \\ -\sqrt{2}\mu F\_z \le F\_x + F\_y \le \sqrt{2}\mu F\_z \\ -\sqrt{2}\mu F\_z \le F\_x - F\_y \le \sqrt{2}\mu F\_z \end{cases} \tag{28}$$

For the four tires of the vehicle, the matrix form is expressed as follows:

$$F = \begin{bmatrix} F\_{\mathbf{x},fl} & F\_{\mathbf{x},fr} & F\_{\mathbf{x},rl} & F\_{\mathbf{x},rr} & F\_{\mathbf{y},fl} & F\_{\mathbf{y},fr} & F\_{\mathbf{y},rl} & F\_{\mathbf{y},rr} \end{bmatrix}^T \tag{29}$$

#### 3.2.2. Objective Function

The tire adhesion margin is defined as follows [24]:

$$\varepsilon = 1 - \frac{F}{\mu F\_z} \tag{30}$$

The adhesion margin of the tire represents the ratio of the remaining utilization force of the tire to the maximum force provided by the tire. The value range is 0 ∼ 1 and *ε* = 0 represents the tire reaching the adhesion limit. Thus, the objective function is defined to maximize the attachment margin. The tire usage rate is defined as follows:

$$\eta = 1 - \varepsilon = \frac{F}{\mu F\_z} \tag{31}$$

Based on Equation (31), the objective function is defined as the minimum sum of the tire usage rates of the four tires, which is as follows:

$$J = \sum\_{i=1}^{4} \frac{F\_{xi}^2 + F\_{yi}^2}{\left(\mu F\_{zi}\right)^2} \tag{32}$$

In this article, the main object of coordination control is the lateral and longitudinal tire forces, which are related to vertical force. However, the vertical forces are not the control object of the control strategy, so the vertical force is regarded as directly available through sensors or other means.

#### 3.2.3. Tire Force Distribution Algorithm

From Equation (32), the tire force distribution problem is regarded as an optimization problem with constraints:

$$\begin{array}{l}\underset{\{F\_{x,ij},F\_{y,ij}\}}{\min} \\ \text{s.t.} \quad Ax \le b \end{array} \tag{33}$$

Based on the objective function in Equation (32) and the optimization variable in Equation (29), the optimization problem is a quadratic programming problem, which can be solved quickly.

#### *3.3. Actuator Control Layer*

X-by-wire electric vehicles follow a desired trajectory by having independent drive/brake and steering control of each wheel. The actuator control layer converts the longitudinal force and lateral force of each wheel under the vehicle coordinate system obtained by the tire force distribution layer into the tire coordinate system. Figure 3 presents the tire force diagram.

**Figure 3.** Tire force of a single wheel.

The tire force relationship is as follows:

$$
\begin{bmatrix} F\_{c, ij} \\ F\_{l, ij} \end{bmatrix} = \begin{bmatrix} \cos \delta\_{ij} & \sin \delta\_{ij} \\ -\sin \delta\_{ij} & \cos \delta\_{ij} \end{bmatrix} \begin{bmatrix} F\_{x, ij} \\ F\_{y, ij} \end{bmatrix} \tag{34}
$$

*δij* is the wheel steer angle, *<sup>α</sup>ij* is the tire slip angle, and *θij* is the angle between the driving direction of the wheel and the longitudinal axis of the wheel coordinate system.

$$
\theta\_{i\bar{j}} = \delta\_{i\bar{j}} + \alpha\_{i\bar{j}} \tag{35}
$$

The tire side slip angle *<sup>α</sup>ij* is generally small. Equation (35) is approximated as follows:

$$\theta\_{ij} \approx \delta\_{ij}$$

The force in the tire coordinate system can be obtained as follows:

$$
\begin{bmatrix}
\hat{F}\_{\mathbf{c},ij} \\
\hat{F}\_{l,ij}
\end{bmatrix} = \begin{bmatrix}
\cos\theta\_{ij} & \sin\theta\_{ij} \\
\end{bmatrix} \begin{bmatrix}
F\_{\mathbf{x},ij} \\
F\_{y,jj}
\end{bmatrix} \tag{36}
$$

The lateral and longitudinal velocity of each wheel are as follows:

$$\begin{cases} v\_{y,fl} = v\_{y,fr} = v\_y + l\_f \dot{\phi} \\ v\_{y,rl} = v\_{y,rr} = v\_y - l\_f \dot{\phi} \\ v\_{x,fl} = v\_{x,fr} = v\_x - \frac{R}{2} \dot{\phi} \\ v\_{x,fr} = v\_{x,rr} = v\_x + \frac{R}{2} \dot{\phi} \end{cases} \tag{37}$$

It can be obtained from Equation (37) that the angle between the motion direction of each wheel and the longitudinal direction of the body coordinate system is as follows:

$$\begin{cases} \theta\_{fl} = \tan^{-1} \frac{\overline{v}\_{y,fl}}{\overline{v}\_{x,fl}} = \tan^{-1} \frac{\overline{v}\_{y} + l\_f \dot{\varphi}}{\overline{v}\_{x} - \frac{\dot{\theta}}{\dot{\theta}} \dot{\varphi}}\\ \theta\_{fr} = \tan^{-1} \frac{\overline{v}\_{y,fr}}{\overline{v}\_{x,fr}} = \tan^{-1} \frac{\overline{v}\_{y} + l\_f \dot{\varphi}}{\overline{v}\_{x} + \frac{\dot{\theta}}{\dot{\theta}} \dot{\varphi}}\\ \theta\_{rl} = \tan^{-1} \frac{\overline{v}\_{y,rl}}{\overline{v}\_{x,rl}} = \tan^{-1} \frac{\overline{v}\_{y} - l\_f \dot{\varphi}}{\overline{v}\_{x} - \frac{\dot{\theta}}{\dot{\theta}} \dot{\varphi}}\\ \theta\_{rr} = \tan^{-1} \frac{\overline{v}\_{y,rr}}{\overline{v}\_{x,rr}} = \tan^{-1} \frac{\overline{v}\_{y} + l\_f \dot{\varphi}}{\overline{v}\_{x} + \frac{\dot{\theta}}{\dot{\theta}} \dot{\varphi}} \end{cases} \tag{38}$$

The actuator control layer solves the relationship between tire slip angle and lateral force based on the arctangent model proposed by Sakai et al. [25], the function of which could fit the magic tire formula. The tire side slip model can be expressed as follows:

$$\begin{cases} F\_y = -C G\_x \frac{\mu}{k} \tan^{-1} \left( \frac{k}{\mu} a \right) \\\ G\_x = \sqrt{1 - \left( \frac{F\_y}{\mu F\_z} \right)^2} \\\ k = \frac{C \pi}{2F\_z} \end{cases} \tag{39}$$

*k*and *Gx* are factors and *C* is the tire side slip stiffness.

The desired *F l*,*ij* is brought into Equation (39) and the tire slip angle *α*ˆ*ij* can be obtained, then the wheel angle is as follows:

$$
\delta\_{\rm ij} = \theta\_{\rm ij} - \hbar\_{\rm ij} \tag{40}
$$

The longitudinal moment of the wheel is as follows:

$$T\_{i\dot{j}} = \frac{1}{1 + \tau s} F\_{i\dot{j}} r\_{i\dot{j}} \tag{41}$$

*τ*is the time constant and *rij* is the radius of rotation of the wheel.

#### **4. Simulation Test**

ˆ

On the one hand, real vehicle experiments have a certain risk owing to a variety of conditions. On the other hand, in order to ensure the vehicle in the experiment has the expected trajectory and speed and to reduce the influence of the experimenter's subjective control on the experimental results, focus is placed upon the influence of the control strategy on the vehicle driving stability in the trajectory tracking control. In this paper, the effectiveness of the control strategy is verified by simulation experiment.

Through MATLAB/Simulink and CarSim simulation, the hierarchical trajectory tracking control algorithm is verified by the double line change (DLC) test. The vehicle parameters are shown in Table 1.


**Table 1.** Vehicle parameters.

#### *4.1. Variable Velocity DLC Condition*

On a road with a friction coefficient of 0.85, the trajectory tracking of an x-by-wire electric vehicle under variable speed is simulated. The initial velocity of the vehicle is 20 km/h and the simulation time is set to 20 s.

As shown in Figure 4, the control algorithm can follow the changing vehicle velocity well and accurately track the desired trajectory; the lateral deviation and heading angle deviation of the trajectory tracking are kept within a small range.

**Figure 4.** Vehicle state and tracking error: (**a**) velocity; (**b**) vehicle location; (**c**) lateral error; (**d**) heading error.

As shown in Figure 5, the values of the vehicle's center of mass slip angle and yaw angular velocity all change within a small range and the changes are relatively stable; the change range in the vehicle's lateral speed and lateral acceleration is also small, which proves the vehicle has lateral stability.

As shown in Figure 5a, the driving torque has a jitter before 5 s, which is due to the large deviation between the real vehicle velocity and the reference velocity at this time. In order to reduce this deviation, the driving torque is increased. When the deviation is reduced to a reasonable range, the driving torque returns to the normal value. The change in velocity is shown in Figure 4a.

**Figure 5.** Vehicle state: (**a**) driving torque; (**b**) steer angle; (**c**) side slip angle; (**d**) yaw rate; (**e**) lateral velocity; (**f**) lateral acceleration.

As shown in Figure 6, the trajectory tracking control algorithm can reasonably distribute the load of the four tires while following the desired trajectory to ensure the stability of the vehicle.

**Figure 6.** Tire state: (**a**) tire vertical load; (**b**)tire utilization.

#### *4.2. Low Road Friction Coefficient DLC Condition*

In order to further verify that the control algorithm has good following ability and stability, a pavement with friction coefficient of 0.35 was selected for the DLC test. The initial speed is 40 km/h.

As shown in Figure 7, the control algorithm can also follow the desired vehicle speed well and accurately track the desired trajectory in the low friction coefficient road environment.

**Figure 7.** Vehicle state and tracking error: (**a**) velocity; (**b**) vehicle location; (**c**) lateral error; (**d**) heading error.

As shown in Figure 8, the lateral deviation and heading angle deviation of trajectory tracking are kept within a small range. The change in the control amount is relatively stable and there is no major fluctuation. The side slip angle and yaw rate also kept fluctuating within a relatively stable range.

**Figure 8.** Vehicle state: (**a**) driving torque; (**b**) steer angle; (**c**) side slip angle; (**d**) yaw rate; (**e**) lateral velocity; (**f**) lateral acceleration.

As shown in Figure 9, it can be seen from the tire utilization coefficient of each wheel that the vehicle runs stably and safely on a low-adhesion road surface. It is proved that the algorithm proposed in this paper is suitable for various working conditions. Based on the curve of tire force and tire utilization coefficient in the simulation experiment, it conforms to the change in tire performance of each wheel in the real driving process of the vehicle. On the other hand, the vehicle state parameters also conform to the actual driving performance of the vehicle. Therefore, the results of the simulation experiment are considered to be reliable.

**Figure 9.** Tire state: (**a**) tire vertical load; (**b**) tire utilization.
