Path Tracking Control of Autonomous Vehicle Based on Nonlinear Tire Model

Abstract

The tire forces of vehicles will fall into the non-linear region under extreme handling conditions, which cause poor path tracking performance. In this paper, a model predictive controller based on a nonlinear tire model is designed. The tire forces are characterized with nonlinear composite functions of the magic formula instead of a simple linear relation model. Taylor expansion is used to linearize the controller, the first-order difference quotient method is used for discretization, and the partial derivative of the composite function is used for matrix transformation. Constant velocity and variable velocity conditions are selected to compare the designed controller with the conventional controller in Carsim/Simulink. The results show that when the tire forces fall in the nonlinear region, two controllers have good stability, and the tracking accuracy of the controller designed in this paper is slightly better. However, after the tire forces become nonlinear, the controller with linear tire force becomes worse, the tracking accuracy is far worse than the controller with the nonlinear tire model, and the vehicle stability is also degraded. In addition, an active steering test platform based on LabVIEW-RT is established, and hardware-in-the-loop tests are carried out. The effectiveness of the designed controller is verified.

1. Introduction

Autonomous vehicles are the focus of attention in the world today. Many car manufacturers and even some Internet companies have invested a lot of energy in the research and development of autonomous vehicles. As early as 2008, Urmson C believed that intelligent systems could be applied to the automotive industry. He introduced the advantages of autonomous vehicles and the challenges faced by cities in the direction of autonomous driving [1]. Then, Lee J analyzed the needs of autonomous vehicles, designed an agent that can simulate autonomous vehicles, and developed a prototype for it [2]. Obinata G believed that it is necessary to solve the problems of machine control, road safety and automation technology reliability of autonomous vehicles, and explained the current situation and challenges of unmanned operation technology [3]. More and more researchers are advancing the development of autonomous cars from different angles and within different fields. From the perspective of data storage, Yeshodara NS proposed a new idea to reduce the data storage problems of autonomous vehicles, only relying on cloud infrastructure to drive cars [4]. From the perspective of obstacle recognition, Hane C proposed a method of extracting static obstacles from depth maps calculated from multiple consecutive images [5]. Kumar GA proposed a method to estimate the distance between autonomous vehicles and other vehicles, objects and signs using precise fusion methods [6]. From the perspective of computer vision perception, Dongpu Cao used a direct perception method to train a deep neural network in terms of increasing learning ability in relation to autonomous driving [7]. From the perspective of path planning, Li CC proposed a model-based dynamic environment path planning algorithm [8], and Jinkai Y proposed a hybrid path planning algorithm based on a simulated annealing algorithm and a particle swarm algorithm [9]. Different research angles have brought hope to the real realization of autonomous cars to some extent; however, they must be implemented with consideration of the control of the vehicle, so path tracking control will eventually grow in importance.

In the control of autonomous cars, Jain AK installed a Raspberry Pi and cameras on the top of the car. An image was sent to the convolutional neural network for manipulation prediction, and then the prediction signal was sent to the controller to make the vehicle move in the intended direction without human intervention [10]. Chen H proposed a new path-following control method and designed a three-step controller with output feedback to realize the coordinated control of the lateral and longitudinal directions without measuring the lateral speed [11]. Li L proposed a high-performance automatic steering control strategy, using the system state equation of direction angle and lateral offset deviation to describe the tracking accuracy, and the control method adopted the sliding mode variable structure control method [12]. Zhang H proposed the use of the Takagi–Sugeno model to deal with the time-varying characteristics of vehicle speed, as well as the use of Lyapunov stability parameters to improve the transient stability of the vehicle system [13]. Celentano L proposed a simple and robust smooth path tracking control method that allows unmanned electric vehicles to track a smooth reference signal [14]. Hain H proposed the use of a fuzzy controller to collect driver data while controlling an off-road vehicle. The data were used to develop the membership function of the input and output of the fuzzy controller, thereby improving the autonomous vehicle’s path tracking ability [15]. Taking into account the conflict between tracking accuracy and stability under extreme conditions, Liang YX proposed a new path tracking controller based on yaw rate tracking [16]. Hang P proposed a control algorithm to establish a linear parameter change system model to achieve different longitudinal speeds and different road friction coefficients, and combined the use of feedforward control with linear quadratic regulators to increase tracking stability [17]. There are various vehicle motion control methods, but they are all aimed at improving control accuracy and stability. Among these control methods, model predictive control takes place earlier, and it uses non-minimized description models to improve the robustness of the system. It uses a rolling optimization strategy to make up for the uncertainty and has better dynamic performance [18], and the algorithm has the ability to deal with constraints systematically, so it is widely used to solve the path tracking problems of autonomous vehicles.

Falcone P first proposed a new automatic steering system method [19]. The designed model predictive controller emphasizes the trade-off between vehicle speed and required road preview to stabilize the vehicle. Then, Falcone P performed continuous online linearization of the complex nonlinear vehicle dynamics model [20], and also proposed a low-complexity linear time-varying model, which considers the state of the vehicle and input constraints to improve vehicle stability on slippery roads and at high speeds, and used a variety of optimization algorithms to extend the range of conditions for real-time MPC implementation. According to the literature [19,20], Falcone proposed a predictive model for the joint control of braking and steering of autonomous vehicles. On high-speed snow-covered roads, the vehicle can be better controlled to change lanes on two lanes [21]. For the problem of model mismatch caused by different road conditions and small angles of vehicle models, Guo HY adopted the form of measurable disturbance, and used differential evolution algorithm for the optimization problem of path tracker [22]. Yin GD designed a decentralized MPC controller, which considers the rear-end collision avoidance of the control area and the collision avoidance at intersections to improve traffic efficiency and safety [23]. He ZW proposed a two-layer controller for accurate and robust lateral path tracking control of highly automated vehicles. The upper controller is realized by assuming linear time-varying MPC at the small angle, and the slip angle is limited to ensure stability. The lower controller is a radial basis function neural network PID controller [24]. Peng HN proposed a robust model predictive control method with a finite time domain to realize the path tracking and direct yaw moment control of an electric vehicle independently driven by an autonomous four-wheel motor [25]. Based on the analysis of the sideslip angle, Tang LQ proposed a vehicle sideslip angle compensator that can correct the motion model prediction, which significantly improves the path tracking performance and the control ability at high speed [26]. Xu SB proposed the preview steering control algorithm and the concept of its closed-loop system analysis, and compared the model method with the MPC model based on the bicycle model [27]. Yao proposed an MPC model with longitudinal velocity compensation in the prediction range, and analyzed the longitudinal velocity change mechanism and the stability within the prediction range [28]. Sara M introduced a Tube-based MPC to consider the dynamic difference between the real vehicle and this constant nominal model, not only considering the lateral error, but also the directional error [29]. Michael J extended constraint removal to Tube-based MPC, detecting inactive constraints in advance, deleting invalid constraints, and reducing calculation time [30]. Shilp D used a combination of potential fields such as vehicle functions and reachability sets to identify safe zones on the road, and provided them to the Tube-based MPC, which can be free from non-convex collision avoidance constraints and ensure the feasibility of lateral movement during acceleration and deceleration [31]. Among the above MPC controller design processes, most chose the bicycle model to simplify the calculation of the algorithm, and used the small angle assumption to represent the tire force in the vehicle model. When the tire force enters the nonlinear region, errors will inevitably occur, which will affect the performance of the controller. Therefore, an MPC algorithm based on a nonlinear tire model is proposed: Replace the conventional tire force linearly expressed by cornering stiffness with the tire force expressed by the magic formula composite function, so that even if the tire force enters the nonlinear region, the tire force calculated by the magic formula is close to the tire force under real conditions. Avoid unnecessary errors caused by tire force in path tracking.

The main arrangement of this work is as follows. First, in the first section, a bicycle model and a non-linear magic formula tire model are established. In the second section, the complex partial derivative processing and matrix transformation related to the nonlinear tire model, as well as the prediction model, objective function and constraint conditions of the controller are introduced in detail. In the third section, a variety of constant working conditions and variable working conditions are selected, and the effectiveness of the controller is verified based on Carsim/Simulink simulation. In the fourth section, an active steering test platform based on LabVIEW-RT is established, and hardware-in-the-loop testing is carried out to verify the actual control effect of the controller. The fifth section describes the conclusion.

Only the replacement of the linearly expressed tire force with the tire force expressed by a composite function is considered to verify the performance of the controller, so as to promote the improvement and research of the controller on this basis in the future.

2. Vehicle Dynamics Model

2.1. Bicycle Model

Establishing a suitable vehicle model is the prerequisite and basis for realizing model predictive control. Since the main research goal is to make the car follow the desired path quickly and stably, the vehicle model needs to be simplified to a certain extent, and the bicycle model is used to transfer left and right without load, assuming that the vehicle has a good anti-lock braking system.

The plane motion of the vehicle is mainly considered as lateral, longitudinal and yaw motion. The vehicle is set as front-wheel drive and front-wheel steering, and the vehicle model is selected as a three-degree-of-freedom bicycle model [32]. In the plane rectangular coordinate system, the vehicle model is shown in Figure 1.

According to Newton’s second law, the force balance equations along the x-axis, y-axis and z-axis can be obtained.

In the x-axis direction:

$$m\ddot{x}=2F_{xf}+2F_{xr}+m\dot{y}\dot{\phi }$$

(1)

In the y-axis direction:

$$m\ddot{y}=2F_{yf}+2F_{yr}-m\dot{x}\dot{\phi }$$

(2)

Around the z axis:

$$I_z\ddot{\phi }=2F_{yf}\cdot l_f-2F_{yr}\cdot l_r$$

(3)

where $m$ is the curb weight of the vehicle, $\dot{x}$ and $\dot{y}$ are the longitudinal and lateral speeds, $F_{xf}$ and $F_{yf}$ are the forces in the x and y directions received by the front tires, $F_{xr}$ and $F_{yr}$ are the forces in the x and y directions received by the rear tires, $\dot{\phi }$ is the yaw rate, $I_z$ is the moment of inertia of the vehicle around the z-axis, and $l_f$ and $l_r$ are the distances from the center of mass to the front and rear axles.

The force received by the front and rear tires in the x and y directions is related to the longitudinal force and lateral force received by the front and rear tires. The relationship is as follows [32]:

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

(4)

In the formula, $F_{lf}$ and $F_{lr}$ are the longitudinal forces received by the front and rear wheels, $F_{cf}$ and $F_{cr}$ are the lateral forces received by the front and rear wheels, ${\delta }_f$ and ${\delta }_r$ are the front and rear steering angles. Since the selected vehicle model is front-wheel steering, ${\delta }_r=0$.

The longitudinal force and lateral force of the tire can be expressed as a complex function of the tire side slip angle $\alpha$, tire slip rate $s$, vertical load $F_z$ and other parameters:

$$\begin{array}{l}{F}_l=f_l(\alpha ,s,F_z)\\ F_c=f_c(\alpha ,s,F_z)\end{array}$$

(5)

Among them, the side slip angle of the tire is related to the longitudinal speed and the lateral speed, as shown in Equation (6) [32]:

$$\begin{array}{l}{\alpha }_f=\arctan \frac{{\nu }_{cf}}{{\nu }_{lf}}\\ {\alpha }_r=\arctan \frac{{\nu }_{cr}}{{\nu }_{lr}}\end{array}$$

(6)

In the formula, ${\alpha }_f$ and ${\alpha }_r$ are the front wheel slip angle and the rear wheel slip angle, ${\nu }_{lf}$ and ${\nu }_{lr}$ are the longitudinal speeds of the front and rear wheels, and ${\nu }_{cf}$ and ${\nu }_{cr}$ are the lateral speeds of the front and rear wheels.

The longitudinal and lateral speed of each tire can be represented by the speed in the x direction, the speed in the y direction, and the front and rear steering angles [33]:

$$\begin{array}{l}{\nu }_{lf}={\nu }_{yf}\sin {\delta }_f+v_{xf}\cos {\delta }_f\\ {\nu }_{lr}={\nu }_{yr}\sin {\delta }_r+v_{xr}\cos {\delta }_r\\ {\nu }_{cf}={\nu }_{yf}\cos {\delta }_f-v_{xf}\sin {\delta }_f\\ {\nu }_{cr}={\nu }_{yr}\cos {\delta }_r-v_{xr}\sin {\delta }_r\end{array}$$

(7)

In the formula, ${\nu }_{xf}$ and ${\nu }_{xr}$ are the speeds of the front and rear wheels in the x direction, and ${\nu }_{yf}$ and ${\nu }_{yr}$ are the speeds of the front and rear wheels in the y direction.

The speed of the tires in the x and y directions is generally difficult to obtain, and can be calculated from the vehicle speed [33]:

$$\begin{array}{l}{\nu }_{yf}=\dot{y}+l_f\dot{\phi }\\ {\nu }_{yr}=\dot{y}-l_r\dot{\phi }\\ {\nu }_{xf}=\dot{x}\\ {\nu }_{xr}=\dot{x}\end{array}$$

(8)

Therefore, the front and rear wheel slip angles can be obtained from Equations (7) and (8).

Assuming that the vehicle has a good anti-lock braking system, the adhesion coefficient reaches the maximum when the tire slip rate $s$ is 15–20%. In order to achieve the best braking effect, the slip rate needs to be controlled at 15–20% Therefore, the slip rate is a constant value of 20%.

Since the load transfer of the vehicle is not considered, assuming that the vehicle speed changes slowly and there is no load transfer of the front and rear axles, the vertical loads $F_{zf}$ and $F_{zr}$ of the front and rear wheels can be expressed as [33]:

$$\begin{array}{l}{F}_{zf}=\frac{l_{r}mg}{2\left(l_f+l_r\right)}\\ F_{zr}=\frac{l_{f}mg}{2\left(l_f+l_r\right)}\end{array}$$

(9)

The above formulas are derived from the body coordinate system, and the relationship conversion between the body coordinate system and the inertial coordinate system needs to be considered [33]:

$$\begin{array}{l}\dot{Y}=\dot{x}\sin \phi +\dot{y}\cos \phi \\ \dot{X}=\dot{x}\cos \phi -\dot{y}\sin \phi \end{array}$$

(10)

In the formula, $\phi$ is the yaw angle, $X$ is the longitudinal position in the inertial coordinate system, and $Y$ is the horizontal position in the inertial coordinate system.

2.2. Magic Formula Tire Model

In order to reduce the amount of calculation in most of the current path tracking algorithm models, the small angle assumption is used in Equation (5) to approximate the longitudinal force and lateral force of the tire. When the tire slip angle and longitudinal slip rate are small, the tire force can be linearly approximated as:

$$F_l=C_{l}s,F_c=C_c\alpha$$

(11)

In the formula, $C_l$ is the tire longitudinal stiffness, and $C_c$ is the tire cornering stiffness.

The use of Equation (11) has a certain limit range, and only when the lateral acceleration $a_y$ is less than 0.4 g can it have a higher fitting accuracy for conventional tires. If you consider that at a high speed or when the road adhesion coefficient is low, when the lateral acceleration of the vehicle is greater than 0.4 g at a certain moment, the tire longitudinal force and lateral force at that moment are substituted into the control algorithm model, and there will be a certain error in the calculation result. Because the increase in tire force at this time becomes slower with the increase in slip rate and cornering angle, the tire force calculated with the small-angle linear assumption is obviously greater than the tire force under real conditions. In order to make the model more accurate, the path tracking algorithm model is based on the magic formula tire model [34], replacing the small angle hypothesis of Equation (11) with the magic formula.

When the vehicle is running, the camber angle is generally small. The influence of the camber angle on the tire longitudinal force can be ignored, and the empirical parameters related to the camber in the magic formula can be ignored. Ignore the parameters related to zero-point drift, ignore the empirical parameters that cause asymmetry in the formula, and ignore the higher-order terms more than 3 times; thus, the simplified magic formula is as follows:

In the purely longitudinal case, the longitudinal force is [34]:

$$F_{x0}=D_x\sin \left[C_x\arctan \left\{B_x{s}_x-E_x\left(B_x{s}_x-\arctan \left(B_x{s}_x\right)\right)\right\}\right]+S_{vx}$$

(12)

In the formula:

$$\begin{array}{l}{s}_x=s+S_{Hx},S_{Hx}=P_{Hx1}+P_{Hx2}d{f}_z,C_x=P_{Cx1},D_x=\left(P_{Dx1}+P_{Dx2}\cdot d{f}_z\right)\cdot F_z,\\ E_x=P_{Ex1}+P_{Ex2}\cdot d{f}_z+P_{Ex3}\cdot d{f}_z{}^2,S_{vx}=F_z\cdot \left(P_{Vx1}+P_{Vx2}\cdot d{f}_z\right),d{f}_z=\frac{F_Z-F_{Z0}}{F_{Z0}},\\ B_x=\left(P_{Kx1}+P_{Kx2}\cdot d{f}_z\right)\cdot e^{P_{Kx3}\cdot d{f}_z}/\left[P_{Cx1}\cdot \left(P_{Dx1}+P_{Dx2}\cdot d{f}_z\right)\right].\end{array}$$

In the case of pure sideslip, the lateral force is [34]:

$$F_{y0}=D_y\sin \left[C_y\arctan \left\{B_y{\alpha }_y-E_y\left(B_y{\alpha }_y-\arctan \left(B_y{\alpha }_y\right)\right)\right\}\right]+S_{vy}$$

(13)

In the formula:

$$\begin{array}{l}{\alpha }_y=\alpha +S_{Hy},S_{Hy}=P_{Hy1}+P_{Hy2}d{f}_z,C_y=P_{Cy1},D_y=\left(P_{Dy1}+P_{Dy2}\cdot d{f}_z\right)/\left(1+P_{Dy3}\cdot {\gamma }^{*}\right)\cdot F_z,{\gamma }^{*}\approx 1^{°},\\ E_y=P_{Ey1}+P_{Ey2}\cdot d{f}_z,S_{vy}=F_z\cdot \left(P_{Vy1}+P_{Vy2}\cdot d{f}_z\right),d{f}_z=\frac{F_Z-F_{Z0}}{F_{Z0}},\\ B_y=P_{Ky1}\cdot F_{z0}\cdot \sin \left[2\arctan \left(\frac{F_z}{P_{Ky2}\cdot F_{z0}}\right)\right]/\left[P_{Cy1}\cdot \left(P_{Dy1}+P_{Dy2}\cdot d{f}_z\right)\cdot F_z\right].\end{array}$$

In the mixed case of longitudinal sideslip [34]:

The longitudinal force is:

$$F_x=G_{x\alpha }\cdot F_{x0}$$

(14)

$$G_{x\alpha }=\frac{\cos \left[C_{x\alpha }\arctan \left\{B_{x\alpha }{\alpha }_s-E_{x\alpha }\left(B_{x\alpha }{\alpha }_s-\arctan \left(B_{x\alpha }{\alpha }_s\right)\right)\right\}\right]}{\cos \left[C_{x\alpha }\arctan \left\{B_{x\alpha }{S}_{Hx\alpha }-E_{x\alpha }\left(B_{x\alpha }{S}_{Hx\alpha }-\arctan \left(B_{x\alpha }{S}_{Hx\alpha }\right)\right)\right\}\right]}$$

(15)

In the formula:

$$C_{x\alpha }=r_{Cx1},E_{x\alpha }=r_{Ex1}+r_{Ex2}\cdot d{f}_z,{\alpha }_s=\alpha +S_{Hx\alpha },S_{Hx\alpha }=r_{Hx1},B_{x\alpha }=r_{Bx1}\cdot \cos \left(\arctan \left(r_{Bx2}\cdot s\right)\right).$$

The lateral force is:

$$F_y=G_{ys}\cdot F_{y0}+S_{vys}$$

(16)

$$G_{ys}=\frac{\cos \left[C_{ys}\arctan \left\{B_{ys}{s}_s-E_{ys}\left(B_{ys}{s}_s-\arctan \left(B_{ys}{s}_s\right)\right)\right\}\right]}{\cos \left[C_{x\alpha }\arctan \left\{B_{ys}{S}_{Hys}-E_{x\alpha }\left(B_{ys}{S}_{Hys}-\arctan \left(B_{ys}{S}_{Hys}\right)\right)\right\}\right]}$$

(17)

In the formula:

$$\begin{array}{l}{s}_s=s+S_{Hys},S_{Hys}=r_{Hy1}+r_{Hy2}\cdot d{f}_z,C_{ys}=r_{Cy1},E_{ys}=r_{Ey1}+r_{Ey2}\cdot d{f}_z,B_{ys}=r_{By1}\cdot \cos \left(\arctan \left\{r_{By2}\left(\alpha -r_{By3}\right)\right\}\right),\\ S_{vys}={\mu }_y\cdot F_z\cdot \left(r_{vy1}+r_{vy2}\cdot d{f}_z\right)\cdot \cos \left[\arctan \left(r_{vy4}\alpha \right)\right]\cdot \sin \left[r_{vy5}\arctan \left(r_{vy6}\cdot s\right)\right],\\ {\mu }_y=\left(P_{Dy1}+P_{Dy2}\cdot d{f}_z\right)/\left(1+P_{Dy3}\cdot r^{{*}^2}\right), r^{*}\approx 1^{°}.\end{array}$$

Combining Equations (12)–(17) with the parameters in Table 1, the longitudinal force and lateral force of the front and rear wheels can be calculated. Equation (5) can be expressed as follows:

$$\begin{array}{l}{F}_{lf}=4724\cos \left[1.125\arctan \left\{4.47\left({\alpha }_f-0.03\right)-0.0015\left(4.47\left({\alpha }_f-0.03\right)-\arctan \left(4.47\left({\alpha }_f-0.03\right)\right)\right)\right\}\right]\\ F_{lr}=4018\cos \left[1.125\arctan \left\{4.47\left({\alpha }_r-0.03\right)-0.0007\left(4.47\left({\alpha }_r-0.03\right)-\arctan \left(4.47\left({\alpha }_r-0.03\right)\right)\right)\right\}\right]\\ F_{cf}=\frac{\cos \left[1.1\arctan \left\{0.2036B_{ys1}-0.2792\left(0.2036B_{ys1}-\arctan \left(0.2036B_{ys1}\right)\right)\right\}\right]}{\cos \left[1.1\arctan \left\{0.0036B_{ys1}-0.2792\left(0.0036B_{ys1}-\arctan \left(0.0036B_{ys1}\right)\right)\right\}\right]}\cdot F_{yof}\\ F_{cr}=\frac{\cos \left[1.1\arctan \left\{0.1993B_{ys2}-0.2054\left(0.1993B_{ys2}-\arctan \left(0.1993B_{ys2}\right)\right)\right\}\right]}{\cos \left[1.1\arctan \left\{-0.00074B_{ys2}-0.2054\left(-0.00074B_{ys2}-\arctan \left(-0.00074B_{ys2}\right)\right)\right\}\right]}\cdot F_{yor}\end{array}$$

(18)

In the formula:

$$\begin{array}{l}{B}_{ys1}=6.38\cos \left[\arctan \left\{7.95\left({\alpha }_f+0.06\right)\right\}\right],B_{ys2}=6.38\cos \left[\arctan \left\{7.95\left({\alpha }_r+0.06\right)\right\}\right],\\ F_{yof}=-4037\sin \left[1.29\arctan \left\{9.33\left({\alpha }_f+0.00314\right)+0.988\left(\begin{array}{l}9.33\left({\alpha }_f+0.00314\right)\\ -\arctan \left(9.33\left({\alpha }_f+0.00314\right)\right)\end{array}\right)\right\}\right]+190,\\ F_{yor}=-3513\sin \left[1.29\arctan \left\{9.87\left({\alpha }_r+0.00368\right)+1.111\left(\begin{array}{l}9.87\left({\alpha }_r+0.00368\right)\\ -\arctan \left(9.87\left({\alpha }_r+0.00368\right)\right)\end{array}\right)\right\}\right]+180\end{array}$$

It can be seen from Equation (18) that the longitudinal force and lateral force of front and rear tires are only related to their respective slip angles. By substituting the above equation into Equation (4), the forces received by the front and rear tires in the x direction and the y direction are:

$$\begin{array}{l}{F}_{xf}=F_{lf}\cos {\delta }_f-F_{cf}\sin {\delta }_f\\ F_{xr}=F_{lr}\\ F_{yf}=F_{lf}\sin {\delta }_f+F_{cf}\cos {\delta }_f\\ F_{yr}=F_{cr}\end{array}$$

(19)

Substituting Equation (19) into Equations (1)–(3), and combining Equation (10), we can obtain the vehicle dynamics nonlinear model based on the magic formula tire model:

$$\{\begin{array}{l}m\ddot{y}=-m\dot{x}\dot{\phi }+2F_{lf}\sin {\delta }_f+2F_{cf}\cos {\delta }_f+2F_{cr}\\ m\ddot{x}=m\dot{y}\dot{\phi }+2F_{lf}\cos {\delta }_f-2F_{cf}\sin {\delta }_f+2F_{lr}\\ \dot{\phi }=\dot{\phi }\\ I_z\ddot{\phi }=2l_f{F}_{lf}\sin {\delta }_f+2l_f{F}_{cf}\cos {\delta }_f-2l_r{F}_{cr}\\ \dot{Y}=\dot{x}\sin \phi +\dot{y}\cos \phi \\ \dot{X}=\dot{x}\cos \phi -\dot{y}\sin \phi \end{array}$$

(20)

The selected tire type is 175/70 R13 (asymmetric). The specific parameters of the tire based on the magic formula are as follows:

Table 1. 175/70 R13 (asymmetric) tire parameters.

Parameter Name	Parameter Value	Parameter Name	Parameter Value
$P_{Cx1}$	1.62	$P_{Cy1}$	1.29
$P_{Dx1}$	1.035	$P_{Dy1}$	−0.9
$P_{Dx2}$	−0.0487	$P_{Dy2}$	0.18
$P_{Ex1}$	0.5	$P_{Dy3}$	−4.5
$P_{Ex2}$	−0.122	$P_{Ey1}$	−1.07
$P_{Ex3}$	−0.0063	$P_{Ey2}$	0.68
$P_{Kx1}$	19.4	$P_{Ky1}$	−12.95
$P_{Kx2}$	−0.13	$P_{Ky2}$	1.72
$P_{Kx3}$	0.171	$P_{Hy1}$	0.0035
$P_{Hx1}$	−0.0005	$P_{Hy2}$	−0.003
$P_{Hx2}$	8.42 × 10−5	$P_{Vy1}$	0.0045
$P_{Vx1}$	0	$P_{Vy2}$	−0.03
$P_{Vx2}$	0	$F_{z0}$	4100
$r_{Cx1}$	1.125	$r_{Cy1}$	1.1
$r_{Ex1}$	0.078	$r_{Ey1}$	0.23
$r_{Ex2}$	−0.16	$r_{Ey2}$	0.41
$r_{Hx1}$	−0.03	$r_{By1}$	6.38
$r_{Bx1}$	9	$r_{By2}$	7.95
$r_{Bx2}$	−8.75	$r_{By3}$	−0.06
$r_{Hy1}$	0.0007	$r_{Vy4}$	10
$r_{Hy2}$	0.024	$r_{Vy5}$	1.95
$r_{Vy1}$	0	$r_{Vy6}$	−50
$r_{Vy2}$	0

Table 1. 175/70 R13 (asymmetric) tire parameters.

Parameter Name	Parameter Value	Parameter Name	Parameter Value
$P_{Cx1}$	1.62	$P_{Cy1}$	1.29
$P_{Dx1}$	1.035	$P_{Dy1}$	−0.9
$P_{Dx2}$	−0.0487	$P_{Dy2}$	0.18
$P_{Ex1}$	0.5	$P_{Dy3}$	−4.5
$P_{Ex2}$	−0.122	$P_{Ey1}$	−1.07
$P_{Ex3}$	−0.0063	$P_{Ey2}$	0.68
$P_{Kx1}$	19.4	$P_{Ky1}$	−12.95
$P_{Kx2}$	−0.13	$P_{Ky2}$	1.72
$P_{Kx3}$	0.171	$P_{Hy1}$	0.0035
$P_{Hx1}$	−0.0005	$P_{Hy2}$	−0.003
$P_{Hx2}$	8.42 × 10−5	$P_{Vy1}$	0.0045
$P_{Vx1}$	0	$P_{Vy2}$	−0.03
$P_{Vx2}$	0	$F_{z0}$	4100
$r_{Cx1}$	1.125	$r_{Cy1}$	1.1
$r_{Ex1}$	0.078	$r_{Ey1}$	0.23
$r_{Ex2}$	−0.16	$r_{Ey2}$	0.41
$r_{Hx1}$	−0.03	$r_{By1}$	6.38
$r_{Bx1}$	9	$r_{By2}$	7.95
$r_{Bx2}$	−8.75	$r_{By3}$	−0.06
$r_{Hy1}$	0.0007	$r_{Vy4}$	10
$r_{Hy2}$	0.024	$r_{Vy5}$	1.95
$r_{Vy1}$	0	$r_{Vy6}$	−50
$r_{Vy2}$	0

3. Model Predictive Controller Design

3.1. Prediction Model

Convert Equation (20) in vehicle dynamics model into state space form [33]. Select lateral velocity $\dot{y}$, longitudinal velocity $\dot{x}$, yaw angle $\phi$, yaw rate $\dot{\phi }$, lateral offset $Y$, and longitudinal position $X$ as the state variable $\mathit{x}={\left[\dot{y},\dot{x},\phi ,\dot{\phi },Y,X\right]}^T$. In Equation (20), $F_{lf},F_{cf},F_{lr},F_{cr}$ are related to ${\alpha }_f,{\alpha }_r$, and ${\alpha }_f,{\alpha }_r$ are related to the front steering angle ${\delta }_f$ from Equations (6)–(8). Therefore, in the control system, the control variable is selected as $\mathit{u}={\delta }_f$. The general form of the state equation is as follows:

$$\dot{\mathit{x}}=f\left(\mathit{x},\mathit{u}\right)$$

(21)

The general vehicle path tracking control will give a reference path to realize the control tracking of the actual vehicle. The state variable and control variable of the reference vehicle at each time satisfy the above equation, and $d$ represents the reference variable, so the form becomes:

$${\dot{x}}_d=f\left(x_d,u_d\right)$$

(22)

In the formula: $x_d={\left[{\dot{y}}_d,{\dot{x}}_d,{\phi }_d,{\dot{\phi }}_d,Y_d,X_d\right]}^T,u_d={\delta }_{fd}$.

Perform a Taylor expansion of Equation (21) at any reference trajectory point, ignoring higher-order terms, and obtain [33]:

$$\dot{x}=f\left(x_d,u_d\right)+{\frac{\partial f\left(x,u\right)}{\partial x}|}_{\begin{array}{l}x=x_d\\ u=u_d\end{array}}\left(x-x_d\right)+{\frac{\partial f\left(x,u\right)}{\partial u}|}_{\begin{array}{l}x=x_d\\ u=u_d\end{array}}\left(u-u_d\right)$$

(23)

Subtract Equations (22) and (23) to obtain:

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

(24)

In the formula, the $A$ matrix is the first order partial derivative matrix of $f$ to $x$, and the $B$ matrix is the first order partial derivative matrix of $f$ to $u$. $\tilde{x}=x-x_d$, $\tilde{u}=u-u_d$.

Equation (24) is the state equation, but the state equation is continuous. The model predictive controller is a discrete-time control method, so Equation (24) needs to be discretized [33]:

$$\begin{array}{l}A\left(k\right)=I+TA\\ B\left(k\right)=TB\end{array}$$

(25)

In the formula, $T$ is the sampling time, and $I$ is the $6\times 6$ unit matrix.

The new discretized state space model is [33]:

$$\tilde{x}\left(k+1\right)=A\left(k\right)\tilde{x}\left(k\right)+B\left(k\right)\tilde{u}\left(k\right)$$

(26)

In the formula:

$$\begin{array}{l}A\left(k\right)=\left[\begin{array}{cccccc}1+T\cdot \frac{\partial \ddot{y}}{\partial \dot{y}}& T\cdot \frac{\partial \ddot{y}}{\partial \dot{x}}& 0& T\cdot \frac{\partial \ddot{y}}{\partial \dot{\phi }}& 0& 0\\ T\cdot \frac{\partial \ddot{x}}{\partial \dot{y}}& 1+T\cdot \frac{\partial \ddot{x}}{\partial \dot{x}}& 0& T\cdot \frac{\partial \ddot{x}}{\partial \dot{\phi }}& 0& 0\\ 0& 0& 1& T& 0& 0\\ T\cdot \frac{\partial \ddot{\phi }}{\partial \dot{y}}& T\cdot \frac{\partial \ddot{\phi }}{\partial \dot{x}}& 0& 1+T\cdot \frac{\partial \ddot{\phi }}{\partial \dot{\phi }}& 0& 0\\ T\cdot \cos \phi & T\cdot \sin \phi & T\cdot \left(\dot{x}\cos \phi -\dot{y}\sin \phi \right)& 0& 1& 0\\ T\cdot \left(-\sin \phi \right)& T\cdot \cos \phi & T\cdot \left(-\dot{x}\sin \phi -\dot{y}\cos \phi \right)& 0& 0& 1\end{array}\right], \\ B\left(k\right)=\left[\begin{array}{cccccc}T\cdot \frac{\partial \ddot{y}}{\partial {\delta }_f}& T\cdot \frac{\partial \ddot{x}}{\partial {\delta }_f}& 0& T\cdot \frac{\partial \ddot{\phi }}{\partial {\delta }_f}& 0& 0\end{array}\right].\end{array}$$

among them:

$$\begin{array}{l}\frac{\partial \ddot{y}}{\partial \dot{y}}=\frac{2\sin {\delta }_f}{m}\cdot \frac{\partial F_{lf}}{\partial {\alpha }_f}\cdot \frac{\partial {\alpha }_f}{\partial \dot{y}}+\frac{2\cos {\delta }_f}{m}\cdot \frac{\partial F_{cf}}{\partial {\alpha }_f}\cdot \frac{\partial {\alpha }_f}{\partial \dot{y}}+\frac{2}{m}\cdot \frac{\partial F_{cr}}{\partial {\alpha }_r}\cdot \frac{\partial {\alpha }_r}{\partial \dot{y}}\\ \frac{\partial \ddot{y}}{\partial \dot{x}}=-\dot{\phi }+\frac{2\sin {\delta }_f}{m}\cdot \frac{\partial F_{lf}}{\partial {\alpha }_f}\cdot \frac{\partial {\alpha }_f}{\partial \dot{x}}+\frac{2\cos {\delta }_f}{m}\cdot \frac{\partial F_{cf}}{\partial {\alpha }_f}\cdot \frac{\partial {\alpha }_f}{\partial \dot{x}}+\frac{2}{m}\cdot \frac{\partial F_{cr}}{\partial {\alpha }_r}\cdot \frac{\partial {\alpha }_r}{\partial \dot{x}}\\ \frac{\partial \ddot{y}}{\partial \dot{\phi }}=-\dot{x}+\frac{2\sin {\delta }_f}{m}\cdot \frac{\partial F_{lf}}{\partial {\alpha }_f}\cdot \frac{\partial {\alpha }_f}{\partial \dot{\phi }}+\frac{2\cos {\delta }_f}{m}\cdot \frac{\partial F_{cf}}{\partial {\alpha }_f}\cdot \frac{\partial {\alpha }_f}{\partial \dot{\phi }}+\frac{2}{m}\cdot \frac{\partial F_{cr}}{\partial {\alpha }_r}\cdot \frac{\partial {\alpha }_r}{\partial \dot{\phi }}\\ \frac{\partial \ddot{x}}{\partial \dot{y}}=\dot{\phi }+\frac{2\cos {\delta }_f}{m}\cdot \frac{\partial F_{lf}}{\partial {\alpha }_f}\cdot \frac{\partial {\alpha }_f}{\partial \dot{y}}-\frac{2\sin {\delta }_f}{m}\cdot \frac{\partial F_{cf}}{\partial {\alpha }_f}\cdot \frac{\partial {\alpha }_f}{\partial \dot{y}}+\frac{2}{m}\cdot \frac{\partial F_{lr}}{\partial {\alpha }_r}\cdot \frac{\partial {\alpha }_r}{\partial \dot{y}}\\ \frac{\partial \ddot{x}}{\partial \dot{x}}=\frac{2\cos {\delta }_f}{m}\cdot \frac{\partial F_{lf}}{\partial {\alpha }_f}\cdot \frac{\partial {\alpha }_f}{\partial \dot{x}}-\frac{2\sin {\delta }_f}{m}\cdot \frac{\partial F_{cf}}{\partial {\alpha }_f}\cdot \frac{\partial {\alpha }_f}{\partial \dot{x}}+\frac{2}{m}\cdot \frac{\partial F_{lr}}{\partial {\alpha }_r}\cdot \frac{\partial {\alpha }_r}{\partial \dot{x}}\\ \frac{\partial \ddot{x}}{\partial \dot{\phi }}=\dot{y}+\frac{2\cos {\delta }_f}{m}\cdot \frac{\partial F_{lf}}{\partial {\alpha }_f}\cdot \frac{\partial {\alpha }_f}{\partial \dot{\phi }}-\frac{2\sin {\delta }_f}{m}\cdot \frac{\partial F_{cf}}{\partial {\alpha }_f}\cdot \frac{\partial {\alpha }_f}{\partial \dot{\phi }}+\frac{2}{m}\cdot \frac{\partial F_{lr}}{\partial {\alpha }_r}\cdot \frac{\partial {\alpha }_r}{\partial \dot{\phi }}\\ \frac{\partial \ddot{\phi }}{\partial \dot{y}}=\frac{2l_f\sin {\delta }_f}{I_z}\cdot \frac{\partial F_{lf}}{\partial {\alpha }_f}\cdot \frac{\partial {\alpha }_f}{\partial \dot{y}}+\frac{2l_f\cos {\delta }_f}{I_z}\cdot \frac{\partial F_{cf}}{\partial {\alpha }_f}\cdot \frac{\partial {\alpha }_f}{\partial \dot{y}}-\frac{2l_r}{I_z}\cdot \frac{\partial F_{cr}}{\partial {\alpha }_r}\cdot \frac{\partial {\alpha }_r}{\partial \dot{y}}\\ \frac{\partial \ddot{\phi }}{\partial \dot{x}}=\frac{2l_f\sin {\delta }_f}{I_z}\cdot \frac{\partial F_{lf}}{\partial {\alpha }_f}\cdot \frac{\partial {\alpha }_f}{\partial \dot{x}}+\frac{2l_f\cos {\delta }_f}{I_z}\cdot \frac{\partial F_{cf}}{\partial {\alpha }_f}\cdot \frac{\partial {\alpha }_f}{\partial \dot{x}}-\frac{2l_r}{I_z}\cdot \frac{\partial F_{cr}}{\partial {\alpha }_r}\cdot \frac{\partial {\alpha }_r}{\partial \dot{x}}\\ \frac{\partial \ddot{\phi }}{\partial \dot{\phi }}=\frac{2l_f\sin {\delta }_f}{I_z}\cdot \frac{\partial F_{lf}}{\partial {\alpha }_f}\cdot \frac{\partial {\alpha }_f}{\partial \dot{\phi }}+\frac{2l_f\cos {\delta }_f}{I_z}\cdot \frac{\partial F_{cf}}{\partial {\alpha }_f}\cdot \frac{\partial {\alpha }_f}{\partial \dot{\phi }}-\frac{2l_r}{I_z}\cdot \frac{\partial F_{cr}}{\partial {\alpha }_r}\cdot \frac{\partial {\alpha }_r}{\partial \dot{\phi }}\\ \frac{\partial \ddot{y}}{\partial {\delta }_f}=\frac{2\sin {\delta }_f}{m}\cdot \frac{\partial F_{lf}}{\partial {\alpha }_f}\cdot \frac{\partial {\alpha }_f}{\partial {\delta }_f}+\frac{2F_{lf}\cos {\delta }_f}{m}+\frac{2\cos {\delta }_f}{m}\cdot \frac{\partial F_{cf}}{\partial {\alpha }_f}\cdot \frac{\partial {\alpha }_f}{\partial {\delta }_f}-\frac{2F_{cf}\sin {\delta }_f}{m}\\ \frac{\partial \ddot{y}}{\partial \dot{x}}=\frac{2\cos {\delta }_f}{m}\cdot \frac{\partial F_{lf}}{\partial {\alpha }_f}\cdot \frac{\partial {\alpha }_f}{\partial {\delta }_f}-\frac{2F_{lf}\sin {\delta }_f}{m}-\frac{2\sin {\delta }_f}{m}\cdot \frac{\partial F_{cf}}{\partial {\alpha }_f}\cdot \frac{\partial {\alpha }_f}{\partial {\delta }_f}-\frac{2F_{cf}\cos {\delta }_f}{m}\\ \frac{\partial \ddot{y}}{\partial \dot{\phi }}=\frac{2l_f\sin {\delta }_f}{I_z}\cdot \frac{\partial F_{lf}}{\partial {\alpha }_f}\cdot \frac{\partial {\alpha }_f}{\partial {\delta }_f}+\frac{2l_f{F}_{lf}\cos {\delta }_f}{I_z}+\frac{2l_f\cos {\delta }_f}{I_z}\cdot \frac{\partial F_{cf}}{\partial {\alpha }_f}\cdot \frac{\partial {\alpha }_f}{\partial {\delta }_f}-\frac{2l_f{F}_{cf}\sin {\delta }_f}{I_z}\end{array}$$

In the above formula, the partial derivatives of the side slip angle of the front and rear wheels ${\alpha }_f,{\alpha }_r$ with respect to the longitudinal speed $\dot{x}$, lateral speed $\dot{y}$, yaw rate $\dot{\phi }$ and front wheel rotation angle ${\delta }_f$ can be obtained by calculating the partial derivatives of Equations (6)–(8).

Since the control variable is the front steering angle, in order to ensure the normal running of the vehicle and prevent sudden changes in the front steering angle, it is necessary to limit the front steering angle increment, so the control variable in Equation (26) is converted into a control increment [33]:

$$\overline{x}\left(k+1\right)=\tilde{\mathit{A}}\overline{\mathit{x}}\left(k\right)+\tilde{B}\Delta \tilde{u}\left(k\right)$$

(27)

In the formula: $\overline{x}\left(k\right)=\left[\begin{array}{c}\tilde{x}\left(k\right)\\ \tilde{u}\left(k-1\right)\end{array}\right],\tilde{A}=\left[\begin{array}{cc}{A}_{6\times 6}& B_{6\times 1}\\ 0_{1\times 6}& I_{1\times 1}\end{array}\right],\tilde{B}=\left[\begin{array}{c}{B}_{6\times 1}\\ I_{1\times 1}\end{array}\right],\Delta \tilde{u}\left(k\right)=\tilde{u}\left(k\right)-\tilde{u}\left(k-1\right)$.

The target path tracking needs to select the output of the state equation. The selected output is the yaw angle and the lateral offset. The former can be used to track and detect the forward direction of the vehicle, and the latter can be used to track and detect the forward position of the vehicle, and the output equation of the prediction model [33] is as follows [33]:

$$y\left(k\right)=\tilde{C}\overline{x}\left(k\right)$$

(28)

In the formula: $\tilde{C}=\left[\begin{array}{ccccccc}0& 0& 1& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 1& 0\end{array}\right]$.

Combining Equations (27) and (28) as a complete prediction model, Equation (29):

$$\{\begin{array}{c}\overline{x}\left(k+1\right)=\tilde{A}\overline{x}\left(k\right)+\tilde{B}\Delta \tilde{u}\left(k\right)\\ y\left(k\right)=\tilde{C}\overline{x}\left(k\right)\end{array}$$

(29)

3.2. Output in the Forecast Time Domain

It is assumed that the prediction time domain is $N_p$, the control time domain is $N_c$, and the prediction time domain is greater than the control time domain. Based on Equation (29), it is predicted that $N_c$ control increments in the control time domain will act on the system. When time $k$ is in the control time domain between the time domain and the prediction time domain, the input control amount of the system remains unchanged, and the control increment is 0. After derivation, the system’s predicted output expression can be obtained [33]:

$$\tilde{y}\left(k\right)=\overline{A}\left(k\right)\overline{x}\left(k|k\right)+\overline{B}\left(k\right)\Delta \tilde{U}\left(k\right)$$

(30)

In the formula:

$$\begin{array}{l}\tilde{y}\left(k\right)=\left[\begin{array}{c}y\left(k+1|k\right)\\ ⋮\\ y\left(k+N_p|k\right)\end{array}\right]\Delta \tilde{U}\left(k\right)=\left[\begin{array}{c}\Delta \tilde{u}\left(k|k\right)\\ ⋮\\ \Delta \tilde{u}\left(k+N_c-1|k\right)\end{array}\right]\\ \overline{A}\left(k\right)=\left[\begin{array}{c}\tilde{C}\tilde{A}\\ ⋮\\ \tilde{C}{\tilde{A}}^{N_p}\end{array}\right]\overline{B}\left(k\right)=\left[\begin{array}{ccc}\tilde{C}\tilde{B}& 0& 0\\ ⋮& ⋮& 0\\ \tilde{C}{\tilde{A}}^{N_c-1}\tilde{B}& \cdots & \tilde{C}\tilde{B}\\ ⋮& & ⋮\\ \tilde{C}{\tilde{A}}^{N_p-1}\tilde{B}& \cdots & \tilde{C}{\tilde{A}}^{N_p-N_c}\tilde{B}\end{array}\right]\end{array}$$

3.3. Objective Function Design

The objective function is the performance standard of the system. To ensure that the unmanned vehicle can quickly and smoothly follow the target path and reflect its tracking ability, it is also necessary to optimize the deviation and control of the system state variable. It is necessary to ensure that the control increment in each sampling period is constrained to avoid sudden changes. Therefore, referring to the soft constraint method used in the literature [35], the selected objective function form is as follows:

$$\begin{array}{ll}J=& {\displaystyle \sum _{i=1}^{N_p}{\Vert \tilde{y}(k+i|k)-y_{ref}(k+i|k)\Vert }_Q^2+{\displaystyle \sum _{i=0}^{N_c-1}{\Vert \Delta u(k+i|k)\Vert }_R^2}}+\rho {\epsilon }^2\\ & ={\displaystyle \sum _{i=1}^{N_p}q{[\tilde{y}(k+i|k)-y_{ref}(k+i|k)]}^2}\hspace{0.17em}\\ & +{\displaystyle \sum _{i=0}^{N_c-1}r{\left[u\left(k+i|k\right)-u\left(k+i-1|k\right)\right]}^2}+\rho {\epsilon }^2\end{array}$$

(31)

In the formula, $\rho$ is the weight coefficient, $\epsilon$ is the relaxation factor, and $Q,R$ are the output weighting matrix and the control weighting matrix, respectively. The first item represents the accumulation of deviations between the output variable and the reference value of the output variable, which reflects the ability to follow the reference trajectory, and the second item is the accumulation of deviations of the control variable, which reflects the requirement for the steady change of the control variable. The third item ensures that there is a feasible optimal solution when changing in real time. $\tilde{y}(k+i|k)$ is the predicted value of the output at time $k$ in the future at time $k+i$, $y_{ref}(k+i|k)$ is the reference value of the output at time $k$ in the future at time $k+i$, and $\Delta u\left(k+i|k\right)$ is the future control variable sequence.

The output weighting matrix: $Q={\left[\begin{array}{cccc}q& & & \\ & q& & \\ & & \ddots & \\ & & & q\end{array}\right]}_{\left(2\times N_p\right)\times \left(2\times N_p\right)}, q={\left[\begin{array}{cc}{q}_1& \\ & q_2\end{array}\right]}_{2\times 2}$.

The control weight matrix: $R={\left[\begin{array}{cccc}r& & & \\ & r& & \\ & & \ddots & \\ & & & r\end{array}\right]}_{N_c\times N_c}$.

In order to simplify the calculation, the control increment matrix A is taken as a single design variable, so the objective function [36] after matrix transformation is as follows:

$$J=\frac{1}{2}\xi {\left(k\right)}^{T}H\left(k\right)\xi \left(k\right)+f\left(k\right)\xi \left(k\right)+P\left(k\right)$$

(32)

In the formula:

$$\begin{array}{l}H\left(k\right)=\left[\begin{array}{cc}2\left(\Theta {\left(k\right)}^{T}Q\Theta \left(k\right)+R\right)& 0\\ 0& 2\rho \end{array}\right], \xi \left(k\right)=\left[\begin{array}{c}\Delta U\left(k\right)\\ \epsilon \end{array}\right], f\left(k\right)=\left[\begin{array}{cc}2E{\left(k\right)}^{T}Q\Theta \left(k\right)& 0\end{array}\right],\\ P\left(k\right)=E{\left(k\right)}^{T}QE\left(k\right), E\left(k\right)=\overline{A}\left(k\right)\overline{x}\left(k|k\right)-y_{ref}(k|k).\end{array}$$

3.4. Design of Constraints

The constraint conditions of model predictive control generally include control variable constraint, control increment constraint and output variable constraint. In real vehicle driving scenarios, the front steering angle has a certain range limit, so the actual situation should be considered when designing the controller, and the control variable front steering angle should be restricted. At the same time, in order to consider the stability of the vehicle, the constraint on the control increment is also indispensable, and the change range of the control increment cannot be too large. In addition, in order to avoid excessive deviation from the target path, it is necessary to refer to the target path to constrain the output yaw angle and lateral offset [36].

Impose constraints on the control amount:

$$U_{\min }\left(k\right)\le U\left(k\right)\le U_{\max }\left(k\right)$$

(33)

In the formula, $U_{\min }\left(k\right)$ is the minimum value of the control variable, and $U_{\max }\left(k\right)$ is the maximum value of the control variable.

Impose constraints on the control increment:

$$\Delta U_{\min }\left(k\right)\le \Delta U\left(k\right)\le \Delta U_{\max }\left(k\right)$$

(34)

In the formula, $\Delta U_{\min }\left(k\right)$ is the minimum value of the control increment, and $\Delta U_{\max }\left(k\right)$ is the maximum value of the control increment.

Impose constraints on the output:

$${\tilde{y}}_{\min }\left(k\right)\le \tilde{y}\left(k\right)\le {\tilde{y}}_{\max }\left(k\right)$$

(35)

In the formula, ${\tilde{y}}_{\min }\left(k\right)$ is the minimum value of the output and ${\tilde{y}}_{\max }\left(k\right)$ is the maximum value of the output.

4. Simulation

In order to verify the performance of the proposed MPC controller based on the non-linear magic formula tires, a simulation test was carried out using the control function of Matlab/Simulink and Carsim’s vehicle dynamics model and vehicle operating environment platform.

4.1. Parameter Selection

Select the E-Class Sedan model in Carsim for simulation, and the main parameters are shown in

Table 2. Vehicle parameters.

Parameter Name	Parameter Value	Parameter Name	Parameter Value
Vehicle quality/kg	$m=1723$	Distance from the center of mass of the vehicle to the front axle/m	$l_f=1.232$
Moment of inertia around the z axis/(kg·m2)	$I_z=4175$	Distance from the center of mass of the vehicle to the rear axle/m	$l_r=1.468$
Vehicle wheelbase/m	$L=2.7$

.

This paper mainly focuses on the lateral force of the tire. According to Equations (16)–(18) and the data in Table 1 and Table 2, the relationship between the lateral force of the tire and the tire slip angle is shown in Figure 2.

In line with [33], in the lateral tracking ability experiment of unmanned vehicles, the front wheel of the vehicle is continuously turned from the left extreme position to the right extreme position. Record the corresponding time and position, and then you can measure the limit value and time of the front steering angle in clockwise and counterclockwise directions. Therefore, the constraint on the front steering angle and the front steering angle increment and the constraint on the output based on the double-shifting reference target path are shown in Equation (36):

$$\begin{array}{l}-{10}^{°}\le {\delta }_f\le {10}^{°}\\ -{0.85}^{°}\le \Delta {\delta }_f\le {0.85}^{°}\\ -0.3\le \phi \le 0.21\\ -3\le Y\le 5\end{array}$$

(36)

4.2. Constant Speed Simulation

4.2.1. Reference Path Selection

In the test of vehicle handling and stability, the double lane change is used widely [37,38], so the selected reference path equation is as follows [21]

$$Y\left(X\right)=\frac{d_{y1}}{2}\left(1+\tan \left(z_1\right)\right)-\frac{d_{y2}}{2}\left(1+\tan \left(z_2\right)\right)$$

(37)

In the formula:

$$z_1=\frac{2.4}{25}\left(X-27.19\right)-1.2,z_2=\frac{2.4}{21.95}\left(X-56.46\right)-1.2,d_{y1}=4.05,d_{y2}=5.7.$$

The reference path is shown in Figure 3:

4.2.2. Simulation Results of Constant Speed Conditions

The controller, designed based on the magic formula tire 175/70 R13 (asymmetric), takes into account the situation in which the tire part enters the nonlinear region. It is referred to as controller A below. In order to reflect the performance of controller A, it is compared with the conventional model predictive controller B. Controller B does not consider the tire entering the nonlinear region, and is only designed based on the simplification of the tire’s longitudinal stiffness and cornering stiffness. The stiffness parameters of the magic formula tire are: $C_{cf}=\mathrm{48,400},C_{cr}=\mathrm{44,800},C_{lf}=\mathrm{90,800},C_{lr}=\mathrm{76,200}$.

(1)

First select a good road with a road adhesion coefficient of 0.8; the vehicle speed is 36 km/h, and the controller parameters are:

$$T=0.05 {\mathrm{s},\mathrm{N}}_p=10,{\mathrm{N}}_c=3, Q=\left[\begin{array}{cc}2000& 0\\ 0& 1000\end{array}\right], R=5\times {10}^5,\rho =1000, \epsilon =1000.$$

The simulation results are shown in Figure 4a–f, Figure 4a is the comparison chart of the tracking reference path, Figure 4b is the comparison chart of the front steering angle, Figure 4c is the comparison chart of the yaw rate, Figure 4d is the comparison chart of the side slip angle of the center of mass, Figure 4e is the comparison chart of lateral acceleration, and Figure 4f is the comparison chart of the tire lateral force.

First of all, it can be seen from Figure 4e that the lateral accelerations of controller A and controller B are both within 0.4 g when the vehicle speed is 36 km/h and the adhesion coefficient is 0.8. Furthermore, it can be seen from Figure 2 that the saturation value of the tire force is about 3.6 kN; combined with Figure 4f, it can be seen that the tire lateral force is in the linear region. At this time, it can be seen from Figure 4a that in the interval of 0–70 m, the tracking effects of controller A and controller B are similar, but in the interval of 70–120 m, it is obvious that the tracking effect of controller A is better than that of controller B. Although the tire force is in the linear region at this time, controller A can more accurately describe the tire forces, while controller B describes the simplified linear tire forces, so the accuracy of controller A is better than that of controller B. From Figure 4b–d, it can be seen that the front steering angle, yaw rate and side slip angle of controller A and controller B are relatively close, and the magnitude and shape of the change are similar. Although controller A generates two small fluctuations in the range of 6–10 s, controller A’s is smaller than controller B’s in terms of peak value. In summary, it can be seen that when the lateral acceleration is less than 0.4 g and the tire lateral force is in the linear region, the performance of the two controllers is close, and they can track the reference path well and smoothly, and the overall performance of controller A is slightly better.

(2)

Increase the vehicle speed from 36 km/h to 72 km/h, the road adhesion coefficient is still 0.8, and the parameters of the controller are:

$$T=0.05 {\mathrm{s},\mathrm{N}}_p=16,{\mathrm{N}}_c=3, Q=\left[\begin{array}{cc}2000& 0\\ 0& 1000\end{array}\right], R=5\times {10}^5,\rho =1000, \epsilon =1000.$$

The simulation results are shown in Figure 5a–f, Figure 5a is the comparison chart of the tracking reference path, Figure 5b is the comparison chart of the front steering angle, Figure 5c is the comparison chart of the yaw rate, Figure 5d is a comparison chart of the side slip angle of the center of mass, Figure 5e is the comparison chart of lateral acceleration, and Figure 5f is the comparison chart of tire lateral force.

It can be seen from Figure 5a that when the vehicle speed increases to 72 km/h, the tracking effects of controller A and controller B are similar in the 0–50 m interval. In the interval of 70–85 m, the tracking effect of controller A is not as good as that of controller B, but at 60 m, the peak of controller A is closer to the reference path, and in the interval of 85−120 m, controller A can control the vehicle to track the reference path earlier. On the whole, the tracking accuracy of controller A is better than that of controller B. From Figure 5b–d, the change amplitude and shape of the front steering angle, yaw rate and side slip angle of controller A and controller B are still similar. Combining the three figures and Figure 5a, it can be seen that the tracking accuracy of the two controllers is similar within 0–2 s, but the front steering angle, yaw rate and side slip angle of controller A are all smaller than that of controller B. In the 4–6 s interval, controller A is inferior to controller B. This is because during this period, controller A sacrifices some stability to increase the tracking accuracy and can control the vehicle closer to the reference path more quickly and effectively. It can be seen from Figure 5e that the lateral acceleration of the vehicle under the action of controller B is greater than 0.4 g for some time. In combination with Figure 5f, it can be seen that the tire lateral force has been saturated within 2~4 s, which indicates that the vehicle model has entered the nonlinear region during this period. The linear tire model of controller B cannot accurately describe the tire force, which affects the tracking accuracy and stability of the vehicle. The lateral acceleration of the vehicle in controller A is less than 0.4 g, and the tire force is still in the linear region. The tracking accuracy and vehicle stability are better than that of controller B. In summary, controller B tends to be unstable under this working condition, and controller A can still control the vehicle to track the reference path more smoothly and effectively.

(3)

The vehicle speed is maintained at 72 km/h, the road adhesion coefficient is 0.3, and the parameters of the controller are:

$$T=0.05 {\mathrm{s},\mathrm{N}}_p=24,{\mathrm{N}}_c=4, Q=\left[\begin{array}{cc}2000& 0\\ 0& 1000\end{array}\right], R=5\times {10}^5,\rho =1000, \epsilon =1000.$$

The simulation results are shown in Figure 6a–f, Figure 6a is the tracking reference path diagram of controller A and controller B. Figure 6b is a diagram of the front steering angle of controller A. Figure 6c is the yaw rate diagram of controller A. Figure 6d is a diagram of the side slip angle of the centroid of controller A. Figure 6e is the lateral acceleration diagram of controller A. Figure 6f is the tire lateral force diagram of controller A.

When the vehicle speed remains at 72 km/h and the road adhesion coefficient decreases from 0.8 to 0.3, it is clear that controller B in Figure 6a cannot keep up with the reference path, so only the performance effect of controller A needs to be analyzed.

From Figure 6e, it can be seen that, due to the reduction in the road adhesion coefficient, the peak value of the lateral acceleration of the vehicle exceeds 0.4 g. According to Figure 6f, within 2–6 s, the tire lateral force approaches and partially exceeds the saturation value, and the vehicle model enters the nonlinear region. In Figure 6b, the maximum front wheel rotation angle has reached 12 degrees. The maximum yaw rate in Figure 6c and the maximum lateral deflection angle of the center of mass in Figure 6d undergo a large increase as compared to Figure 5c,d. Driving at high speed on low-attachment roads, the stability of the car decreases, but from Figure 6a, it can be seen that the vehicle under the action of controller A can still keep up with the reference path more effectively as a whole.

4.3. Variable Speed Simulation

4.3.1. Reference Path Selection

The vehicle speed is changed from a constant speed to a variable speed, and a snake path is selected for the test. The reference path equation is as follows:

$$y=\{\begin{array}{l} 0\hspace{0.17em}\hspace{0.17em} ,0\le x\le 20\\ 3.5\ast \sin (\mathrm{pi}/50\ast (x-20),20<x\le 220\\ 0 ,220<x\le 300\end{array}$$

(38)

The reference path diagram is as follows:

4.3.2. Simulation Results of Variable Speed Conditions

In order to investigate the robustness of the controller to the change of vehicle speed, the vehicle speed is set to increase from 36 km/h to 72 km/h within 0 to 10 s, and then maintain a constant speed of 72 km/h. In order to prevent the vehicle from losing control at high speeds, the road adhesion coefficient is set to 0.8, the reference path is the snake path in Figure 7, and the speed change is shown in Figure 8.

The parameters of the controller are as follows:

$$T=0.05 {\mathrm{s},\mathrm{N}}_p=17,{\mathrm{N}}_c=2, Q=\left[\begin{array}{cc}\mathrm{10,000}& 0\\ 0& 2000\end{array}\right], R=5\times {10}^5,\rho =1000, \epsilon =1000.$$

The simulation results are shown in Figure 9a–e: Figure 9a is the comparison chart of the tracking reference path, Figure 9b is the comparison chart of the front steering angle, Figure 9c is the comparison chart of the yaw rate, Figure 9d is a comparison chart of the side slip angle of the center of mass, Figure 9e is the comparison chart of the lateral acceleration, and Figure 9f is the comparison chart of the tire lateral force.

It can be seen from Figure 9a that the tracking accuracy of the vehicle is better in the low speed stage of 0–80 m. As the vehicle speed increases, the tracking accuracy of controller B gradually decreases; while controller A can still maintain good tracking accuracy, several peaks of controller A are closer to the reference path, and the overall tracking effect of controller A is significantly better than that of controller B. Figure 9b–d show that at 12 s, the peak value of controller B is higher than that of controller A, indicating that controller B is inferior to controller A in controlling the vehicle stability at high speeds. At 14 s, the peak value of controller A is larger, because the control adaptability of controller A is stronger. With reference to Figure 9a, in the 200–250 m interval, the vehicle speed has reached 72 km/h at this time. Under the action of controller A, the vehicle can track the reference path more quickly and accurately. It can be seen from Figure 9e that the lateral accelerations of both controllers exceed 0.4 g at 10 s. Furthermore, it can be seen from Figure 9f that at about 10 s, part of the tire lateral force has reached the saturation value, and after 14 s, the tire lateral force of controller B has chattering. The linear calculation of the tire force of controller B is quite different from the real tire force. The accumulated error will affect the later tracking control effect. In contrast to controller A, the tire force calculation is closer to the real tire force, so it can still maintain a better tracking effect. In summary, in the variable speed condition, controller A can still maintain better tracking accuracy and vehicle stability than controller B.

5. Hardware in the Loop Experiment

The simulation analysis in the third section preliminarily verifies the performance effect of controller A under constant speed and variable speed conditions. In order to verify the actual effect of the controller, a hardware-in-the-loop test platform based on Labview-RT is established in this section. In the hardware-in-the-loop test, although commercial simulation software is used to replace the real vehicle, the steering actuator and control system are complete physical objects. Therefore, hardware-in-the-loop testing is used by many automobile manufacturers to verify the logic of vehicle path tracking control and vehicle stability control. The test plan is shown in Figure 10 [39,40,41]. The simulation software Carsim is embedded in the automation platform PXI real-time system, and the vehicle status and road information are transmitted to the dSPACE controller through the PXI CAN card. The designed control scheme is stored in dSPACE, and the active steering angle is calculated in real time based on the vehicle status and road information transmitted by PXI. The steering signal is transmitted to the underlying controller, and the underlying controller drives the steering motor and the resistance loading system to execute the angle signal. The steering wheel is controlled to complete the angle change. The angle sensor attached to the steering wheel mechanism collects the actual steering steering angle, and then transmits the angle signal back to Carsim through the PXI CAN card. The computer collects the motor current graph and the vehicle dynamic response graph in Carsim. Through the above steps, the hardware-in-the-loop test is completed.

The hardware-in-the-loop test condition is selected as the constant speed condition ($\mu =0.8, u=72 \mathrm{km}/\mathrm{h}$) under a high adhesion coefficient, and the tracking path is still double lance change. The basic parameter settings of the controller are the same as the simulation condition. The test results are shown in Figure 11a–g; the results verify the effectiveness of the designed controller.

Comparing Figure 5a and Figure 11a, the path tracking effects of the two controllers in simulation and experiment are similar. Controller A is closer to the reference path at 60 m, while controller B is slightly better at 60–80 m, but controller A can keep up with the reference path more quickly and effectively after 80 m. The amplitudes and shapes of the experimental graphs 11b–d and simulation graphs 5b–d are also similar. The peak value of controller A is smaller than that of controller B. It can be seen from Figure 11e that the lateral acceleration of controller A is still within 0.4 g, while the upper and lower peaks of controller B exceed 0.4 g. Furthermore, it can be seen from Figure 11f that the tire lateral force reaches its saturation value within 2–3 s. After the tire force enters the nonlinear region, the tire force calculated by controller B has errors in terms of the real tire force, and the accumulation of errors affects the later tracking accuracy and vehicle stability.

Figure 11g is the collected motor current graph. It can be clearly seen that the current of controller B has reached ±30 A, while the current of controller A is only within ±20 A. After 8 s, the motor current of controller B is still fluctuating, while the motor current of controller A tends to be stable. The overall noise fluctuation of controller A is small, the current is more stable, and the response speed is faster. In summary, considering the path tracking accuracy and vehicle stability, the effect of controller A is better.

6. Conclusions

In this paper, a model predictive algorithm controller, based on a nonlinear tire model, is designed. It is used to reduce the calculation error of the tire force of the vehicle in the case of poor road conditions, high speed, etc., to avoid the deterioration of the path tracking effect due to excessive error. In the controller, the tire forces are characterized with nonlinear composite functions of the magic formula instead of a simple linear relation model. In the design of the controller, Taylor expansion is used for linearization, the first-order difference quotient method is used for discretization, and the partial derivative of the composite function is used for matrix transformation. Constant velocity and variable velocity conditions are selected to compare the designed controller with the conventional linear tire force controller in Carsim/Simulink. The results show that when the tire force does not fall into the nonlinear region, the stability of the two controllers is similar. The tracking accuracy of the controller designed in this paper is slightly better. However, after the tire force enters the nonlinear region, the effect of the controller, when using linear tire force, becomes worse, the tracking accuracy is far worse than the controller when using the magic formula tire design, and the vehicle stability is also degraded. In addition, a semi-physical test platform based on the LabVIEW-RT system was established, and hardware-in-the-loop test experiments were carried out. The test results verified the effectiveness of the designed controller. In order to further improve the tracking performance and stability of unmanned vehicles, the MPC model, based on nonlinear tires, that is proposed in this study will continue to be expanded in the future, so that autonomous vehicles can effectively operate in more complex and realistic traffic scenarios.