Neural Network and Extended State Observer-Based Model Predictive Control for Smooth Braking at Preset Points in Autonomous Vehicles

Abstract

In this paper, we explore the problem of smooth braking at preset points in autonomous vehicles using model predictive control (MPC) with a receding horizon extended state observer (RHESO) and a neural network (NN). An NN-based modeling method is proposed to intuitively describe the relationship between vehicle speed and the vehicle controllers (brake and throttle), and establish a dynamic model of autonomous vehicles. A sufficient condition is put forward to guarantee the convergence of the proposed NN. Furthermore, a composite MPC strategy based on RHESO is designed, which optimizes a given cost function over the receding horizon while mitigating the effects of modeling inaccuracies and disturbances. Additionally, easily verifiable conditions are provided to ensure the autonomous driving vehicles’ uniform boundedness. Numerically illustrative examples are given to demonstrate the effectiveness of the proposed approach.

1. Introduction

Generally, non-stationary braking in vehicles poses major safety hazards, which may cause vehicle imbalance and accidents [1], resulting in uncomfortable riding experiences for passengers [2], damaging brake system components, increasing repair and maintenance costs [3], and even posing a threat to pedestrians on the road [4].

Unlike manual braking, the emergence of autonomous vehicles (AVs) makes it easier for vehicles to brake smoothly and accurately. Accuracy and smoothness are the two key performance criteria of smooth braking in AVs, meaning low tracking errors can be achieved without aggressive braking/throttle operations. Furthermore, smooth and precise braking in AVs has potential applications in many fields, such as bus passenger pick-up and drop-off [5], truck loading and unloading of goods at preset points [6], and passenger elevator docking with airplanes [7]. In these application scenarios, AVs must stop smoothly and precisely at preset points to ensure comfortable riding experiences for passengers and avoid collisions. Moreover, smooth braking at preset points in AVs significantly improves both parking efficiency and parking space utilization.

Smooth braking at preset points in AVs is a typical position and speed (PS) tracking control problem, which is also a dynamic control problem [8,9]. Existing PS tracking methods can be broadly categorized into three groups. The first group of methods does not consider vehicle dynamics or disturbance models. Control commands are solely generated based on tracking errors, employing designs such as proportional–integral–derivative (PID) controllers and their variants [10,11]. The disadvantage is that control precision heavily depends on adjustable parameters with high conservation and poor portability. The second method type relies on precise dynamic models that overlook disturbances to better predict PS information [12]. In [13], the authors analyzed a smooth braking torque control strategy for a brushless DC motor. In [14], a preview of a servo-loop speed control algorithm was proposed to achieve smooth, accurate, and computationally inexpensive speed tracking for AVs. Furthermore, cooperative braking control was introduced to ensure that AVs maintained safe spacing and stopped rapidly, smoothly, and accurately at the designated target positions with zero velocity [15]. These studies range from torque control based on brushless DC motor models to designing smooth and accurate braking methods using brushless DC motor models to designing smooth and accurate braking methods based on vehicle kinematic models, achieving various conditions for smooth and accurate braking of autonomous vehicles from hardware to software levels. However, precise dynamic models of AVs are difficult to obtain, and the real-time performance of control calculations using these models is often poor. A common drawback of the first two method types is their inability to handle the inherent constraints of AVs, such as speed limits, front-wheel deviation limits, and control input constraints

Different from the first two groups of methods, the third type of control strategy is based on model predictive control (MPC), which effectively handles various constraints. By minimizing the difference between the actual and predicted PS, the optimal control signal can be generated through repeated online optimization [16,17]. Under the MPC framework, there are schemes that use dynamic models of AVs [18,19], as well as schemes that use kinematic models of AVs [20,21]. Note that schemes using dynamic models for tracking PS constitute a servo control problem [22,23], with the advantage of directly tracking speed as a state variable. However, the disadvantage is that accurate model parameters are hard to obtain, and the real-time solution of the control signal cannot be guaranteed due to strong nonlinearity. In contrast, the scheme using kinematic models of AVs essentially addresses a speed planning problem [24], with the advantage of easily obtaining precise model parameters and ensuring real-time control signal solutions; the disadvantage is that the planned speed must be tracked by the servo layer, which leads to delays in speed tracking.

To maintain the high tracking accuracy of servo control while reducing the optimization load, the best approach is to combine the strengths of both the kinematic and dynamic models to establish a new dynamic model that can effectively handle the servo control of PS. One of the most direct methods is to add the speed variable as a state variable into the kinematic model. Correspondingly, throttle and brake openings are used as control inputs to manage the AVs’ speed servo control. However, due to unknown nonlinearities and disturbances, the relationship between the throttle (or brake) opening and AVs’ speed is challenging to determine. As an excellent tool for approximating nonlinear functions, neural networks (NNs) can describe the relationship between AVs’ speed and the throttle (or brake) opening [25,26]. There is significant research on using NNs to approximate unknown nonlinear system models, but few studies have applied NNs to autonomous vehicle modeling.

In addition, distractions such as pedestrians in vehicles, uneven road surfaces, and the coupling between equipment inside the vehicle often affect the accurate tracking of AVs’ PS. To minimize the impact of these disturbances, it is necessary to design disturbance observers for autonomous vehicles. The existing observers are primarily divided into two categories. The first type of observer is designed using the dynamic characteristics of the disturbance, which usually requires prior knowledge of the dynamic model parameters of the disturbance [27,28,29,30]. Another type of disturbance observer is the extended state observer (ESO); it can estimate the disturbance with a bounded rate of change and it has a wider application range [31,32,33]. However, these studies mainly focus on the final observation accuracy of the disturbance observer, while neglecting the stationarity of the disturbance observation process (i.e., the disturbance observation process may have a significant overshoot). One of the goals of this paper is to comprehensively consider the disturbance estimation process and smooth braking, as well as realize smooth and accurate braking at preset points in AVs while canceling the disturbance effects.

This paper aims to address the problem of smooth braking at preset points in AVs with a receding horizon extended state observer (RHESO)-based composite MPC strategy. The main contributions of this paper are as follows:

(1)

A novel extended kinematic model of AVs is established using a proposed NN such that the relationship between the AVs’ speed and throttle (and brake) opening is more clearly described;

(2)

A novel composite MPC strategy based on the RHESO is put forward, which is optimized over the receding horizon while eliminating the effects of the model’s inaccuracy and disturbance;

(3)

Easily verifiable sufficient conditions are established to ensure the recursive feasibility of the MPC scheme and the stability of the closed-loop system.

Finally, the notation used in this article are presented in

Table 1. Notation used in this paper.

$\mathbb{R}$	Set of Real Numbers
${\mathbb{Z}}_{\ge 0}$	Set of nonnegative integers
${\mathbb{Z}}^{+}$	Set of positive integers
$\mathbb{Z}[a,b]$	The positive integer set $\{a,a+1,\cdots ,b\}$
$A\in {\mathbb{R}}^{m\times n}$	The dimensions of A are $m\times n$
$A>0$	A is a positive definite matrix
$A<0$	A is a negative definite matrix
I	The unit matrix with appropriate dimensions
$\parallel ·\parallel$	The Euclidean norm
$\mathbf{0}$	The matrix full of 0 with appropriate dimensions
$\overline{\lambda }\left(A\right)$	The maximum eigenvalue of A
$\underline{\lambda }\left(A\right)$	The minimum eigenvalue of A
$x(k+s|k)$	The prediction of x at the time instant $k+s$ from the time instant k

.

2. Problem Formulations and Preliminaries

2.1. System Model

A general vehicle kinematic model is given as follows:

$$\begin{array}{c}\hfill \dot{{\xi }_c}\left(t\right)=\left[\begin{array}{c}\dot{x}\left(t\right)\\ \dot{y}\left(t\right)\\ \dot{\theta }\left(t\right)\end{array}\right]=\left[\begin{array}{c}v\left(t\right)cos\left(\theta \right(t\left)\right)\\ v\left(t\right)sin\left(\theta \right(t\left)\right)\\ \omega \left(t\right)\end{array}\right]\end{array}$$

(1)

where ${\xi }_c\left(t\right)\in {\mathbb{R}}^3$ indicates the state vector, $x\left(t\right)\in \mathbb{R}$ and $y\left(t\right)\in \mathbb{R}$ are the reference coordinates on the rear axle center in a geodetic coordinate system, $\theta \left(t\right)\in \mathbb{R}$ denotes the heading angle, $u_c\left(t\right)={\left[\omega \left(t\right)v\left(t\right)\right]}^T\in {\mathbb{R}}^2$ is the control input with the linear speed $v\left(t\right)\in \mathbb{R}$, the angular speed is $\omega \left(t\right)\triangleq v\left(t\right)\frac{tan\left(\delta \right(t\left)\right)}{L}\in \mathbb{R}$, where $L\in \mathbb{R}$ is the vehicle wheelbase, and $\delta \left(t\right)\in \mathbb{R}$ is the wheel deflection angle. The vehicle model is presented in Figure 1.

Generally, the inner and outer ring structures are used to control the intelligent driving vehicle. The function of the outer loop is to generate the desired longitudinal speed $v\left(t\right)$ and lateral speed $\omega \left(t\right)$ according to the desired trajectory, while the function of the inner loop is to obtain the real control input signal, i.e., the throttle opening (or brake opening) and steering wheel angle, by tracking the desired longitudinal speed $v\left(t\right)$ and lateral speed $\omega \left(t\right)$, respectively.

Note that if the desired longitudinal speed $v\left(t\right)$ and lateral speed $\omega \left(t\right)$ change rapidly, the tracking of $v\left(t\right)$ and $\omega \left(t\right)$ in the inner loop may experience hysteresis and inaccuracies. Under such conditions, it becomes challenging to achieve smooth braking at preset points for AVs. Since the variation in longitudinal velocity during braking typically exceeds that of lateral velocity, it can be assumed that $\omega \left(t\right)=0$ and $\theta \left(t\right)=\overline{\theta }$ with $\overline{\theta }$ being a known parameter. Considering these characteristics, this paper proposes a neural network (NN)-based longitudinal kinematic model, as follows:

$$\begin{array}{c}\hfill \dot{\xi }\left(t\right)=\left[\begin{array}{c}\dot{x}\left(t\right)\\ \dot{y}\left(t\right)\\ \dot{v}\left(t\right)\end{array}\right]=\left[\begin{array}{c}v\left(t\right)cos\left(\overline{\theta }\right)\\ v\left(t\right)sin\left(\overline{\theta }\right)\\ W^T\sigma \left(u_{\mathrm{tb}}\left(t\right)\right)+\epsilon \left(t\right)\end{array}\right],\end{array}$$

(2)

where the extended dynamic model $\dot{v}\left(t\right)=W^T\sigma \left(u_{\mathrm{tb}}\left(t\right)\right)+\epsilon \left(t\right)$ is based on an NN to be designed; $W\in {\mathbb{R}}^{3\times 1}$ is an unknown ideal matrix consisting of weight and bias; $\sigma (·)$ is a vector of smooth activation functions; $\epsilon \left(t\right)$ is the external disturbance satisfying $\parallel \epsilon \left(t\right)\parallel \le \overline{\epsilon }$ and $\parallel \dot{\epsilon }\left(t\right)\parallel \le \tilde{\epsilon }$ with $\overline{\epsilon }$, and $\tilde{\epsilon }$ being known positive constants; $u_{\mathrm{tb}}\left(t\right)\triangleq {\left[v\left(t\right)u_{\mathrm{o}}\left(t\right)1\right]}^T$ is the input signal of the NN, and $u_{\mathrm{o}}\left(t\right)$ is the longitudinal speed control variable, where

$$\begin{array}{c}\hfill u_{\mathrm{o}}\left(t\right)=\left\{\begin{array}{cc}{u}_{\mathrm{to}}\left(t\right),& \mathrm{if}{u}_{\mathrm{to}}\left(t\right)\in (0,1],u_{\mathrm{bo}}\left(t\right)=0\\ u_{\mathrm{bo}}\left(t\right),& \mathrm{if}{u}_{\mathrm{bo}}\left(t\right)\in (0,9],u_{\mathrm{to}}\left(t\right)=0\\ 0,& \mathrm{if}{u}_{\mathrm{to}}\left(t\right)=0,u_{\mathrm{bo}}\left(t\right)=0\end{array}\right.\end{array}$$

with $u_{\mathrm{to}}\left(t\right)$ and $u_{\mathrm{bo}}\left(t\right)$ being the throttle opening and brake opening, respectively.

As stated in [34], it is assumed that there exists a known constant, ${ϵ}_W$, such that $\parallel W\parallel \le {ϵ}_W$. In addition, according to [35], the gradient of $\sigma (·)$ is assumed to satisfy $\parallel \nabla \sigma (·)\parallel \le {ϵ}_{\sigma }$ with ${ϵ}_{\sigma }$ being a positive constant. For simplicity, we denote $\sigma \left(t\right)\triangleq \sigma \left(u_{\mathrm{tb}}\left(t\right)\right)$.

Note that the AV is driven by the control input $u\left(t\right)\triangleq u_o\left(t\right)$ in the form of a segment constant in each sampling interval $[iT,(i+1\left)T\right)$, $i\in {\mathbb{Z}}_{\ge 0}$, with  sampling period T (in seconds). Based on (2), a discrete-time AV kinematic model is given as follows:

$$\begin{array}{c}\hfill \xi (k+1)=\left[\begin{array}{c}x\left(k\right)\\ y\left(k\right)\\ v\left(k\right)\end{array}\right]+\left[\begin{array}{c}Tv\left(k\right)cos\left(\overline{\theta }\right)\\ Tv\left(k\right)sin\left(\overline{\theta }\right)\\ T{W}^T\sigma \left(k\right)+T\epsilon \left(k\right)\end{array}\right].\end{array}$$

(3)

In the following, the control framework diagram of this paper is presented in Figure 2. In Figure 2, the control framework is divided into offline and online parts. The offline part is used to collect data on the speed and acceleration of AVs under different throttle and brake openings and to train NN parameters using the collected data to prepare for subsequent system control. The online part is a closed-loop control loop, including the RHESO and MPC that need to be designed later.

2.2. Estimation of the Unknown NN Parameter

The ideal NN parameter W is an unknown matrix that exists in theory, but an estimated matrix ${\widehat{W}}^{*}$ that is close enough to W can be obtained through estimation methods. Therefore, an adaptive estimator is designed to obtain the estimation ${\widehat{W}}^{*}$.

Consider the following adaptive estimator:

$$\begin{array}{c}\hfill \widehat{v}(k+1)=v\left(k\right)+T{\widehat{W}}^T\left(k\right)\sigma \left(k\right),\end{array}$$

where $\widehat{v}\left(k\right)$ is the estimation of $v\left(k\right)$; $\widehat{W}\left(k\right)$ is the estimation of W at time instant k. Subtracting $v(k+1)$ from $\widehat{v}(k+1)$ leads to the following:

$$\begin{array}{c}\hfill \widehat{v}(k+1)-v(k+1)=T({\widehat{W}}^T\left(k\right)-W^T)\sigma \left(k\right)-T\epsilon \left(k\right),\end{array}$$

which yields the following:

$$\begin{array}{c}\hfill \frac{1}{T}(\widehat{v}(k+1)-v(k+1))=({\widehat{W}}^T\left(k\right)-W^T)\sigma \left(k\right)-\epsilon \left(k\right).\end{array}$$

(4)

Denote

$$\begin{array}{c}\hfill e(k+1)\triangleq \frac{1}{T}(\widehat{v}(k+1)-v(k+1)),\tilde{W}\left(k\right)\triangleq \widehat{W}\left(k\right)-W,\end{array}$$

then, (4) can be rewritten as follows:

$$\begin{array}{c}\hfill e(k+1)={\tilde{W}}^T\left(k\right)\sigma \left(k\right)-\epsilon \left(k\right).\end{array}$$

To minimize $E\left(k\right)=\frac{1}{2}{e}^T(k+1)e(k+1)$, using the Levenberg–Marquardt algorithm [36], the tuning law of $\widehat{W}\left(k\right)$ is designed as follows:

$$\begin{array}{cc}\hfill \widehat{W}(k+1)=& \widehat{W}\left(k\right)-\frac{\Lambda }{1+{\sigma }^T\left(k\right)\sigma \left(k\right)}{\left(\frac{\partial E\left(k\right)}{\partial \widehat{W}\left(k\right)}\right)}^T=\widehat{W}\left(k\right)-\frac{\Lambda \sigma \left(k\right)e^T(k+1)}{1+{\sigma }^T\left(k\right)\sigma \left(k\right)},\hfill \end{array}$$

(5)

where $\Lambda$ is a positive constant to be designed.

In the following, a theorem is given to demonstrate that the proposed tuning law (5) can make the estimation $\widehat{W}\left(k\right)$ converge to a neighborhood of ideal value W.

Theorem 1. 

If $\sigma \left(k\right)$ in (5) is persistently exciting and $\mathbf{0}<\Lambda <1$ holds, then $\tilde{W}\left(k\right)$ is uniformly bounded to the set $\mathbb{W}$ where $\mathbb{W}\triangleq \left\{\tilde{W}\right|\parallel \tilde{W}{\parallel }^2\le \overline{\kappa }\}$ with the following:

$$\begin{array}{c}\hfill \overline{\kappa }\triangleq \frac{1}{\kappa }{\left(1+\frac{1}{\sqrt{1-\Lambda }}\right)}^2{\overline{\epsilon }}^2,\kappa \triangleq \underline{\lambda }\left(\sigma \left({\left[001\right]}^T\right){\sigma }^T\left({\left[001\right]}^T\right)\right).\end{array}$$

Proof. 

We select a Lyapunov candidate function, as follows:

$$\begin{array}{c}\hfill V\left(\tilde{W}\left(k\right)\right)={\tilde{W}}^T\left(k\right){\Lambda }^{-1}\tilde{W}\left(k\right),\end{array}$$

thus, we have the following:

$$\begin{array}{cc}\hfill & V\left(\tilde{W}(k+1)\right)-V\left(\tilde{W}\left(k\right)\right)\hfill \\ \hfill =& -\frac{e(k+1){\sigma }^T\left(k\right)\tilde{W}\left(k\right)}{1+{\sigma }^T\left(k\right)\sigma \left(k\right)}-\frac{{\tilde{W}}^T\left(k\right)\sigma \left(k\right)e^T(k+1)}{1+{\sigma }^T\left(k\right)\sigma \left(k\right)}+\frac{e(k+1){\sigma }^T\left(k\right)\Lambda \sigma \left(k\right)e^T(k+1)}{{\left(1+{\sigma }^T\left(k\right)\sigma \left(k\right)\right)}^2}\hfill \\ \hfill =& -\frac{2e(k+1)e^T(k+1)}{1+{\sigma }^T\left(k\right)\sigma \left(k\right)}-\frac{e(k+1){\epsilon }^T\left(k\right)}{1+{\sigma }^T\left(k\right)\sigma \left(k\right)}-\frac{\epsilon \left(k\right)e^T(k+1)}{1+{\sigma }^T\left(k\right)\sigma \left(k\right)}\hfill \\ \hfill & +\frac{{\sigma }^T\left(k\right)\Lambda \sigma \left(k\right)e(k+1)e^T(k+1)}{{\left(1+{\sigma }^T\left(k\right)\sigma \left(k\right)\right)}^2}\hfill \\ \hfill \le & -\frac{e(k+1)e^T(k+1)}{1+{\sigma }^T\left(k\right)\sigma \left(k\right)}+\frac{\epsilon \left(k\right){\epsilon }^T\left(k\right)}{1+{\sigma }^T\left(k\right)\sigma \left(k\right)}+\frac{{\sigma }^T\left(k\right)\Lambda \sigma \left(k\right)e(k+1)e^T(k+1)}{{\left(1+{\sigma }^T\left(k\right)\sigma \left(k\right)\right)}^2}\hfill \\ \hfill \le & -\left(1-\frac{{\sigma }^T\left(k\right)\Lambda \sigma \left(k\right)}{1+{\sigma }^T\left(k\right)\sigma \left(k\right)}\right)\frac{e(k+1)e^T(k+1)}{1+{\sigma }^T\left(k\right)\sigma \left(k\right)}+\frac{\epsilon \left(k\right){\epsilon }^T\left(k\right)}{1+{\sigma }^T\left(k\right)\sigma \left(k\right)}\hfill \\ \hfill \le & -(1-\Lambda )\frac{e(k+1)e^T(k+1)}{1+{\sigma }^T\left(k\right)\sigma \left(k\right)}+\frac{\epsilon \left(k\right){\epsilon }^T\left(k\right)}{1+{\sigma }^T\left(k\right)\sigma \left(k\right)}.\hfill \end{array}$$

(6)

Select a proper $\Lambda$, such that $0<\Lambda <1$. By letting $V\left(\tilde{W}(k+1)\right)-V\left(\tilde{W}\left(k\right)\right)<0$, we can obtain the following:

$$\begin{array}{c}\hfill \frac{\epsilon \left(k\right){\epsilon }^T\left(k\right)}{1+{\sigma }^T\left(k\right)\sigma \left(k\right)}<(1-\Lambda )\frac{e(k+1)e^T(k+1)}{1+{\sigma }^T\left(k\right)\sigma \left(k\right)},\end{array}$$

which indicates the following:

$$\begin{array}{c}\hfill {\epsilon }^T\left(k\right)\epsilon \left(k\right)<(1-\Lambda )e^T(k+1)e(k+1).\end{array}$$

(7)

This means that $V\left(\tilde{W}\left(k\right)\right)$ will remain convergent until condition (7) is not satisfied. Due to $e(k+1)\in \mathbb{R}$ and $\epsilon \left(k\right)\in \mathbb{R}$, we have the following:

$$\begin{array}{c}\hfill \frac{1}{1-\Lambda }{\epsilon }^2\left(k\right)<{({\tilde{W}}^T\left(k\right)\sigma \left(k\right)-\epsilon \left(k\right))}^2,\end{array}$$

which yields the following:

$$\begin{array}{c}\hfill {\left(1+\frac{1}{\sqrt{1-\Lambda }}\right)}^2{\epsilon }^2\left(k\right)<{\tilde{W}}^T\left(k\right)\sigma \left(k\right){\sigma }^T\left(k\right)\tilde{W}\left(k\right).\end{array}$$

(8)

To ensure that inequality (8) always holds, we require the following:

$$\begin{array}{c}\hfill {\left(1+\frac{1}{\sqrt{1-\Lambda }}\right)}^2{\overline{\epsilon }}^2<\underline{\lambda }\left(\sigma \left(k\right){\sigma }^T\left(k\right)\right){\tilde{W}}^T\left(k\right)\tilde{W}\left(k\right).\end{array}$$

(9)

Considering $\sigma \left(k\right)=\sigma \left(u_{\mathrm{tb}}\left(k\right)\right)$, we have the following:

$$\begin{array}{c}\hfill \kappa \triangleq \underline{\lambda }\left(\sigma \left(k\right){\sigma }^T\left(k\right)\right)=\underline{\lambda }\left(\sigma \left({\left[001\right]}^T\right){\sigma }^T\left({\left[001\right]}^T\right)\right).\end{array}$$

Then, (9) can be rewritten as follows:

$$\begin{array}{c}\hfill \frac{1}{\kappa }{\left(1+\frac{1}{\sqrt{1-\Lambda }}\right)}^2{\overline{\epsilon }}^2<{\tilde{W}}^T\left(k\right)\tilde{W}\left(k\right).\end{array}$$

Denote $\mathbb{W}\triangleq \left\{\tilde{W}\left(k\right)\right|{\parallel \tilde{W}\left(k\right)\parallel }^2\le \overline{\kappa }\}$ with $\overline{\kappa }\triangleq \frac{1}{\kappa }{\left(1+\frac{1}{\sqrt{1-\Lambda }}\right)}^2{\overline{\epsilon }}^2$. Based on the fact that $\sigma \left(k\right)$ is persistently exciting [37], it can be obtained that $\tilde{W}\left(k\right)$ would converge to the set $\mathbb{W}$. Therefore, $\tilde{W}\left(k\right)$ is uniformly bounded. This completes the proof.    □

From Theorem 1, it can be obtained that $\tilde{W}\left(k\right)$ will finally converge to the set $\mathbb{W}$. Therefore, if there exists $\tilde{W}\left(k_s\right)\in \mathbb{W}$ for $k_s\in {\mathbb{Z}}_{\ge 0}$, then $\tilde{W}(k_s+1)\in \mathbb{W}$ holds. Hence, it is easy to obtain the following:

$$\begin{array}{cc}\hfill \parallel \tilde{W}(k_s+1)-\tilde{W}\left(k_s\right)\parallel =& \parallel \widehat{W}(k_s+1)-\widehat{W}\left(k_s\right)\parallel \le \parallel \tilde{W}(k_s+1)\parallel + \parallel \tilde{W}\left(k_s\right)\parallel \le 2\sqrt{\overline{\kappa }}.\hfill \end{array}$$

(10)

Inequality (10) can be used as a termination condition of the tuning law in (5). When (10) is satisfied, ${\widehat{W}}^{*}=\widehat{W}\left(k_s\right)$ can be obtained.

In the following, the process of approximating AVs’ dynamics using a NN is given in Algorithm 1.

Algorithm 1 first collects data on changes in AV velocity and acceleration under different throttle and brake opening conditions. Secondly, these data are divided into inputs and outputs to train the NN. Thirdly, the NN structure, loss function, and activation function are designed. Finally, the NN is trained until its parameters converge. At this stage, ${\widehat{W}}^{*}$ is obtained and the approximate dynamic model of the AV can be obtained.

Remark 1. 

The process of training $\widehat{W}\left(k\right)$ involves first tuning the throttle and brake, that is, recording the impact of the throttle and brake on speed. Then, the collected data are adopted for training to obtain ${\widehat{W}}^{*}$. Therefore, the process of training the NN parameter is completed offline.

2.3. Linearization Model and Constraints of the AV

Note that (3) can be rewritten in the following linear form:

$$\begin{array}{cc}\hfill \xi (k+1)=& A\xi \left(k\right)+Bu\left(k\right)+Dd\left(k\right),\hfill \end{array}$$

(11)

where

$$\begin{array}{cc}\hfill A=& \left[\begin{array}{ccc}1& 0& Tcos\left(\overline{\theta }\right)\\ 0& 1& Tsin\left(\overline{\theta }\right)\\ 0& 0& 1+T{\widehat{W}}_1^{*}{a}_{\sigma }\end{array}\right],B=\left[\begin{array}{c}0\\ 0\\ T{\widehat{W}}_2^{*}{b}_{\sigma }\end{array}\right]=D,\hfill \\ \hfill d\left(k\right)=& \frac{1}{{\widehat{W}}_2^{*}{b}_{\sigma }}\left({\widehat{W}}^{*T}\sigma \left(k\right)-e(k+1)-{\widehat{W}}_1^{*}{a}_{\sigma }v\left(k\right)\right)-u\left(k\right),\hfill \end{array}$$

where ${\widehat{W}}_i^{*}$$(i=1,2,3)$ indicates the i-th element of ${\widehat{W}}^{*}$; and $a_{\sigma }$ and $b_{\sigma }$ represent the linearization coefficients of $\sigma \left(k\right)$ with $u_o\left(k\right)=0$. Assume that $\left|d\left(k\right)\right|\le \overline{d}$ with $\overline{d}>0$ being a known parameter.

Algorithm 1 Algorithm for approximating AVs’ dynamics using a NN

Brake data collection: 1. Divide $V_{max}$ into n equal parts for AVs with the maximum speed of $V_{max}$. Denote $a_i=i\frac{V_{max}}{m}$, where $(i=1,2,\cdots ,m)$ is the step size for each increase in velocity. Divide the maximum brake opening, 9, equally into n parts, and the size of each brake opening is $\frac{9}{n}$. Denote $b_j=j\frac{9}{n}$, where $(j=1,2,\cdots ,n)$ denotes the step size for each increase in the brake opening value. Set $i=1$ and $j=1$. Go to step 2. 2. Set the AV to travel at velocity $a_i$. Then, apply the AV brake with the brake opening of $b_j$. Go to step 3. 3. Collect the data of the velocity $v(a_i,b_j)$ and acceleration $a(a_i,b_j)$ from the start of braking until the AV comes to a complete stop. Go to step 4. 4. If $i+1\le m$, $i=i+1$. Go to step 2. Otherwise, go to step 5. 5. If $j+1\le \frac{9}{n}$, set $i=1$ and $j=j+1$. Go to step 2. Otherwise, go to step 6. Throttle data collection: 6. Divide $V_{max}$ into p equal parts for AVs with the maximum speed of $V_{max}$. Denote $c_i=i\frac{V_{max}}{p}$, $(i=0,1,2,\cdots ,p)$ as the step size for each increase in velocity. Divide the maximum throttle opening 1 equally into q parts, and the size of each brake opening is $\frac{1}{q}$. Denote $d_j=j\frac{1}{q}$, $(j=1,2,\cdots ,q)$ as the step size for each increase in the brake opening value. Set $i=0$ and $j=1$. Go to step 7. 7. Set the AV to travel at a fixed velocity $c_i$ and accelerate at a throttle opening of $d_j$. Go to step 8. 8. Collect the data of the velocity $v(c_i,d_j)$ and acceleration $a(c_i,d_j)$ from the beginning of acceleration until the AV reaches $V_{max}$. Go to 9. 9. If $i+1\le p$, $i=i+1$. Go to step 7. Otherwise, go to step 10. 10. If $j+1\le \frac{1}{q}$, set $i=1$ and $j=j+1$. Go to step 7. Otherwise, go to 11. Training and testing the NN: 11. Set the input of the NN as $b_j$ and $v(a_i,b_j)$, with the output being $a(a_i,b_j)$. Alternatively, the input of NN can be set as $c_i$ and $v(c_i,d_j)$, with the output being $a(c_i,d_j)$. Go to step 12. 12. Extract $70\%$ of these data as training data, and the remaining data as test data for backup. Go to step 13. 13. Define the NN structure, loss function, and activation function $\sigma (·)$. Using the training data obtained in step 12, the NN is trained until the network parameters converge. Go to step 13. 14. Using the NN obtained in step 12. If the loss meets the requirements, it indicates successful training. Terminate the algorithm and export the NN parameters ${\widehat{W}}^{*}$. If the loss does not meet the requirements, reset the training parameters and retrain. Go to step 13.

Considering the driving safety of the AV, its state and control input are subject to the following constraints:

$$\begin{array}{cc}\hfill \xi \left(k\right)\in & \mathcal{X}\triangleq \{\xi \left(k\right):{\tilde{h}}_{i,min,\xi }\le {\xi }_i\left(k\right)\le {\tilde{h}}_{i,max,\xi }\},\hfill \end{array}$$

(12)

$$\begin{array}{cc}\hfill u\left(k\right)\in & \mathcal{U}\triangleq \{u\left(k\right):{\tilde{h}}_{min,u}\le u\left(k\right)\le {\tilde{h}}_{max,u}\},\hfill \end{array}$$

(13)

where ${\xi }_i\left(k\right)$$(i=1,2,3)$ indicates the i-th element of $\xi \left(k\right)$; ${\tilde{h}}_{i,min,\xi }$ and ${\tilde{h}}_{i,max,\xi }$ represent the given minimum and maximum values of ${\xi }_i\left(k\right)$, respectively; and ${\tilde{h}}_{min,u}$ and ${\tilde{h}}_{max,u}$ denote the given minimum and maximum values of $u\left(k\right)$, respectively. Then, (12) and (13) can be rewritten as follows:

$$\begin{array}{cc}\hfill \xi \left(k\right)\in & \mathcal{X}\triangleq \{\xi \left(k\right):b_{\xi }^T\xi \left(k\right)\le h_{\xi }\},\hfill \end{array}$$

(14)

$$\begin{array}{cc}\hfill u\left(k\right)\in & \mathcal{U}\triangleq \{u\left(k\right):b_u^{T}u\left(k\right)\le h_u\},\hfill \end{array}$$

(15)

where $b_{\xi }$, $b_u$, $h_{\xi }$, and $h_u$ are matrices obtained from (12) and (13).

2.4. The Design of RHESO

In this subsection, a RHESO is designed to accurately estimate the unknown disturbance $d\left(k\right)$, as follows:

$$\begin{array}{cc}\hfill m(k+1)=& (I-L_{d}D)(m\left(k\right)+L_d\xi \left(k\right))-L_d(A\xi \left(k\right)+Bu\left(k\right))+b\left(k\right),\hfill \end{array}$$

(16a)

$$\begin{array}{cc}\hfill \widehat{d}\left(k\right)=& m\left(k\right)+L_d\xi \left(k\right),\hfill \end{array}$$

(16b)

where $m\left(k\right)\in \mathbb{R}$ is an intermediate variable; $\widehat{d}\left(k\right)\in \mathbb{R}$ is an estimation of $d\left(k\right)$; $L_d\in {\mathbb{R}}^{1\times 3}$ is an observer gain; and $b\left(k\right)\in \mathbb{R}$ is a decision variable to be designed. As a result, both $L_d$ and $b\left(k\right)$ need to be designed subsequently.

Furthermore, system (11) can be extended to the following system:

$$\begin{array}{cc}\hfill \xi (k+1)=& A\xi \left(k\right)+Bu\left(k\right)+Dd\left(k\right),\hfill \end{array}$$

(17a)

$$\begin{array}{cc}\hfill d(k+1)=& d\left(k\right)+h\left(k\right),\hfill \end{array}$$

(17b)

where $h\left(k\right)\in \mathbb{R}$ is an unknown bounded nonlinear term satisfying $\left|h\right(k\left)\right|\le H$, which is equivalent to $b_h^{T}h\left(k\right)\le h_h$ with $b_h=[1-1]$ and $h_h={\left[HH\right]}^T$.

Considering $e_d\left(k\right)=d\left(k\right)-\widehat{d}\left(k\right)$, we obtain the following:

$$\begin{array}{cc}\hfill e_d(k+1)=& (\mathbf{I}-L_{d}D)e_d\left(k\right)-b\left(k\right)+h\left(k\right).\hfill \end{array}$$

(18)

To restrict the transient performance of the RHESO (16a) and (16b), the following constraint of $e_d\left(k\right)$ is given as follows:

$$\begin{array}{c}\hfill e_d\left(k\right)\in \mathcal{E}\triangleq \{e_d\left(k\right):-{\tilde{h}}_{e_d}\le e_d\left(k\right)\le {\tilde{h}}_{e_d}\}\end{array}$$

(19)

where ${\tilde{h}}_{e_d}>0$ is a known parameter. Then, (19) can be rewritten as follows:

$$\begin{array}{c}\hfill e_d\left(k\right)\in \mathcal{E}=\{e_d\left(k\right):b_{e_d}^T{e}_d\left(k\right)\le h_{e_d}\}\end{array}$$

(20)

where $b_{e_d}=[1-1]$ and $h_{e_d}={\left[{\tilde{h}}_{e_d}{\tilde{h}}_{e_d}\right]}^T$.

Remark 2. 

Different from the traditional ESO, the RHESO proposed in (16a) and (16b) allows us to obtain the optimal estimation by optimizing the decision variable $b\left(k\right)$. Furthermore, it is possible to avoid significant overshoot caused by the estimation error $e_d\left(k\right)$ in terms of the constraint (19), thereby ensuring that system (11) does not experience performance degradation due to a large estimation error $e_d\left(k\right)$ before convergence.

2.5. Control Input and System Reconstruction

In this subsection, the control input is designed for system (11). Firstly, the control input of system (11) is designed as follows:

$$\begin{array}{c}\hfill u\left(k\right)=\left\{\begin{array}{cc}K\xi \left(k\right)+c\left(k\right)-V\widehat{d}\left(k\right),& \xi \left(k\right)\notin {\mathcal{X}}_T\\ K\xi \left(k\right)-V\widehat{d}\left(k\right),& \xi \left(k\right)\in {\mathcal{X}}_T\end{array}\right.\end{array}$$

(21)

where $V\triangleq {\left[01\right]}^T$; ${\mathcal{X}}_T$ is a terminal set to be designed; $K\in {\mathbb{R}}^{2\times 3}$ is the feedback gain to be designed; and $c\left(k\right)\in {\mathbb{R}}^2$ is the decision variable to be designed.

Secondly, based on (11), (18) and (21), the extended state system is obtained as follows:

$$\begin{array}{c}\hfill \check{\eta }(k+1)={\tilde{A}}_{\eta }\check{\eta }\left(k\right)+{\tilde{B}}_{\eta }{u}_{\eta }\left(k\right)+d_{\eta }\left(k\right)+\varpi \left(k\right),\end{array}$$

(22)

where

$$\begin{array}{cc}\hfill u_{\eta }\left(k\right)& =K_{\eta }\check{\eta }\left(k\right)+c_{\eta }\left(k\right)-{\left[\widehat{d}\left(k\right)0\right]}^T,\hfill \\ \hfill \check{\eta }\left(1\right)& \triangleq \left[\begin{array}{c}\xi \left(1\right)\\ 0\end{array}\right],\check{\eta }\left(k\right)\triangleq \left[\begin{array}{c}\xi \left(k\right)\\ e_d(k-1)\end{array}\right],{\tilde{A}}_{\eta }\triangleq \left[\begin{array}{cc}A& \mathbf{0}\\ \mathbf{0}& I\end{array}\right],\hfill \\ \hfill {\tilde{B}}_{\eta }& \triangleq \left[\begin{array}{cc}B& \mathbf{0}\\ \mathbf{0}& -I\end{array}\right],d_{\eta }\left(k\right)\triangleq \left[\begin{array}{c}Dd\left(k\right)\\ 0\end{array}\right],K_{\eta }\triangleq \left[\begin{array}{cc}K& \mathbf{0}\\ \mathbf{0}& L_{d}D\end{array}\right],\hfill \end{array}$$

(23)

$$\begin{array}{cc}\hfill \varpi \left(k\right)& \triangleq \left[\begin{array}{c}\mathbf{0}\\ h(k-1)\end{array}\right],k\in \mathbb{Z}[2,+\infty ),\hfill \end{array}$$

(24)

$$\begin{array}{cc}\hfill c_{\eta }\left(k\right)& =\left[\begin{array}{c}c\left(k\right)\\ b(k-1)\end{array}\right],k\in \mathbb{Z}[2,+\infty ).\hfill \end{array}$$

(25)

Then, we obtain the following:

$$\begin{array}{cc}\hfill F\check{\eta }(k+1)=& ({\tilde{A}}_{\eta }+{\tilde{B}}_{\eta }{K}_{\eta })\check{\eta }\left(k\right)+{\tilde{B}}_{\eta }{c}_{\eta }\left(k\right)+\varpi \left(k\right)={\tilde{\Phi }}_{\eta }\check{\eta }\left(k\right)+{\tilde{B}}_{\eta }{c}_{\eta }\left(k\right)+\varpi \left(k\right)\hfill \end{array}$$

(26)

where

$$\begin{array}{cc}\hfill F& \triangleq \left[\begin{array}{cc}I& -D\\ \mathbf{0}& I\end{array}\right],{\tilde{\Phi }}_{\eta }\triangleq \left[\begin{array}{cc}\Phi & \mathbf{0}\\ \mathbf{0}& {\Phi }_{e_d}\end{array}\right],\Phi \triangleq A+BK,{\Phi }_{e_d}\triangleq I-L_{d}D.\hfill \end{array}$$

Since F is a full rank matrix, (26) is converted as follows:

$$\begin{array}{cc}\hfill \check{\eta }(k+1)=& (F^{-1}{\tilde{A}}_{\eta }+F^{-1}{\tilde{B}}_{\eta }{K}_{\eta })\check{\eta }\left(k\right)+F^{-1}{\tilde{B}}_{\eta }{c}_{\eta }\left(k\right)+F^{-1}\varpi \left(k\right)\hfill \\ \hfill =& (A_{\eta }+B_{\eta }{K}_{\eta })\check{\eta }\left(k\right)+B_{\eta }{c}_{\eta }\left(k\right)+F^{-1}\varpi \left(k\right)\hfill \\ \hfill =& {\Phi }_{\eta }\check{\eta }\left(k\right)+B_{\eta }{c}_{\eta }\left(k\right)+F^{-1}\varpi \left(k\right),\hfill \end{array}$$

(27)

where ${\Phi }_{\eta }=A_{\eta }+B_{\eta }{K}_{\eta }$, $A_{\eta }=F^{-1}{\tilde{A}}_{\eta }$ and $B_{\eta }=F^{-1}{\tilde{B}}_{\eta }$. According to (14), (15) and (19), the constraints of $\check{\eta }\left(k\right)$ and $c_{\eta }\left(k\right)$ are obtained as follows:

$$\begin{array}{cc}\hfill \check{\eta }\left(k\right)& \triangleq \{\check{\eta }\left(k\right):b_{\eta }^T\check{\eta }\left(k\right)\le h_{\eta }\},\hfill \\ \hfill c_{\eta }\left(k\right)& \triangleq \{c_{\eta }\left(k\right):{\tilde{b}}_u^T\check{\eta }\left(k\right)+b_{c_{\eta }}^T{c}_{\eta }\left(k\right)\le h_{c_{\eta }}+h_d\left(k\right)\}\hfill \end{array}$$

where

$$\begin{array}{cc}\hfill b_{\eta }& \triangleq \left[\begin{array}{cc}{b}_{\xi }& \mathbf{0}\\ \mathbf{0}& b_{e_d}\end{array}\right],{\tilde{b}}_u\triangleq \left[\begin{array}{cc}{K}^T{b}_u& \mathbf{0}\\ \mathbf{0}& \mathbf{0}\end{array}\right],b_{c_{\eta }}\triangleq \left[\begin{array}{cc}{b}_u& \mathbf{0}\\ \mathbf{0}& \mathbf{0}\end{array}\right],\hfill \\ \hfill h_{\eta }& \triangleq \left[\begin{array}{c}{h}_{\xi }\\ h_{e_d}\end{array}\right],h_{c_{\eta }}\triangleq \left[\begin{array}{c}{h}_u\\ \mathbf{0}\end{array}\right],h_d\left(k\right)\triangleq \left[\begin{array}{c}{b}_u^T\widehat{d}\left(k\right)\\ \mathbf{0}\end{array}\right].\hfill \end{array}$$

Based on $\left|h\right(k\left)\right|\le H$, the constraint of $\varpi \left(k\right)$ is obtained as follows:

$$\begin{array}{c}\hfill \varpi \left(k\right)\in \mathcal{W}\triangleq \{\varpi \left(k\right):b_{\varpi }^T\varpi \left(k\right)\le h_{\varpi }\},\end{array}$$

where

$$\begin{array}{c}\hfill b_{\varpi }\triangleq \left[\begin{array}{cc}\mathbf{0}& \mathbf{0}\\ \mathbf{0}& b_h\end{array}\right],h_{\varpi }\triangleq \left[\begin{array}{c}\mathbf{0}\\ {\overline{h}}_h\end{array}\right],{\overline{h}}_h\triangleq \left[\begin{array}{c}H\\ H\end{array}\right].\end{array}$$

The nominal system of (27) is denoted as follows:

$$\begin{array}{c}\hfill \eta (k+1)={\Phi }_{\eta }\eta \left(k\right)+B_{\eta }{c}_{\eta }\left(k\right),\end{array}$$

where $\eta \left(k\right)$ is the nominal state of $\check{\eta }\left(k\right)$.

Denoting $e_{\eta }\left(k\right)\triangleq \check{\eta }\left(k\right)-\eta \left(k\right)$, we have the following:

$$\begin{array}{c}\hfill e_{\eta }(k+1)={\Phi }_{\eta }{e}_{\eta }\left(k\right)+F^{-1}\varpi \left(k\right).\end{array}$$

(28)

Then, constraints of $\eta \left(k\right)$ and $c_{\eta }\left(k\right)$ are obtained as follows:

$$\begin{array}{cc}\hfill \eta \left(k\right)& \in \mathcal{M}\triangleq \{\eta \left(k\right):b_{\eta }^T\eta \left(k\right)\le h_{\eta }-b_{\eta }^T{e}_{\eta }\left(k\right)\},\hfill \\ \hfill c_{\eta }\left(k\right)& \in \mathcal{N}\triangleq \{c_{\eta }\left(k\right):{\tilde{b}}_u^T\eta \left(k\right)+b_{c_{\eta }}^T{c}_{\eta }\left(k\right)\le h_{c_{\eta }}-({\tilde{b}}_u^T{e}_{\eta }\left(k\right)-h_d\left(k\right))\}.\hfill \end{array}$$

2.6. The MPC Scheme

To achieve smooth braking at preset points in AVs under constraints, an MPC scheme is proposed in this subsection. Firstly, the cost function of system (27) is denoted as follows:

$$\begin{array}{cc}\hfill & J(\widehat{\mathbf{ȷ}}\left(k\right),{\widehat{\mathbf{c}}}_{\eta }\left(k\right))\hfill \\ \hfill \triangleq & \sum _{s=0}^{N-1}\left(\parallel \widehat{\eta }{(k+s|k)\parallel }_Q^2+{\parallel K_{\eta }\widehat{\eta }(k+s|k)+{\widehat{c}}_{\eta }(k+s|k)\parallel }_R^2\right)+{\parallel \widehat{\eta }(k+N|k)\parallel }_P^2,\hfill \end{array}$$

where $\widehat{\eta }(k+s|k)$ is the predictive state with $\widehat{\eta }\left(k\right|k)=\check{\eta }\left(k\right)$; ${\widehat{c}}_{\eta }(k+s|k)$ is the predictive decision variable; $\widehat{\mathbf{ȷ}}\left(k\right)\triangleq {[{\widehat{\eta }}^T\left(k\right|k){\widehat{\eta }}^T(k+1|k)\cdots {\widehat{\eta }}^T(k+N|k)]}^T$ is the predictive state sequence; ${\widehat{\mathbf{c}}}_{\eta }\left(k\right)\triangleq {[{\widehat{c}}_{\eta }^T\left(k\right|k){\widehat{c}}_{\eta }^T(k+1|k)\cdots {\widehat{c}}_{\eta }^T(k+N-1|k)]}^T$ is the predictive decision variable sequence with ${\widehat{c}}_{\eta }\left(k\right|k)\triangleq c_{\eta }\left(k\right)$; $Q>0$ and $R>0$ are the weighting matrices with appropriate dimensions; $P>0$ is the terminal penalty matrix to be designed; and N is the predictive horizon. Then, the optimization problem of MPC is denoted as follows:

Prob 1:

$$\begin{array}{c}\hfill ({\widehat{\mathbf{ȷ}}}^{*}\left(k\right),{\widehat{\mathbf{c}}}_{\eta }^{*}\left(k\right))\triangleq argminJ(\widehat{\mathbf{ȷ}}\left(k\right),{\widehat{\mathbf{c}}}_{\eta }\left(k\right))\end{array}$$

subject to

$$\begin{array}{cc}\hfill & \widehat{\eta }(k+s+1|k)={\Phi }_{\eta }\widehat{\eta }(k+s|k)+B_{\eta }{\widehat{c}}_{\eta }(k+s|k),\hfill \end{array}$$

(29)

$$\begin{array}{cc}\hfill & b_{\eta }^T\widehat{\eta }(k+s|k)\le h_{\eta }-\varphi (k+s|k),s\in \mathbb{Z}[0,N-1],\hfill \end{array}$$

(30)

$$\begin{array}{cc}\hfill & {\tilde{b}}_u^T\widehat{\eta }(k+s|k)+b_{c_{\eta }}^T{\widehat{c}}_{\eta }(k+s|k)\le h_{c_{\eta }}-|{\tilde{h}}_{e_d}+\overline{d}|{\left[b_{u}0\right]}^T-\chi (k+s|k),s\in \mathbb{Z}[0,N-1],\hfill \end{array}$$

(31)

$$\begin{array}{cc}\hfill & \widehat{\eta }(k+N|k)\in {\mathcal{X}}_T(k+N|k),\hfill \end{array}$$

(32)

where

$$\begin{array}{cc}\hfill & \varphi (k+s|k)\triangleq \sum _{j=0}^{s-1}\underset{\varpi (k+s-1-j|k)\in \mathcal{W}}{max}{b}_{\eta }^T{\Phi }_{\eta }^j{F}^{-1}\varpi (k+s-1-j|k),\varphi \left(k\right|k)\triangleq \mathbf{0},\hfill \\ \hfill & \chi (k+s|k)\triangleq \sum _{j=0}^{s-1}\underset{\varpi (k+s-1-j|k)\in \mathcal{W}}{max}{\tilde{b}}_u^T{\Phi }_{\eta }^j{F}^{-1}\varpi (k+s-1-j|k),\chi \left(k\right|k)\triangleq \mathbf{0},\hfill \\ \hfill & e_{\eta }(k+s+1|k)={\Phi }_{\eta }{e}_{\eta }(k+s|k)+F^{-1}\varpi (k+s|k),e_{\eta }\left(k\right|k)=\mathbf{0},\hfill \\ \hfill & {\mathcal{X}}_T(k+N|k)\triangleq \widehat{\mathcal{X}}(k+N|k)\cap \widehat{\mathcal{C}}(k+N|k),\hfill \\ \hfill & \widehat{\mathcal{X}}(k+N|k)\triangleq \{\widehat{\eta }(k+N|k):b_{\eta }^T\widehat{\eta }(k+N|k)\le h_{\eta }-\varphi (k+N|k)\},\hfill \\ \hfill & \widehat{\mathcal{C}}(k+N|k)\triangleq \{\widehat{\eta }(k+N|k):{\tilde{b}}_u^T\widehat{\eta }(k+N|k)\le h_{c_{\eta }}-|{\tilde{h}}_{e_d}+\overline{d}|{\left[b_{u}0\right]}^T-\chi (k+N|k)\},\hfill \\ \hfill & {\widehat{\mathbf{ȷ}}}^{*}\left(k\right)\triangleq {[{\widehat{\eta }}^{*T}\left(k\right|k){\widehat{\eta }}^{*T}(k+1|k)\cdots {\widehat{\eta }}^{*T}(k+N|k)]}^T,\hfill \\ \hfill & {\widehat{\mathbf{c}}}_{\eta }^{*}\left(k\right)\triangleq {[{\widehat{c}}^{*T}\left(k\right|k){\widehat{c}}^{*T}(k+1|k)\cdots {\widehat{c}}^{*T}(k+N-1|k)]}^T.\hfill \end{array}$$

Remark 3. 

The polynomial complexity of the proposed RHESO-based MPC consists of three parts: (1) the polynomial complexity of the RHESO; (2) the polynomial complexity of Prob 1; and (3) the polynomial complexity of the input signal. The polynomial complexity of the RHESO is obtained as follows:

$$\begin{array}{cc}\hfill & \mathcal{O}\left((I-L_{d}D)(m\left(k\right)+L_d\xi \left(k\right))-L_d(A\xi \left(k\right)+Bu\left(k\right))+b\left(k\right)\right)+\mathcal{O}\left(m\left(k\right)+L_d\xi \left(k\right)\right)\hfill \\ \hfill =& \mathcal{O}\left(24\right)\hfill \end{array}$$

The polynomial complexity of Prob 1 is obtained as follows:

$$\begin{array}{cc}\hfill & \mathcal{O}\left(\sum _{s=0}^{N-1}\left(\parallel \widehat{\eta }{(k+s|k)\parallel }_Q^2+{\parallel K_{\eta }\widehat{\eta }(k+s|k)+{\widehat{c}}_{\eta }(k+s|k)\parallel }_R^2\right)+{\parallel \widehat{\eta }(k+N|k)\parallel }_P^2\right)\hfill \\ \hfill =& \mathcal{O}\left(52N+31\right)\hfill \end{array}$$

The polynomial complexity of the input signal is obtained as follows:

$$\begin{array}{cc}\hfill & \mathcal{O}\left(K\xi \left(k\right)+c\left(k\right)-V\widehat{d}\left(k\right)\right)=\mathcal{O}\left(13\right)\hfill \end{array}$$

Then the polynomial complexity of the proposed RHESO-based MPC is obtained as follows:

$$\begin{array}{c}\hfill \mathcal{O}\left(24\right)+\mathcal{O}\left(52N+31\right)+\mathcal{O}\left(13\right)=\mathcal{O}\left(52N+31\right)\end{array}$$

Definition 1 

([38]). If for all ${\overline{\Xi }}_1>0$, there exists a constant ${\overline{\Xi }}_2\left({\overline{\Xi }}_1\right)$ such that we have the following:

$$\begin{array}{c}\hfill \parallel \check{\eta }\left(0\right)\parallel \le {\overline{\Xi }}_1\Rightarrow \parallel \check{\eta }\left(k\right)\parallel \le {\overline{\Xi }}_2\left({\overline{\Xi }}_1\right)\end{array}$$

holds for all $\check{\eta }\left(0\right)\in {\mathbb{R}}^n$ and $k\in {\mathbb{Z}}_{\ge 0}$, then system (27) is said to be uniformly bounded.

3. Results

In this part, a lemma is presented firstly to demonstrate the feasibility of Prob 1, and then a theorem is presented to verify the uniformly bounded stability of the closed-loop system.

Lemma 1. 

If a feasible solution exists for Prob 1 at time instant k, then it also has at least one feasible solution at time instant $k+1$.

Proof. 

Denote a candidate decision variable sequence ${\tilde{\mathbf{c}}}_{\eta }(k+1)$ as

$$\begin{array}{cc}\hfill & {\tilde{\mathbf{c}}}_{\eta }(k+1)\hfill \\ \hfill \triangleq & \left[\begin{array}{ccc}{\tilde{c}}_{\eta }(k+1|k+1)& \cdots & {\tilde{c}}_{\eta }(k+N|k+1)\end{array}\right]=\left[\begin{array}{cccc}{\widehat{c}}_{\eta }^{*}(k+1|k)& \cdots & {\widehat{c}}_{\eta }^{*}(k+N-1|k)& \mathbf{0}\end{array}\right],\hfill \end{array}$$

where ${\tilde{c}}_{\eta }(k+1+s|k+1)$ is a candidate decision variable for $s\in \mathbb{Z}[0,N-1]$. Then, there exists

$$\begin{array}{cc}\hfill \widehat{\eta }(k+1|k+1)& ={\Phi }_{\eta }\widehat{\eta }\left(k\right|k)+B_{\eta }{\widehat{c}}_{\eta }^{*}\left(k\right|k)+F^{-1}\varpi \left(k\right).\hfill \end{array}$$

(33)

Considering the predictive dynamics of $\check{\eta }(k+1)$ using ${\tilde{\mathbf{c}}}_{\eta }(k+1)$ at time instant $k+1$, we have the following:

$$\begin{array}{cc}\hfill \tilde{\eta }(k+1+s+1|k+1)=& {\Phi }_{\eta }\tilde{\eta }(k+1+s|k+1)+B_{\eta }{\tilde{c}}_{\eta }(k+1+s|k+1)\hfill \\ \hfill & +F^{-1}\varpi (k+1+s|k+1),\hfill \end{array}$$

(34)

where $\tilde{\eta }(k+1+s+1|k+1)$ is the feasible state for $s\in \mathbb{Z}[0,N]$.

Combining (29), (33) and (34), we have the following:

$$\begin{array}{c}\hfill \tilde{\eta }(k+s|k+1)={\widehat{\eta }}^{*}(k+s|k)+{\Phi }_{\eta }^{s-1}{F}^{-1}\varpi \left(k\right).\end{array}$$

(35)

Denote

$$\begin{array}{cc}\hfill & \tilde{\mathbf{ȷ}}(k+1)=\left[\begin{array}{ccc}\tilde{\eta }(k+1|k+1)& \cdots & \tilde{\eta }(k+1+N|k+1)\end{array}\right]\hfill \end{array}$$

as the feasible predictive state sequence. To guarantee that Prob 1 is feasible at time instant $k+1$, the following conditions should be satisfied:

$$\begin{array}{cc}\hfill & b_{\eta }^T\tilde{\eta }(k+1+s|k+1)\le h_{\eta }-\varphi (k+1+s|k+1),s\in \mathbb{Z}[0,N]\hfill \end{array}$$

(36)

$$\begin{array}{cc}\hfill {\tilde{b}}_u^T\tilde{\eta }(k+1+s|k+1)& + b_{c_{\eta }}^T{\tilde{c}}_{\eta }(k+1+s|k+1)\le h_{c_{\eta }}-|{\tilde{h}}_{e_d}+\overline{d}|{\left[b_{u}0\right]}^T\hfill \\ \hfill & -\chi (k+1+s|k+1),s\in \mathbb{Z}[0,N-1]\hfill \end{array}$$

(37)

$$\begin{array}{cc}\hfill & \tilde{\eta }(k+1+N|k+1)\in {\mathcal{X}}_T(k+1+N|k+1)\hfill \end{array}$$

(38)

In the following, conditions (36)–(38) are proved one by one. Firstly, the proof of condition (36) is given. Assume that Prob 1 is feasible at time instant k. Then, the constraint (30) holds, which indicates the following:

$$\begin{array}{c}\hfill b_{\eta }^T{\widehat{\eta }}^{*}(k+1+s|k)\le h_{\eta }-\varphi (k+1+s|k)\end{array}$$

(39)

for $s\in \mathbb{Z}[0,N-2]$. Combining (35) and (39), the following exists:

$$\begin{array}{cc}\hfill & b_{\eta }^T(\tilde{\eta }(k+1+s|k+1)-{\Phi }_{\eta }^s{F}^{-1}\varpi \left(k\right))\le h_{\eta }-\varphi (k+1+s|k),\hfill \end{array}$$

which yields

$$\begin{array}{cc}\hfill & b_{\eta }^T\tilde{\eta }(k+1+s|k+1)\hfill \\ \hfill & \le h_{\eta }-\varphi (k+1+s|k)+b_{\eta }^T{\Phi }_{\eta }^s{F}^{-1}\varpi \left(k\right)\hfill \\ \hfill & =h_{\eta }-\sum _{j=0}^{s-1}\underset{\varpi (k+s-1-j|k)\in \mathcal{W}}{max}{b}_{\eta }^T{\Phi }_{\eta }^j{F}^{-1}\varpi (k+s-1-j|k)\hfill \\ \hfill & +b_{\eta }^T{\Phi }_{\eta }^s{F}^{-1}\left(\varpi \left(k\right)-\underset{\varpi (k+s-j|k)\in \mathcal{W}}{max}\varpi (k+s-j|k)\right)\hfill \\ \hfill & \le h_{\eta }-\sum _{j=0}^{s-1}\underset{\varpi (k+s-1-j|k)\in \mathcal{W}}{max}{b}_{\eta }^T{\Phi }_{\eta }^j{F}^{-1}\varpi (k+s-1-j|k)=h_{\eta }-\varphi (k+1+s|k+1)\hfill \end{array}$$

for $s\in \mathbb{Z}[0,N]$. Therefore, condition (36) is satisfied for $s\in \mathbb{Z}[0,N]$.

Secondly, the proof of condition (37) is given. By using the candidate variable sequence ${\tilde{\mathbf{c}}}_{\eta }(k+1)$, we have that ${\tilde{c}}_{\eta }(k+1+s|k+1)={\widehat{c}}_{\eta }(k+1+s|k)$. Based on (31) and (35), we have the following:

$$\begin{array}{cc}\hfill & {\tilde{b}}_u^T(\tilde{\eta }(k+1+s|k+1)-{\Phi }_{\eta }^s{F}^{-1}\varpi \left(k\right))+b_{c_{\eta }}^T{\widehat{c}}_{\eta }(k+1+s|k)\hfill \\ \hfill \le & h_{c_{\eta }}-|h_{e_d}+\overline{d}|{\left[b_{u}0\right]}^T-\chi (k+1+s|k),\hfill \end{array}$$

which indicates the following:

$$\begin{array}{cc}\hfill & {\tilde{b}}_u^T\tilde{\eta }(k+1+s|k+1)+b_{c_{\eta }}^T{\tilde{c}}_{\eta }(k+1+s|k+1)\hfill \\ \hfill & \le h_{c_{\eta }}-|h_{e_d}+\overline{d}|{\left[b_{u}0\right]}^T-\chi (k+1+s|k)+{\tilde{b}}_u^T{\Phi }_{\eta }^s{F}^{-1}\varpi \left(k\right)\hfill \\ \hfill & =h_{c_{\eta }}-|{\tilde{h}}_{e_d}+\overline{d}|{\left[b_{u}0\right]}^T-\sum _{j=0}^{s-1}\underset{\varpi (k+s-1-j|k)\in \mathcal{W}}{max}{\tilde{b}}_u^T{\Phi }_{\eta }^j{F}^{-1}\varpi (k+s-1-j|k)\hfill \\ \hfill & +{\tilde{b}}_u^T{\Phi }_{\eta }^s{F}^{-1}\left(\varpi \left(k\right)-\underset{\varpi (k+s-j|k)\in \mathcal{W}}{max}\varpi (k+s-j|k)\right)\hfill \\ \hfill & \le h_{c_{\eta }}-|{\tilde{h}}_{e_d}+\overline{d}|{\left[b_{u}0\right]}^T-\sum _{j=0}^{s-1}\underset{\varpi (k+s-1-j|k)\in \mathcal{W}}{max}{\tilde{b}}_u^T{\Phi }_{\eta }^j{F}^{-1}\varpi (k+s-1-j|k)\hfill \\ \hfill & =h_{c_{\eta }}-|{\tilde{h}}_{e_d}+\overline{d}|{\left[b_{u}0\right]}^T-\chi (k+1+s|k+1)\hfill \end{array}$$

for $s\in \mathbb{Z}[0,N-2]$. For $s=N-1$, there exists ${\tilde{c}}_{\eta }(k+N|k+1)=\mathbf{0}$. Based on (32), (35) and ${\tilde{c}}_{\eta }(k+N|k+1)=\mathbf{0}$, it is easy to obtain the following:

$$\begin{array}{cc}\hfill & {\tilde{b}}_u^T(\tilde{\eta }(k+N|k+1)-{\Phi }_{\eta }^{N-1}{F}^{-1}\varpi \left(k\right))\le h_{c_{\eta }}-|{\tilde{h}}_{e_d}+\overline{d}|{\left[b_{u}0\right]}^T-\chi (k+N|k),\hfill \end{array}$$

which means that we have the following:

$$\begin{array}{cc}\hfill & {\tilde{b}}_u^T\tilde{\eta }(k+N|k+1)\hfill \\ \hfill \le & h_{c_{\eta }}-|{\tilde{h}}_{e_d}+\overline{d}|{\left[b_{u}0\right]}^T-\chi (k+N|k)+{\tilde{b}}_u^T{\Phi }_{\eta }^{N-1}{F}^{-1}\varpi \left(k\right)\hfill \\ \hfill =& h_{c_{\eta }}-|{\tilde{h}}_{e_d}+\overline{d}|{\left[b_{u}0\right]}^T-\sum _{j=0}^{N-2}\underset{\varpi (k+N-2-j|k)\in \mathcal{W}}{max}{\tilde{b}}_u^T{\Phi }_{\eta }^j{F}^{-1}\varpi (k+N-2-j|k)\hfill \\ \hfill & -{\tilde{b}}_u^T{\Phi }_{\eta }^{N-1}{F}^{-1}\underset{\varpi (k+N-1-j|k)\in \mathcal{W}}{max}\varpi (k+N-1-j|k)+{\tilde{b}}_u^T{\Phi }_{\eta }^{N-1}{F}^{-1}\varpi \left(k\right)\hfill \\ \hfill \le & h_{c_{\eta }}-|{\tilde{h}}_{e_d}+\overline{d}|{\left[b_{u}0\right]}^T-\sum _{j=0}^{N-2}\underset{\varpi (k+N-2-j|k)\in \mathcal{W}}{max}{\tilde{b}}_u^T{\Phi }_{\eta }^j{F}^{-1}\varpi (k+N-2-j|k)\hfill \\ \hfill =& h_{c_{\eta }}-|{\tilde{h}}_{e_d}+\overline{d}|{\left[b_{u}0\right]}^T-\chi (k+N|k+1).\hfill \end{array}$$

Therefore, condition (37) is satisfied for $s\in \mathbb{Z}[0,N-1]$.

Thirdly, the proof of condition (38) is proved, as ${\mathcal{X}}_T(k+1+N|k+1)\triangleq \widehat{\mathcal{X}}(k+1+N|k)\cap \widehat{\mathcal{C}}(k+1+N|k+1)$. In terms of (36), it follows that $\tilde{\eta }(k+1+N|k+1)\in \widehat{\mathcal{X}}(k+1+N|k)$. Furthermore, based on (32), it is easy to obtain that $\tilde{\eta }(k+1+N|k+1)\in \widehat{\mathcal{C}}(k+1+N|k)$. Therefore, $\tilde{\eta }(k+1+N|k+1)\in {\mathcal{X}}_T(k+1+N|k+1)$ is satisfied, such that condition (38) holds.

As a result, if Prob 1 is feasible at time instant k, then Prob 1 is also feasible at time instant $k+1$. This completes the proof. □

Theorem 2. 

If $\Xi >0$ and ${{\rm Y}}_{\eta }$ exist, which make the linear matrix inequality (LMI),

$$\begin{array}{cc}\hfill & \left[\begin{array}{cccc}-\Xi & \Xi A_{\eta }^T+{{\rm Y}}_{\eta }{B}_{\eta }^T& {{\rm Y}}_{\eta }& \Xi \\ A_{\eta }{\Xi }^T+B_{\eta }{{\rm Y}}_{\eta }^T& -\Xi & \mathbf{0}& \mathbf{0}\\ {{\rm Y}}_{\eta }^T& \mathbf{0}& -R^{-1}& \mathbf{0}\\ {\Xi }^T& \mathbf{0}& \mathbf{0}& -Q^{-1}\end{array}\right]<0\hfill \end{array}$$

(40)

then system (22) is uniformly bounded and $K_{\eta }={{\rm Y}}_{\eta }^T{\Xi }^{-1}$ with ${\Xi }^{-1}=P$. Furthermore, the state of system (22) converges to the set $\mathbb{S}\triangleq \{\check{\eta }\left(k\right):\parallel \check{\eta }\left(k\right){\parallel }_Q^2\le \overline{\sigma }\}$, where

$$\begin{array}{cc}\hfill \overline{\sigma }\triangleq & \sum _{s=1}^{N-1}\left\{{\overline{\lambda }}^2\left(Q\right){\Omega }^2\left(s\right)+2{\overline{\lambda }}^2\left(Q\right)\Omega \left(s\right)F\left(s\right)+{\overline{\lambda }}^2\left(R\right){\Xi }^2\left(s\right)+2{\overline{\lambda }}^2\left(R\right)\Xi \left(s\right)\Psi \left(s\right)\right\}\hfill \\ \hfill & +{\overline{\lambda }}^2\left(P\right){\Omega }^2\left(N\right)+\overline{\lambda }\left(P\right)\Omega \left(N\right)\pi \left(N\right),\hfill \end{array}$$

with

$$\begin{array}{cc}\hfill & \Omega \left(s\right)\triangleq \underset{\varpi \left(k\right)\in \mathcal{W}}{max}\parallel {\Phi }_{\eta }^{s-1}{F}^{-1}\varpi \left(k\right)\parallel ,F\left(s\right)\triangleq \parallel h_{\eta }-\varphi (k+s|k)\parallel ,\hfill \\ \hfill & \Xi \left(s\right)\triangleq \underset{\varpi \left(k\right)\in \mathcal{W}}{max}\parallel K_{\eta }{\Phi }_{\eta }^{s-1}{F}^{-1}\varpi \left(k\right)\parallel ,\Psi \left(s\right)\triangleq \parallel h_{c_{\eta }}-|h_{e_d}+\overline{d}|{\left[b_{u}0\right]}^T-\chi (k+s|k)\parallel .\hfill \end{array}$$

Proof. 

Considering the optimal solution of Prob 1 is obtained at time instant k, it can be concluded that Prob 1 has at least one feasible solution at time instant $k+1$. Denote $J(\tilde{\mathbf{ȷ}}(k+1),{\tilde{\mathbf{c}}}_{\eta }(k+1))$ as the feasible cost of Prob 1 at time instant $k+1$ with $\tilde{\mathbf{ȷ}}(k+1)$ and ${\tilde{\mathbf{c}}}_{\eta }(k+1)$ given in Lemma 1.

The optimal cost of Prob 1 at time instant k is denoted as $J({\widehat{\mathbf{ȷ}}}^{*}\left(k\right),{\widehat{\mathbf{c}}}_{\eta }^{*}\left(k\right))$. Choose $\Delta J\left(k\right)\triangleq J(\tilde{\mathbf{ȷ}}(k+1),{\tilde{\mathbf{c}}}_{\eta }(k+1))-J({\widehat{\mathbf{ȷ}}}^{*}\left(k\right),{\widehat{\mathbf{c}}}_{\eta }^{*}\left(k\right))$. Split $\Delta J\left(k\right)={\Delta }_1+{\Delta }_2+{\Delta }_3$ with

$$\begin{array}{cc}\hfill {\Delta }_1=& \sum _{s=1}^{N-1}\left\{\parallel \tilde{\eta }{(k+s|k+1)\parallel }_Q^2-\parallel {\widehat{\eta }}^{*}{(k+s|k)\parallel }_Q^2+{\parallel K_{\eta }\tilde{\eta }(k+s|k+1)+{\tilde{c}}_{\eta }(k+s|k+1)\parallel }_R^2\right.\hfill \\ \hfill & \left.-\parallel K_{\eta }{\widehat{\eta }}^{*}(k+s|k)+{\widehat{c}}_{\eta }^{*}{(k+s|k)\parallel }_R^2\right\},\hfill \\ \hfill {\Delta }_2=& \parallel \tilde{\eta }(k+T_s|k+1){\parallel }_Q^2+\parallel K_{\eta }\tilde{\eta }(k+T_s|k+1){\parallel }_R^2+{\parallel \tilde{\eta }(k+1+T_s|k+1)\parallel }_P^2\hfill \\ \hfill & -\parallel {\widehat{\eta }}^{*}(k+T_s|k){\parallel }_P^2,\hfill \\ \hfill {\Delta }_3=& -\parallel {\widehat{\eta }}^{*}{\left(k\right|k)\parallel }_Q^2-{\parallel K_{\eta }{\widehat{\eta }}^{*}\left(k\right|k)+{\widehat{c}}_{\eta }^{*}\left(k\right|k)\parallel }_Q^2.\hfill \end{array}$$

For ${\Delta }_1$, we have the following:

$$\begin{array}{cc}\hfill {\Delta }_1=& \sum _{s=1}^{N-1}\left\{\parallel \tilde{\eta }{(k+s|k+1)\parallel }_Q^2-\parallel {\widehat{\eta }}^{*}(k+s|k){\parallel }_Q^2+\parallel K_{\eta }({\widehat{\eta }}^{*}(k+s|k)+{\Phi }_{\eta }^{s-1}{F}^{-1}\varpi \left(k\right))\right.\hfill \\ \hfill & \left.+{\tilde{c}}_{\eta }{(k+s|k+1)\parallel }_R^2-{\parallel K_{\eta }{\widehat{\eta }}^{*}(k+s|k)+{\widehat{c}}_{\eta }^{*}(k+s|k)\parallel }_R^2\right\}\hfill \\ \hfill =& \sum _{s=1}^{N-1}\left\{\parallel \tilde{\eta }{(k+s|k+1)\parallel }_Q^2-\parallel {\widehat{\eta }}^{*}(k+s|k){\parallel }_Q^2+\parallel K_{\eta }{\widehat{\eta }}^{*}(k+s|k)+{\widehat{c}}_{\eta }^{*}(k+s|k)\right.\hfill \\ \hfill & \left.+K_{\eta }{\Phi }_{\eta }^{s-1}{F}^{-1}{\varpi \left(k\right)\parallel }_R^2-{\parallel K_{\eta }{\widehat{\eta }}^{*}(k+s|k)+{\widehat{c}}_{\eta }^{*}(k+s|k)\parallel }_R^2\right\}\hfill \\ \hfill \le & \sum _{s=1}^{N-1}\left\{\left(\parallel \tilde{\eta }{(k+s|k+1)\parallel }_Q-{\parallel {\widehat{\eta }}^{*}(k+s|k)\parallel }_Q\right)\left(\parallel \tilde{\eta }{(k+s|k+1)\parallel }_Q+{\parallel {\widehat{\eta }}^{*}(k+s|k)\parallel }_Q\right)\right.\hfill \\ \hfill & +\left(\parallel K_{\eta }{\widehat{\eta }}^{*}(k+s|k)+{\widehat{c}}_{\eta }^{*}(k+s|k)\right.+K_{\eta }{\Phi }_{\eta }^{s-1}{F}^{-1}{\varpi \left(k\right)\parallel }_R\hfill \\ \hfill & \left.-\parallel K_{\eta }{\widehat{\eta }}^{*}(k+s|k)+{\widehat{c}}_{\eta }^{*}{(k+s|k)\parallel }_R\right)\times \left(\parallel K_{\eta }{\widehat{\eta }}^{*}(k+s|k)+{\widehat{c}}_{\eta }^{*}(k+s|k)\right.\hfill \\ \hfill & \left.+K_{\eta }{\Phi }_{\eta }^{s-1}{F}^{-1}{\varpi \left(k\right)\parallel }_R\left.+\parallel K_{\eta }{\widehat{\eta }}^{*}(k+s|k)+{\widehat{c}}_{\eta }^{*}{(k+s|k)\parallel }_R\right)\right\}\hfill \\ \hfill \le & \sum _{s=1}^{N-1}\left\{\overline{\lambda }\left(Q\right)\parallel {\Phi }_{\eta }^{s-1}{F}^{-1}\varpi \left(k\right)\parallel \times \left(2\overline{\lambda }\left(Q\right)\parallel {\widehat{\eta }}^{*}(k+s|k)\parallel +\overline{\lambda }\left(Q\right)\parallel {\Phi }_{\eta }^{s-1}{F}^{-1}\varpi \left(k\right)\parallel \right)\right.\hfill \\ \hfill & +\overline{\lambda }\left(R\right)\parallel K_{\eta }{\Phi }_{\eta }^{s-1}{F}^{-1}\varpi \left(k\right)\parallel \left(2\overline{\lambda }\left(R\right)\parallel K_{\eta }{\widehat{\eta }}^{*}(k+s|k)+{\widehat{c}}_{\eta }^{*}(k+s|k)\parallel \right.\hfill \\ \hfill & \left.\left.+\overline{\lambda }\left(R\right)\parallel K_{\eta }{\Phi }_{\eta }^{s-1}{F}^{-1}\varpi \left(k\right)\parallel \right)\right\}\hfill \\ \hfill \le & \sum _{s=1}^{N-1}\left\{\overline{\lambda }\left(Q\right)\parallel {\Phi }_{\eta }^{s-1}{F}^{-1}\varpi \left(k\right)\parallel \left(2\overline{\lambda }\left(Q\right)\parallel h_{\eta }-\varphi (k+s|k)\parallel +\overline{\lambda }\left(Q\right)\parallel {\Phi }_{\eta }^{s-1}{F}^{-1}\varpi \left(k\right)\parallel \right)\right.\hfill \\ \hfill & +\overline{\lambda }\left(R\right)\parallel K_{\eta }{\Phi }_{\eta }^{s-1}{F}^{-1}\varpi \left(k\right)\parallel \left(2\overline{\lambda }\left(R\right)\parallel h_{c_{\eta }}-{\left[({\tilde{h}}_{e_d}+\overline{d})0\right]}^T-\chi (k+s|k)\parallel \right.\hfill \\ \hfill & \left.\left.+\overline{\lambda }\left(R\right)\parallel K_{\eta }{\Phi }_{\eta }^{s-1}{F}^{-1}\varpi \left(k\right)\parallel \right)\right\}\hfill \\ \hfill \le & \sum _{s=1}^{N-1}\left\{{\overline{\lambda }}^2\left(Q\right){\Omega }^2\left(s\right)+2{\overline{\lambda }}^2\left(Q\right)\Omega \left(s\right)F\left(s\right)+{\overline{\lambda }}^2\left(R\right){\Xi }^2\left(s\right)+2{\overline{\lambda }}^2\left(R\right)\Xi \left(s\right)\Psi \left(s\right)\right\},\hfill \end{array}$$

where

$$\begin{array}{cc}\hfill & \Omega \left(s\right)\triangleq \underset{\varpi \left(k\right)\in \mathcal{W}}{max}\parallel {\Phi }_{\eta }^{s-1}{F}^{-1}\varpi \left(k\right)\parallel ,F\left(s\right)\triangleq \parallel h_{\eta }-\varphi (k+s|k)\parallel ,\hfill \\ \hfill & \Xi \left(s\right)\triangleq \underset{\varpi \left(k\right)\in \mathcal{W}}{max}\parallel K_{\eta }{\Phi }_{\eta }^{s-1}{F}^{-1}\varpi \left(k\right)\parallel ,\Psi \left(s\right)\triangleq \parallel h_{c_{\eta }}-|{\tilde{h}}_{e_d}+\overline{d}|{\left[b_{u}0\right]}^T-\chi (k+s|k)\parallel .\hfill \end{array}$$

For ${\Delta }_2$, the following can be obtained:

$$\begin{array}{cc}\hfill {\Delta }_2=& \parallel \tilde{\eta }{(k+N|k+1)\parallel }_Q^2+\parallel K_{\eta }\tilde{\eta }{(k+N|k+1)\parallel }_R^2+{\parallel \tilde{\eta }(k+1+N|k+1)\parallel }_P^2\hfill \\ \hfill & -\parallel {\widehat{\eta }}^{*}{(k+N|k)\parallel }_P^2+\parallel \tilde{\eta }{(k+N|k+1)\parallel }_P^2-{\parallel \tilde{\eta }(k+N|k+1)\parallel }_P^2\hfill \\ \hfill \le & (\parallel \tilde{\eta }{(k+N|k+1)\parallel }_P^2-\parallel {\widehat{\eta }}^{*}{(k+N|k)\parallel }_P^2)-\parallel \tilde{\eta }{(k+N|k+1)\parallel }_{\overline{Q}}^2\hfill \\ \hfill \le & \parallel \tilde{\eta }(k+N|k+1)-{\widehat{\eta }}^{*}{(k+N|k)\parallel }_P(\parallel \tilde{\eta }{(k+N|k+1)\parallel }_P+\parallel {\widehat{\eta }}^{*}(k+N|k){\parallel }_P)\hfill \\ \hfill & -\parallel \tilde{\eta }{(k+N|k+1)\parallel }_{\overline{Q}}^2\hfill \\ \hfill \le & -\parallel \tilde{\eta }{(k+N|k+1)\parallel }_{\overline{Q}}^2+\overline{\lambda }\left(P\right)\parallel {\Phi }_{\eta }^{N-1}{F}^{-1}\varpi \left(k\right)\parallel \hfill \\ \hfill & \times \left(2\overline{\lambda }\left(P\right)\parallel {\widehat{\eta }}^{*}(k+N|k)\parallel +\overline{\lambda }\left(P\right)\parallel {\Phi }_{\eta }^{N-1}{F}^{-1}\varpi \left(k\right)\parallel \right),\hfill \end{array}$$

where $-\overline{Q}={\Phi }_{\eta }^{T}P{\Phi }_{\eta }+Q+K_{\eta }^{T}R{K}_{\eta }-P$ and $\pi \left(N\right)\triangleq min\left\{F\right(N),\Psi (N\left)\right\}$. Letting $-\overline{Q}<0$, we obtain the following:

$$\begin{array}{c}\hfill {\Phi }_{\eta }^{T}P{\Phi }_{\eta }+Q+K_{\eta }^{T}R{K}_{\eta }-P<0\end{array}$$

which can be equivalently converted to LMI (40), such that K and $L_d$ can be obtained by solving LMI (40). Thus, $\parallel \tilde{\eta }{(k+N|k+1)\parallel }_{\overline{Q}}^2>0$, and we have the following:

$$\begin{array}{c}\hfill {\Delta }_2\le {\overline{\lambda }}^2\left(P\right){\Omega }^2\left(N\right)+\overline{\lambda }\left(P\right)\Omega \left(N\right)\pi \left(N\right).\end{array}$$

For ${\Delta }_3$, there exists ${\Delta }_3\le -{\parallel {\widehat{\eta }}^{*}\left(k\right|k)\parallel }_Q^2$. Combining ${\Delta }_1$, ${\Delta }_2$, and ${\Delta }_3$, we obtain the following:

$$\begin{array}{c}\hfill \Delta J\left(k\right)\le -\parallel {\widehat{\eta }}^{*}{\left(k\right|k)\parallel }_Q^2+\overline{\sigma }.\end{array}$$

with

$$\begin{array}{cc}\hfill \overline{\sigma }\triangleq & \sum _{s=1}^{N-1}\left\{{\overline{\lambda }}^2\left(Q\right){\Omega }^2\left(s\right)+2{\overline{\lambda }}^2\left(Q\right)\Omega \left(s\right)F\left(s\right)+{\overline{\lambda }}^2\left(R\right){\Xi }^2\left(s\right)+2{\overline{\lambda }}^2\left(R\right)\Xi \left(s\right)\Psi \left(s\right)\right\}\hfill \\ \hfill & +{\overline{\lambda }}^2\left(P\right){\Omega }^2\left(N\right)+\overline{\lambda }\left(P\right)\Omega \left(N\right)\pi \left(N\right).\hfill \end{array}$$

According to the optimality principle, it is easy to obtain the following:

$$\begin{array}{cc}\hfill J({\widehat{\mathbf{ȷ}}}^{*}(k+1),{\widehat{\mathbf{c}}}_{\eta }^{*}(k+1))-J({\widehat{\mathbf{ȷ}}}^{*}\left(k\right),{\widehat{\mathbf{c}}}_{\eta }^{*}\left(k\right))\le & J(\tilde{\mathbf{ȷ}}(k+1),{\tilde{\mathbf{c}}}_{\eta }(k+1))-J({\widehat{\mathbf{ȷ}}}^{*}\left(k\right),{\widehat{\mathbf{c}}}_{\eta }^{*}\left(k\right))\hfill \\ \hfill \le & -\parallel {\widehat{\eta }}^{*}{\left(k\right|k)\parallel }_Q^2+\overline{\sigma }\hfill \\ \hfill =& -\parallel \check{\eta }{\left(k\right)\parallel }_Q^2+\overline{\sigma }.\hfill \end{array}$$

Therefore, $J({\widehat{\mathbf{ȷ}}}^{*}\left(k\right),{\widehat{\mathbf{c}}}_{\eta }^{*}\left(k\right))$ can be used as the Lyapunov function. As a result, $\check{\eta }\left(k\right)$ can converge to the set $\mathbb{S}$ as $k\to +\infty$, suggesting that system (22) is uniformly bounded, such that the uniform boundedness of system (11) and error system (18) can also be guaranteed. □

Remark 4. 

Compared to existing smooth and precise braking methods [13,14,15], the RHESO-based MPC proposed in this paper has the following advantages: (1) A new neural network-based vehicle kinematic model is proposed, which simulates the dynamic characteristics of vehicles through neural networks to supplement the shortcomings of the kinematic model in accurately describing the throttle and brake control on the vehicle speed; (2) the proposed RHESO can estimate the unmodeled dynamics and potential disturbances in AV kinematic modeling, providing preliminary theoretical support for the AV to achieve smooth and accurate preset point braking; and (3) design an AV speed control method based on MPC, taking into account the inherent speed and control constraints of AV to enable smooth and accurate braking at preset points.

4. Numerical Example

An electric AV model provided by the CarSim software is employed to obtain the associated data of the throttle/brake on vehicle speed and acceleration. Therefore, the corresponding approximate dynamic model can be further obtained by using the proposed NN. Partial parameters of the electric AV used in this example are given in

Table 2. Partial parameters of the electric AV.

Parameter	Value
Sprung mass	1270 kg
Wheelbase	2910 mm
Vehicle height	1730 mm
Vehicle width	2082 mm
Roll inertia	536.6 kg·${\mathrm{m}}^2$
Pitch inertia	1536.7 kg·${\mathrm{m}}^2$
Yaw inertia	1536.7 kg·${\mathrm{m}}^2$
Brake torque at front wheel	250 N·$\mathrm{m}/\mathrm{Mpa}$
Brake torque at rear wheel	150 N·$\mathrm{m}/\mathrm{Mpa}$

.

As this paper focuses on longitudinal control, the driving angle $\overline{\theta }={60}^{°}$ is selected in (2). The external disturbance is selected as $\epsilon \left(t\right)=sin(t+0.2)$. Other simulation parameters are given in

Table 3. Parameters used in the simulation.

Parameter	$\Lambda$	$\kappa$	$\overline{\epsilon }$	H
Value	0.5	1.19	3.5	0
Parameter	${\widehat{W}}_1^{*}$	${\widehat{W}}_2^{*}$	$a_{\sigma }$	$b_{\sigma }$
Value	−0.21	−1.58	1	1
Parameter	${\tilde{h}}_{1,min,\xi }$	${\tilde{h}}_{1,max,\xi }$	${\tilde{h}}_{2,min,\xi }$	${\tilde{h}}_{2,max,\xi }$
Value	−100 m	100 m	−100 m	100 m
Parameter	${\tilde{h}}_{3,min,\xi }$	${\tilde{h}}_{3,max,\xi }$	${\tilde{h}}_{min,u}$	${\tilde{h}}_{max,u}$
Value	0 m/s	50 m/s	0	9
Parameter	N	Q	R	T
Value	10	$150I$	I	0.1 s

.

The activation function vector and the weight matrix for the NN are reconstructed as follows:

$$\begin{array}{cc}\hfill \sigma \left(t\right)={\left[v\left(t\right)u_o\left(t\right)1\right]}^T,{\widehat{W}}^{*}=& {[-0.21-1.581.09]}^T.\hfill \end{array}$$

The matrices $L_d$, K, and P are obtained by solving LMI (40) as follows:

$$\begin{array}{cc}\hfill L_d=& [00-6.1],K=\left[1.42.63.6\right],\hfill \\ \hfill P=& \left[\begin{array}{cccc}2576.1& -1376.5& 18.3& 0\\ -1376.5& 1003.7& 33.9& 0\\ 18.3& 33.9& 43.9& 0\\ 0& 0& 0& 20.9\end{array}\right].\hfill \end{array}$$

The initial state of the electric AV is selected as follows:

$$\begin{array}{c}\hfill \xi \left(0\right)={[-15.1-26.28]}^T.\end{array}$$

To display the effectiveness of the proposed RHESO-based MPC for achieving smooth braking at preset points in AVs, its numerical simulation result is compared with those obtained by the traditional ESO-based MPC, active disturbance rejection control (ADRC), and robust MPC (RMPC). The simulation results are displayed in Figure 3, Figure 4, Figure 5 and Figure 6. Figure 3 shows the position response of the electric AV using the proposed RHESO-based MPC, as well as traditional ESO-based MPC, ADRC, and RMPC. It can be clearly seen that the position components of the AV, i.e., $x\left(t\right)$ and $y\left(t\right)$, enable it to converge well to the origin when the proposed RHESO-based MPC is employed. However, if the traditional ESO-based MPC, ADRC, and RMPC are adopted, the AV’s position components $x\left(t\right)$ and $y\left(t\right)$ cannot converge to the origin. Similarly, Figure 4 exhibits the velocity response of the electric AV obtained by the proposed RHESO-based MPC, traditional ESO-based MPC, ADRC, and RMPC. The simulation results in Figure 3 and Figure 4 reveal that the proposed RHESO-based MPC allows for achieving smooth braking at designated locations in the electric AV. On the contrary, smooth braking at designated positions in the electric AV is difficult to achieve if the traditional ESO-based MPC, ADRC, and RMPC are applied. Figure 5 plots the control input of the electric AV. Furthermore, the history of estimation errors obtained by the proposed RHESO and the traditional ESO are shown in Figure 6. Note that the error convergence speed of the traditional ESO is faster than that of the proposed RHESO, but to some extent, the existence of optimized estimation errors makes the electric AV’s braking process smoother and more accurate, which is another advantage of the proposed RHESO-based MPC.

In the following, simulation results under another set of parameters are shown. The initial conditions are given as follows:

$$\begin{array}{c}\hfill \xi \left(0\right)={[-30-17.318]}^T,\overline{\theta }={30}^{\circ },\epsilon \left(t\right)=5sin\left(0.2t\right)+cos\left(t\right).\end{array}$$

Figure 7, Figure 8, Figure 9 and Figure 10 show the impacts of the proposed RHESO-MPC and traditional ESO-based MPC, ADRC, and RMPC on AV when the initial distance and initial speed are relatively large. From Figure 7, Figure 8, Figure 9 and Figure 10, it is easy to observe that the RHESO-MPC outperforms traditional ESO-based MPC, ADRC, and RMPC in terms of position control accuracy and speed convergence accuracy when both the initial position and speed increase. As a conclusion, from Figure 3, Figure 4, Figure 5, Figure 6, Figure 7, Figure 8, Figure 9 and Figure 10, the effectiveness of the proposed method is verified.

5. Conclusions

This paper analyzed the problem of smooth braking at preset points in AVs. Given the difficulty of obtaining dynamic model parameters for AVs and the limitations of using kinematic models for servo control, an NN was designed to approximate the velocity order model within the dynamic model, enabling the use of throttle/brake openings as control signals. A new extended kinematic model of AVs was proposed by integrating the proposed NN model with the traditional kinematic model, which facilitates longitudinal speed servo control. An RHESO-based MPC method was also proposed to achieve smooth braking at preset points in AVs. Furthermore, sufficient conditions were established to guarantee the recursive feasibility of the optimization problem and the practical stability of the closed-loop system. Finally, the effectiveness of the proposed method was validated through an illustrative numerical simulation, demonstrating that the proposed extended kinematic model can effectively control velocity and position and that the RHESO-based MPC method ensures smooth braking at preset points.