Research on Longitudinal Control of Electric Vehicle Platoons Based on Robust UKF-MPC

1. Introduction

With the popularization of electric vehicles (EVs) and the rapid development of telematics, the intelligent control of EV platoons has become an important application in intelligent transportation systems. This technology can effectively improve transportation safety, enhance efficiency, and reduce energy consumption [1,2].

In EV platoons, vehicles exchange real-time information through wireless communication technologies such as vehicle-to-vehicle (V2V) communication, enabling collaborative control and maintaining desired spacing for safe travel. However, the complexity of traffic environments and the nonlinear nature of vehicle dynamics systems present challenges, including communication delays, data packet loss, and external interference, which can adversely affect vehicle behavior and weaken platoon stability [3,4].

To address these challenges and enhance the stability and control precision of vehicle platoons, numerous researchers have conducted studies. Yang et al. [5] addressed the robustness of cooperative adaptive cruise control (CACC) systems under cyber-attacks by proposing a method to utilize multiple V2V communication networks for data fusion and redundant transmission, and designing an H∞ controller to minimize the impact of sensor and channel noise on vehicle platoon performance. Li et al. [6] proposed a distributed nonlinear vehicle platoon longitudinal controller based on a third-order model for heterogeneous vehicle platoons in communication delay environments. Samii et al. [7] developed a linear predictive feedback Cooperative Adaptive Cruise Control (CACC) controller to address communication and actuator delays in heterogeneous platoons. Halder et al. [8] introduced a discrete-time distributed state feedback control strategy for homogeneous vehicle platoons with an undirected network topology, capable of resisting external disturbances and random continuous network packet losses. Lu et al. [9] created a distributed model predictive control (DMPC) strategy to ensure platoon stability for nonlinear vehicle platoons with a unidirectional communication topology, considering sensor jitter, control delays, and errors in real vehicle conditions. Meng et al. [10] designed a robust MPC controller and a slip ratio controller to address control errors due to parameter uncertainties in platoon modeling and ensure the stability of following vehicles. Tian et al. [11] employed Model Predictive Control (MPC) and Long Short-Term Memory (LSTM) prediction methods to study communication delay compensation in vehicle coordination control under CACC systems. Wang et al. [12] proposed an improved MPC algorithm based on the Kalman Filter (KF) to tackle issues such as environmental disturbances, sensor noise, and poor stability in following time-varying speeds.

In summary, while research on the longitudinal control of EV platoons has progressed in addressing external interference, communication delay, and data packet loss, the randomness of these factors still impacts vehicle state accuracy and platoon stability. To solve these issues, this paper constructs a longitudinal kinematic model of a vehicle platoon and establishes discrete state-space equations. By introducing Huber-M estimation to enhance the UKF algorithm, vehicle state information within the platoon is accurately estimated. Based on this estimated state information, feedback correction compensation is integrated into the MPC algorithm to design a robust UKF-MPC longitudinal controller. Finally, the effectiveness of this controller is verified through CarSim/Simulink co-simulation.

2. Problem Description

Consider a traffic scenario in which N electric vehicles form a platoon travelling in a straight line, as shown in Figure 1, where the leading vehicle is numbered $i=0$ and the following vehicles behind are numbered $i=\mathrm{1,2},...,N$. The following vehicle in the platoon adopts the Predecessor Following (PF) communication topology form through the V2V communication technology to obtain the position, velocity, acceleration and other state information of the vehicle in front of it. The distance between the following vehicle and the front vehicle is obtained through on-board sensors (e.g., radar, camera, etc.).

3. Longitudinal Kinematic Model of the Vehicle Platoon

Based on the consideration of various factors such as stability, safety and the complexity of the use of working conditions, for the desired distance $d_{i-1,i\_des}$ between neighboring vehicles, the desired distance is more easily adjusted with the Constant Time Headway (CTH) as the vehicle platoon spacing strategy.

$$d_{i-1,i\_des}=d_{min}+{\tau }_h{v}_i$$

(1)

where $d_{min}$ represents the minimum safe distance, ${\tau }_h$ represents the headway, $v_i$ represents the velocity of vehicle i.

According to the position information obtained from V2V communication, the workshop distance $d_{i-1,i\_V2V}$ of adjacent vehicles is calculated. Considering the communication delay problem and for the driving safety, the calculated workshop distance is combined with the workshop distance $d_{i-1,i\_sensor}$ obtained from the on-board sensors to obtain the workshop distance $d_{i-1,i}$ of the neighboring vehicles, and the calculation formula is as follows:

$$d_{i-1,i\_V2V}=p_{i-1}-p_i-l_i$$

(2)

$$d_{i-1,i}=(1-\zeta )d_{i-1,i\_V2V}+\zeta d_{i-1,i\_sensor}$$

(3)

where $p_{i-1}$ and $p_i$ represent the position information of neighboring vehicles, $l_i$ represents the length of vehicle i, $\zeta =1-{\displaystyle \frac{v_i}{v_{i\_max}}}$ represents the weight parameter, and $v_i$ represents the maximum velocity of vehicle i.

According to the longitudinal kinematics between vehicle i-1 and vehicle i, the inter-vehicle spacing error $\delta d_{i-1,i}$, and the relative velocity $v_{i-1,i}$, can be calculated as follows:

$$\delta d_{i-1,i}=d_{i-1,i}-d_{i-1,i\_des}$$

(4)

$$v_{i-1,i}=v_{i-1}-v_i$$

(5)

Considering that there is a time delay problem in the vehicle actuation system, the relationship between the actual acceleration and the desired acceleration of vehicle i is described by adding a first-order inertial link, which is expressed by the following relation:

$$a_i={\displaystyle \frac{K}{\tau s+1}}{a}_{i\_des}$$

(6)

where K represents the first-order system gain, which takes the value of 1, and $\tau$ represents the inertial link time constant.

According to the above model and longitudinal kinematics a discrete time domain relational equation can be obtained:

$$\left\{\begin{array}{ll}\hfill d_{i-1,i}(k+1)=& d_{i-1,i}(k)+v_{i-1,i}(k)T_s-{\displaystyle \frac{1}{2}}{a}_i(k)T_s^2+{\displaystyle \frac{1}{2}}{a}_{i-1}(k)T_s^2\hfill \\ \hfill v_{i-1,i}(k+1)=& v_{i-1,i}(k)-a_i(k)T_s+a_{i-1}(k)T_s\hfill \\ \hfill v_i(k+1)=& v_i(k)+a_i(k)T_s\hfill \\ \hfill a_i(k+1)=& (1-{\displaystyle \frac{K{T}_s}{\tau }})a_i(k)+{\displaystyle \frac{K{T}_s}{\tau }}{a}_{i\_des}(k)\hfill \end{array}\right.$$

(7)

where k represents the current moment of the system, k+1 represents the next moment of the system, and $T_s$ represents the sampling time of the system.

Define $x(k)=[d_{i-1,i}(k),v_{i-1,i}(k),v_i(k),a_i(k)]$ as the state variable, $u(k)=a_{i\_des}(k)$ as the control variable, and $w(k)=a_{i-1}(k)$ as the disturbance variable. Rewrite Equation (7) into the form of discrete state space equations:

$$x(k+1)=A_{o}x(k)+B_{o}u(k)+C_{o}w(k)$$

(8)

The coefficient matrix is $A_o=\left[\begin{array}{cccc}1& T_s& 0& -{\displaystyle \frac{T_s^2}{2}}\\ 0& 1& 0& -T_s\\ 0& 0& 1& T_s\\ 0& 0& 0& 1-{\displaystyle \frac{K{T}_s}{\tau }}\end{array}\right]$, $B_o=\left[\begin{array}{c}0\\ 0\\ 0\\ {\displaystyle \frac{K{T}_s}{\tau }}\end{array}\right]$, $C_o=\left[\begin{array}{c}{\displaystyle \frac{T_s}{2}}\\ T_s\\ 0\\ 0\end{array}\right]$.

4. Platoon Longitudinal Controller Design

Design the longitudinal controller of vehicle i as an upper and lower layer. The upper controller obtains the state information of the vehicle i-1 and the state information of the vehicle i through the V2V communication to decide the desired acceleration, The lower controller controls the acceleration and deceleration of the vehicle i according to the desired acceleration. Considering that the V2V communication has the problems of random communication delay and data loss, the robust UKF algorithm is used to estimate the state information of the vehicle under the influence of the communication problem. The specific control framework is shown in Figure 2.

4.1. Robust UKF State Prediction Algorithm

Vehicle nonlinear systems:

$$\left\{\begin{array}{l}{x}_k=F(x_{k-1})+w_k\\ z_k=H(x_k)+v_k\end{array}\right.$$

(9)

where $x_k$ represents the state vector, $z_k$ represents the observation vector, $F(x_{k-1})$ represents the state transfer equation, $H(x_k)$ represents the observation equation, $w_k$ represents the system process noise, and $v_k$ represents the observation noise.

The specific steps of robust UKF are as follows:

1.

Initialize the mean ${\mu }_0$ and covariance $P_0$ of the initial state of the system $x_0$.

$$\left\{\begin{array}{l}{\mu }_0=E(x_0)\\ P_0=E\left[(x_0-\overline{x_0})(x_0-\overline{x_0}{)}^T\right]\end{array}\right.$$

(10)

2.

Calculate the sigma points and construct 2n+1 sigma points along with the corresponding weights.

$${\chi }_{k-1}^j=\left\{\begin{array}{l}{\widehat{x}}_{k-1},j=0\\ {\widehat{x}}_{k-1}+\sqrt{\left(n+\lambda \right)}(\sqrt{P_{k-1}}{)}_j,j=1,\cdots ,n\\ {\widehat{x}}_{k-1}-\sqrt{\left(n+\lambda \right)}(\sqrt{P_{k-1}}{)}_j,j=n+1,\cdots 2n\end{array}\right.$$

(11)

where $P_{k-1}$ represents the state estimation covariance array at k-1 time, $(\sqrt{P_{k-1}}{)}_j$ represents the j-th column of the Cholesky decomposition of the $P_{k-1}$ matrix, $\lambda ={\alpha }^2(n+\kappa )-n$. $\lambda$ represents the scale adjustment factor, n represents the number of dimensions of the system state, $\alpha$ and $\kappa$ represent the distributional scaling parameters for determining the sampling points, and $\alpha =1$, $\kappa =3-n$.

The weighting coefficients for the mean and covariance are, respectively:

$$\left\{\begin{array}{l}{w}_m^0={\displaystyle \frac{\lambda }{(n+\lambda )}}\\ w_c^0={\displaystyle \frac{\lambda }{(n+\lambda )}}+(1+\beta -{\alpha }^2)\\ w_m^j=w_c^j={\displaystyle \frac{1}{2(n+\lambda )}},j=1,\cdots 2n\end{array}\right.$$

(12)

where $w_m$ represents the mean weight factor, $w_c$ represents the covariance weight factor, and $\beta$ represents a non-negative weighting factor to limit the error caused by the higher order terms, taking the value of 2.

3.

Forecast update.

Prediction system sigma point set:

$${\chi }_{k|k-1}^j=F({\chi }_{k-1}^j,u_{k-1})$$

(13)

Estimates of the state of the predictive system:

$${\widehat{\mu }}_{x,k}={\sum }_{i=0}^{2n}{\omega }_m^j{\chi }_{k|k-1}^j$$

(14)

Predicting system state covariance:

$${\stackrel{-}{P}}_{x,k}={\sum }_{i=0}^{2n}{\omega }_c^j({\chi }_{k|k-1}^j-{\widehat{\mu }}_{x,k})\times ({\chi }_{k|k-1}^j-{\widehat{\mu }}_{x,k}{)}^T+Q_k$$

(15)

where $Q_k$ represents the covariance matrix of the system process noise.

4.

Measurement update.

Nonlinear transformation of each Sigma point set using measurement equations.

$${\gamma }_{k|k-1}^j=H({\chi }_{k|k-1}^j)$$

(16)

Estimated value of the measurement:

$${\widehat{\mu }}_{z,k}=\sum _{i=0}^{2n}{\omega }_m^j{\gamma }_{k|k-1}^j$$

(17)

Measured covariance:

$${\stackrel{-}{P}}_{z,k}=\sum _{i=0}^{2n}{\omega }_c^j({\gamma }_{k|k-1}^j-{\widehat{\mu }}_{z,k})\times ({\gamma }_{k|k-1}^j-{\widehat{\mu }}_{z,k}{)}^T+R_k$$

(18)

where $R_k$ represents the covariance matrix of the system measurement noise.

5.

Calculate the gain.

Calculation of cross-covariances:

$${\stackrel{-}{P}}_{xz,k}=\sum _{i=0}^{2n}{\omega }_c^j({\chi }_{k|k-1}^j-{\widehat{\mu }}_{x,k})\times ({\gamma }_{k|k-1}^j-{\widehat{\mu }}_{z,k}{)}^T$$

(19)

Calculate the Kalman gain:

$$K_k={\stackrel{-}{P}}_{z,k}{\stackrel{-}{P}}_{xz,k}^{-1}$$

(20)

6.

Huber-M estimation.

Calculate the residuals:

$$e_{z,k}=z_k-{\widehat{\mu }}_{z,k}$$

(21)

where $z_k$ represents the actual measured value.

Set the threshold parameter:

$$\delta =\eta \sqrt{{{\stackrel{-}{P}}^{jj}}_{z,k}}$$

(22)

where $\eta ={\eta }_1+{\eta }_2‖e_{z,k}‖$, ${\eta }_1$ and ${\eta }_2$ represent the tuning parameters and ${{\stackrel{-}{P}}^{jj}}_{z,k}$ represents the diagonal element of the state estimation covariance matrix.

Calculate Huber weights for M estimation.

$$M(j,j)=\left\{\begin{array}{l}{\displaystyle \frac{\delta }{e_{z,k}}}, |e_{z,k}|\ge \delta \\ 1, |e_{z,k}|<\delta \end{array}\right.$$

(23)

7.

Measurement correction.

Updated estimates:

$${\widehat{\mu }}_k={\widehat{\mu }}_{x,k}+M{K}_k{e}_{z,k}$$

(24)

Update covariates:

$$P_k={\stackrel{-}{P}}_{x,k}-M{K}_k{\stackrel{-}{P}}_{z,k}{K}_k^T$$

(25)

The robustness of the UKF is enhanced by introducing Huber-M estimation to deal with outliers in the observed data.

4.2. Robust MPC Upper Layer Controller

The state space equations for longitudinal kinematics have been obtained above. The workshop distance, workshop distance error, relative velocity between vehicle i-1 and vehicle i, and the velocity and acceleration of vehicle i are chosen as the outputs of the state equations of the system:

$$y(k)=[d_{i-1,i},\delta d_{i-1,i},v_{i-1,i},v_i,a_i{]}^T$$

(26)

Thus, the system state output equation is obtained:

$$y(k)=D_{o}x(k)+E$$

(27)

Its coefficient matrix is $D_o=\left[\begin{array}{cccc}1& 0& 0& 0\\ 1& 0& -{\tau }_h& 0\\ 0& 1& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right]$, $E={\left[\begin{array}{ccccc}0& d_{min}& 0& 0& 0\end{array}\right]}^T$.

Introduce control quantities into the state equation to construct new state quantities:

$$X(k)=\left[\begin{array}{c}x(k)\\ u(k-1)\end{array}\right]$$

(28)

Derive the new state space equation as:

$$X(k+1)=Ax(k)+B\Delta u(k)+Cw(k)$$

(29)

where $A=\left[\begin{array}{cc}{A}_o& B_o\\ 0_{Nu\times Nx}& I_{Nu\times Nu}\end{array}\right]$, $B=\left[\begin{array}{c}B\\ I_{Nu\times Nu}\end{array}\right]$, $C=\left[\begin{array}{c}C\\ 0_{Nu\times 1}\end{array}\right]$,$\mathsf{\Delta }u\left(k\right)$ is the increment of the control quantity, $Nx$ is the dimension of the state quantity, and $Nu$ is the dimension of the control quantity.

The new system state output equation is derived as:

$$Y(k)=DX(k)+E$$

(30)

where $D=\left[D_0,0_{Ny\times 1}\right]$,$Ny$ represents the dimension of the output quantity.

4.2.1. Feedback Correction Compensation

Vehicles in the process of modeling, due to the interference of the external environment and the influence of the vehicle parameter error and other reasons, the established model has a certain mismatch, so the introduction of feedback correction mechanism to compensate for the prediction error of the following model, improve the model’s prediction accuracy of the state of the following system, and improve the robustness of the controller.

Define the state-volume error between the state-volume and the state-prediction at moment k:

$$e(k)=X(k)-X_p(k|k-1)$$

(31)

where $X(k)$ represents the actual state of vehicle i at moment k, and $X_p(k|k-1)$ is the state prediction of the system for moment k at moment k-1.

At moment k-1, the system predicts the amount of state at the moment:

$$X_p(k|k-1)=AX(k-1)+B\Delta u(k-1)+Cw(k-1)$$

(32)

The state prediction is corrected at moment k by the state quantity error:

$$X_p(k+1|k)=AX(k)+B\Delta u(k)+Cw(k)+Fe(k)$$

(33)

where $F=diag(f_{d_{i-1,i}},f_{v_{i-1,i}},f_{v_i},f_{a_i},f_{\Delta u})$ represents the correction matrix, and the value range of each element in F is (0, 1).

4.2.2. Predictive Model Derivation

Assuming that the prediction time domain of the system is $N_p$ and the control time domain is $N_c$, and that the prediction time domain and the control time domain satisfy $N_p\ge N_c$. Due to the introduction of the feedback correction mechanism, based on the discrete longitudinal motion state space Equation (33), a new model prediction equation of state can be obtained:

$$Y_p(k+N_p|k)=M_{DA}X(k)+M_{DB}\Delta U(k+N_c|k)+M_{DC}W(k+N_p|k)+M_{DF}e(k)+M_E$$

(34)

where $Y_p(k+N_p|k)$, $W(k+N_p|k)$ represent the output quantity matrix and disturbance quantity matrix of the system in the prediction time domain, and $\Delta U(k+N_c|k)$ represents the control quantity matrix of the system in the control time domain, as shown below:

$Y_p(k+N_p|k)=\left[\begin{array}{c}y(k+1|k)\\ y(k+2|k)\\ ⋮\\ y(k+N_p|k)\end{array}\right]$, $W(k+N_p)=\left[\begin{array}{c}w(k)\\ w(k+1)\\ ⋮\\ w(k+N_p-1)\end{array}\right]$, $\Delta U(k+N_c|k)=\left[\begin{array}{c}\Delta u(k|k)\\ \Delta u(k+1|k)\\ ⋮\\ \Delta u(k+N_c-1|k)\end{array}\right]$.

$M_{DA}$, $M_{DC}$, $M_{DF}$, $M_E$ is the prediction matrix and $M_{DB}$ is the control matrix as shown below:

$M_{DA}=\left[\begin{array}{c}DA\\ D{A}^2\\ ⋮\\ D{A}^{N_P}\end{array}\right]$, $M_{DC}=\left[\begin{array}{cccc}DC& 0& \cdots & 0\\ DAC& C& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ D{A}^{N_P-1}C& D{A}^{N_P-2}C& \cdots & DC\end{array}\right]$, $M_{DF}=\left[\begin{array}{c}DF\\ D{F}^2\\ ⋮\\ D{F}^{N_P}\end{array}\right]$, $M_E=\left[\begin{array}{c}E\\ E^2\\ ⋮\\ E^{N_P}\end{array}\right]$, $M_{DB}=\left[\begin{array}{cccc}DB& 0& \cdots & 0\\ DAB& DB& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ D{A}^{N_P-1}B& D{A}^{N_P-2}B& \cdots & {\sum }_{j=0}^{N_p-N_c}D{A}^{j}B\end{array}\right]$.

4.2.3. Objective Function

The exponential decay function is introduced as the target function reference trajectory to ensure that the reference trajectory changes more smoothly close to the target value, and the target reference trajectory is:

$$y_{ref}(k+i)={\Lambda }^{i}y(k)$$

(35)

Where $\Lambda =diag({\Lambda }_{d_{i-1,i}},{\Lambda }_{\delta d_{i-1,i}},{\Lambda }_{v_{i-1,i}},{\Lambda }_{v_i},{\Lambda }_{a_i})$, $\Lambda$ represents the coefficient matrix of the reference trajectory, and the value range of each element is (0,1).

The prediction matrix for the target reference trajectory is:

$$Y_{ref}(k+N_p|k)=\left[\begin{array}{c}{y}_{ref}(k+1|k)\\ y_{ref}(k+2|k)\\ ⋮\\ y_{ref}(k+N_p|k)\end{array}\right]$$

(36)

The objective function is designed in terms of followability and comfort requirements:

$$J=\sum _{i=1}^{N_p}{‖Y_p\left(k\right)-Y_{ref}\left(k\right)‖}_Q^2+\sum _{i=1}^{N_c}{‖\Delta U\left(k\right)‖}_R^2$$

(37)

Where $Q=diag(q_{d_{i-1,i}},q_{\delta d_{i-1,i}},q_{v_{i-1,i}},q_{v_i},q_{a_i})$, Q represents the weight matrix of the output quantity, R represents the weight matrix of the control quantity increment.

4.2.4. System State Constraints

In order to satisfy the requirements of traveling safety, following, and ride comfort of the platoon longitudinal traveling system, the MPC algorithm is constrained to be able to satisfy multiple constraints at the same time by imposing constraints on the output, control, and incremental control quantities of the MPC algorithm.

The constraints to be satisfied during the control of the MPC algorithm are as follows:

$$\begin{array}{c}\left\{\begin{array}{c}{d}_{min}\le d_{i-1,i}\left(k\right)\le d_{max}\\ \delta d_{min}\le \delta d_{i-1,i}(k)\le \delta d_{max}\\ \delta v_{min}\le v_{i-1,i}(k)\le \delta v_{max}\\ v_{min}\le v_i(k)\le v_{max}\\ a_{min}\le a_i(k)\le a_{max}\\ u_{min}\le u(k)\le u_{max}\\ \Delta u_{min}\le \Delta u(k)\le \Delta u_{max}\end{array}\right.\end{array}$$

(38)

where the parameters with min and max in the subscripts are the upper and lower bounds of the system constraints, respectively.

Into matrix form:

$$\left[\begin{array}{c}{Y}_{min}\\ u_{min}\\ \Delta u_{min}\end{array}\right]\le \left[\begin{array}{c}Y\\ u\\ \Delta u\end{array}\right]\le \left[\begin{array}{c}{Y}_{max}\\ u_{max}\\ \Delta u_{max}\end{array}\right]$$

(39)

where $Y_{min}=\left[\begin{array}{c}{d}_{min}\\ \delta d_{min}\\ \delta v_{min}\\ v_{min}\\ a_{min}\end{array}\right]$, $Y_{max}=\left[\begin{array}{c}{d}_{max}\\ \delta d_{max}\\ \delta v_{max}\\ v_{max}\\ a_{max}\end{array}\right]$.

For solving the problem of no solution in the rolling optimization process of MPC algorithm, the hard constraints are relaxed by introducing relaxation factors and relaxation coefficients to extend the feasible domain of the solution, as shown below:

$$\left[\begin{array}{c}{Y}_{min}\\ u_{min}\\ {\mathsf{\Delta }u}_{min}\end{array}\right]+\left[\begin{array}{c}{\varphi }_{min}^Y\\ {\varphi }_{min}^u\\ {\varphi }_{min}^{\mathsf{\Delta }u}\end{array}\right]\epsilon \le \left[\begin{array}{c}Y\\ u\\ \Delta u\end{array}\right]\le \left[\begin{array}{c}{Y}_{max}\\ u_{max}\\ {\mathsf{\Delta }u}_{max}\end{array}\right]+\left[\begin{array}{c}{\varphi }_{max}^Y\\ {\varphi }_{max}^u\\ {\varphi }_{max}^{\mathsf{\Delta }u}\end{array}\right]\epsilon$$

(40)

where $\epsilon$ represents the relaxation factor matrix, ${\varphi }_{min}^Y$, ${\varphi }_{max}^Y$, ${\varphi }_{min}^u$, ${\varphi }_{max}^u$, ${\varphi }_{min}^{\Delta u}$, ${\varphi }_{max}^{\Delta u}$ is the relaxation coefficient matrix as shown below:

$$\left[\begin{array}{c}{\varphi }_{min}^Y\\ {\varphi }_{min}^u\\ {\varphi }_{min}^{\mathsf{\Delta }u}\end{array}\right]=diag({\phi }_{min}^d,{\phi }_{min}^{\delta d},{\phi }_{min}^{\delta v},{\phi }_{min}^v,{\phi }_{min}^a,{\phi }_{min}^u,{\phi }_{min}^{\mathsf{\Delta }u}),$$

$$\left[\begin{array}{c}{\varphi }_{max}^Y\\ {\varphi }_{max}^u\\ {\varphi }_{max}^{\mathsf{\Delta }u}\end{array}\right]=diag({\phi }_{max}^d,{\phi }_{max}^{\delta d},{\phi }_{max}^{\delta v},{\phi }_{max}^v,{\phi }_{max}^a,{\phi }_{max}^u,{\phi }_{max}^{\mathsf{\Delta }u}),$$

$${\epsilon =\left[\begin{array}{ccccccc}{\epsilon }_d& {\epsilon }_{\delta d}& {\epsilon }_{\delta v}& {\epsilon }_v& {\epsilon }_a& {\epsilon }_u& {\epsilon }_{\Delta u}\end{array}\right]}^T.$$

4.2.5. Optimization Problem Solving

To avoid the problem of constraint failure due to relaxation factors, a quadratic penalty term for the relaxation factors is introduced into the optimization objective function to limit the range of constraints due to relaxation, then the performance cost function is:

$$J_T=J+{\epsilon }^T\rho \epsilon$$

(41)

Where $\rho =diag({\rho }_d,{\rho }_{\delta d},{\rho }_{\delta v},{\rho }_v,{\rho }_a,{\rho }_u,{\rho }_{\Delta u})$, $\rho$ represents the matrix of penalty coefficients for the relaxation factors.

Expanding and ignoring terms not related to $\mathsf{\Delta }U$, the collation yields:

$$J_T={\displaystyle \frac{1}{2}}\Delta U^{T}H\Delta U+f^T\Delta U+{\epsilon }^T\rho \epsilon$$

(42)

where $\left\{\begin{array}{l}H=2({M_{DB}}^T\tilde{Q}{M}_{DB}+\tilde{R})\\ f^T=-2{E_e}^T\tilde{Q}{M}_{DB}\end{array}\right.$, $\tilde{Q}=diag(Q,Q,\cdots Q)$, $\tilde{R}=diag(R,R,\cdots R)$, $E_e=Y_{ref}-Y_p$.

Transform the performance cost function treatment into a quadratic programming online solution problem with constraints in the form shown below:

$$\left\{\begin{array}{l}\mathit{min}\left\{{\displaystyle \frac{1}{2}}{U}_{\rho }^T{H}_{\rho }{U}_{\rho }+f_{\rho }^T{U}_{\rho }\right\}\\ s.t.A_u{U}_{\rho }^T\le b_u\end{array}\right.$$

(43)

where $U_{\rho }=\left[\begin{array}{c}\Delta U\\ \epsilon \end{array}\right]$, $H_{\rho }=\left[\begin{array}{cc}H& \\ & \rho \end{array}\right]$, ${f_{\rho }}^T=\left[\begin{array}{cc}{f}^T& 0_{1\times (Ny+Nu+1)}\end{array}\right]$.

Substituting the model predicted state equation and the control quantity $U=U_{k-1}+A_u\Delta U$ into the system state constraint (40) yields:

$$A_u=\left[\begin{array}{cc}{M}_{DB}& -M_{max}^{\varphi Y}\\ -M_{DB}& M_{min}^{\varphi Y}\\ A_u& -M_{max}^{\varphi u}\\ -A_u& M_{min}^{\varphi u}\\ I_{Nc\times Nc}& -M_{max}^{\varphi \Delta u}\\ I_{Nc\times Nc}& M_{min}^{\varphi \Delta u}\end{array}\right], b_u=\left[\begin{array}{c}{M}_{max}^{Y_p}-Y_p\\ -M_{min}^{Y_p}+Y_p\\ U_{max}-U_{k-1}\\ -U_{min}+U_{k-1}\\ \Delta U_{max}\\ -\Delta U_{min}\end{array}\right].$$

where $M_{max}^{\varphi Y}$, $M_{min}^{\varphi Y}$, $M_{max}^{\varphi u}$, $M_{min}^{\varphi u}$, $M_{max}^{\varphi \Delta u}$, $M_{min}^{\varphi \Delta u}$ represent the prediction matrix of the relaxation coefficient matrix, and $M_{max}^Y$, $M_{min}^Y$, $U_{max}$, $U_{min}$, $\Delta U_{max}$, $\Delta U_{min}$ represent the prediction matrix of the upper and lower bounds of the output, control, and control increments, respectively.

The upper controller decides the desired acceleration by means of a robust MPC control algorithm.

4.3. Lower Level Controller

4.3.1. Inverse Dynamics Model

According to the vehicle longitudinal kinematics model, the vehicle longitudinal acceleration and deceleration equations of motion are obtained.

Acceleration equations:

$$m_i\delta a_i=F_t-F_w-F_f-F_i$$

(44)

where $m_i$ represents the vehicle mass of vehicle i, $\delta$ represents the vehicle rotating mass conversion factor, $F_t$ represents the driving force, $F_w$ represents the air resistance, $F_f$ represents the rolling resistance and $F_i$ represents the ramp resistance.

Deceleration equations:

$$m_i\delta a_i=F_t-F_b-F_w-F_f-F_i$$

(45)

where $F_b$ represents the desired braking force.

In drive mode, the desired torque of the drive motor is derived from Equation (44):

$$T_{des}=[m_{i}g(f\mathit{cos}\theta +\mathit{sin}\theta )+{\displaystyle \frac{C_{D}A{v}_i^2}{21.15}}+m_i\delta a_i]{\displaystyle \frac{r_i}{i_t{\eta }_t}}$$

(46)

where $g$ represents the acceleration of gravity, $f$ represents the rolling resistance coefficient, $\theta$ represents the ramp angle, is the air resistance coefficient, $C_D$ represents the windward area, $r_i$ represents the rolling radius of the wheels, $i_t$ represents the main gear transmission ratio, ${\eta }_t$ represents the transmission efficiency of the power transmission system.

In braking mode, the desired brake master cylinder pressure is derived from Equation (45):

$$P_{des}=|-m_{i}g(f\mathit{cos}\theta +\mathit{sin}\theta )-{\displaystyle \frac{C_{D}A{v}_i^2}{21.15}}-m_i\delta a_i|{\displaystyle \frac{1}{K_b}}$$

(47)

where $K_b$ represents the coefficient of proportionality between braking force and brake master cylinder pressure.

4.3.2. PID Lower Controller

As shown in Figure 2, the lower controller part adopts the PID control algorithm, utilizes the form of feed-forward plus feedback, and controls the driving torque and braking pressure based on the inverse longitudinal dynamics model of the vehicle by means of the drive-brake switching strategy, which tracks the desired acceleration obtained from the upper controller.

5. Simulation Analysis

For validating the designed robust UKF-MPC longitudinal controller, a joint simulation platform based on CarSim/Simulink is built to simulate and analyze the longitudinal control of vehicle platoon.

The simulation scenario with no communication delay and no data loss and the simulation scenario with a random communication delay of 10~100 ms and a communication data loss rate of 50% are set up in Simulink, and the fleet of vehicles containing one leading vehicle and three following vehicles is set up in CarSim, with the initial speed of the four vehicles being 0 m/s and the initial inter-vehicle spacing of 3 m. The main parameters of the vehicles used are shown in Table 1.

Set the target speed of the leading vehicle, including acceleration, uniform speed, deceleration of three kinds of driving conditions, the specific parameters are set as shown in Table 2.

The robust UKF-MPC longitudinal controller of this paper is compared and analyzed with the MPC longitudinal controller and the improved UKF-MPC longitudinal controller based on UKF in the scenarios with and without communication delay and data packet loss problems, and the simulation results are shown in Figure 3, Figure 4, Figure 5, Figure 6, Figure 7 and Figure 8.

Under the condition of no communication delay and packet loss, as shown in Figure 3 and Figure 4, the speed and acceleration under the three controllers can stably follow the speed and acceleration of the front vehicle to ensure the stability of the longitudinal motion of the vehicle platoon; however, when comparing the inter-vehicle spacing errors of the vehicles controlled by the three controllers, the inter-vehicle spacing error under the MPC longitudinal controller is about 0.1 m, that under the UKF-MPC longitudinal controller is about is 0.3 m, and the workshop distance error under the robust UKF-MPC longitudinal controller is about 0.04 m, as shown in Figure 5. It can be seen that the robust UKF-MPC longitudinal controller in this paper improves the control accuracy of the longitudinal motion of the vehicle platoon compared to the MPC longitudinal controller and the UKF-MPC longitudinal controller under the condition of no communication delay and packet loss.

Under the condition of communication delay and packet loss, as shown in Figure 6, the speeds under the three controllers can basically follow the speed of the front vehicle, but the speed of the vehicle controlled by the MPC longitudinal controller jerks at higher speeds; as shown in Figure 7, the acceleration under the MPC controller oscillates violently, the acceleration under the UKF-MPC longitudinal controller jerks slightly and the acceleration under the robust UKF-MPC longitudinal controller jerks less and is basically stable. As shown in Figure 8, the workshop distance error under the MPC longitudinal controller is about 3 m, and the oscillation is violent, the workshop distance error under the UKF-MPC longitudinal controller is about 0.35 m, and there is a jitter when the speed is higher, and the workshop distance error under the robust UKF-MPC longitudinal controller is about 0.045 m, and it is basically stable. It can be seen that the MPC control algorithm is not effective under the condition of having communication random delay and packet loss, while the UKF-MPC control algorithm and the robust UKF-MPC control algorithm can effectively solve the problem of communication delay and packet loss, but the robust UKF-MPC control algorithm is more stable and has higher accuracy.

6. Conclusions

In this paper, a robust UKF-MPC based longitudinal controller for electric vehicle platoon is designed and validated by simulation comparison with MPC longitudinal controller and UKF-MPC longitudinal controller for continuous acceleration and deceleration conditions with and without communication problems.

The simulation results show that under the conditions of no communication delay and packet loss, the proposed robust UKF-MPC controller reduces the vehicle spacing error compared with the MPC and UKF-MPC longitudinal controllers, and improves the control accuracy. Under the conditions with communication delay and packet loss, the MPC longitudinal controller’s performance significantly degrades, which seriously affects the traveling safety and ride comfort. However, the robust UKF-MPC longitudinal controller still maintains a small inter-vehicle distance error and good stability.