%===============================================================================
% ifacconf.tex 2022-02-11 jpuente  
% 2022-11-11 jpuente change length of abstract
% Template for IFAC meeting papers
% Copyright (c) 2022 International Federation of Automatic Control
%===============================================================================
\documentclass{ifacconf}

%\usepackage[numbers]{natbib}
%\usepackage[authoryear]{natbib}
%\usepackage[authoryear,longnamesfirst]{natbib}
%\usepackage{caption}
%\usepackage{booktabs}
%\usepackage{longtable}
%\usepackage{multirow}
%\usepackage{multicol}
%\usepackage{float}
\usepackage{subfigure}
\usepackage{algorithm}
\usepackage{algpseudocode}
\usepackage{amsmath}
%\usepackage{amsthm}
\usepackage{amssymb}
%\usepackage{amsfonts}
\usepackage{graphicx}
%\usepackage{mathrsfs}
%\usepackage{url}
%\usepackage{amsthm}
%\captionsetup[figure]{labelfont={bf}, labelformat={default}, labelsep=period, name={Fig.}}
%\usepackage{color}
%\usepackage{ctex}
%% hyper link
%\usepackage{hyperref}
%\hypersetup{colorlinks=true, citecolor=blue, linkcolor=blue, urlcolor=blue}
%\usepackage{graphicx}      % include this line if your document contains figures
\usepackage{natbib}        % required for bibliography
%===============================================================================

%\newtheorem{theorem}{Theorem}
\newtheorem{lemma}{Lemma}
%\newtheorem{algorithm}{Algorithm}
\newtheorem{definition}{Definition}
\newtheorem{remark}{Remark}
\newtheorem{problem}{Problem}

\begin{document}
\begin{frontmatter}

\title{Predictive Control via Augmented Lagrange Function for Autonomous Vehicles Trajectory Tracking with Dynamic Quantization}
%\title{Style for IFAC Conferences \& Symposia: Use Title Case for
%  Paper Title\thanksref{footnoteinfo}} 
% Title, preferably not more than 10 words.
%
%\thanks[footnoteinfo]{Sponsor and financial support acknowledgment
%goes here. Paper titles should be written in uppercase and lowercase
%letters, not all uppercase.}

\author[First]{Zhaojin Yu} 
\author[First]{Xiaoming Tang} 
\author[Second,Third]{Xiao Lv}
\author[First]{Yongzhen Cao}

\address[First]{College of Automation, 
   Chongqing University of Posts and Telecommunications, Chongqing 400065, China (e-mail: txmmyeye@126.com).}
\address[Second]{Chongqing Special Equipment Inspection and Research Institute, Chongqing 401121, China}
\address[Third]{Key Laboratory of Electromechanical Equipment Security in Western Complex Environment for State Market Regulation, Chongqing 401121, China}

\begin{abstract}                % Abstract of 50--100 words
This paper investigates the Augmented Lagrange Function (ALM) based model predictive control (MPC) strategy for autonomous vehicle trajectory tracking. Firstly, a linear time-varying error model with measurement uncertainty and dynamic quantization is established. Measurement uncertainty indicates that vehicle state extraction fails due to the jittering of sensors during actual driving. The dynamic quantizer, which can adjust the quantization parameters online, is used to quantize the control input signals. The stability condition of the closed-loop system is obtained through the Lyapunov function and solved using linear matrix inequality technology. Then, according to the stability condition, the model predictive controller is designed by solving a “min-max” optimization problem that is based on a cost function over the finite time horizon. In order to solve the MPC optimization problem, the regularized least-squares method with ALM and an online iterative algorithm are explored, which obtain the analytical solution of the optimization problem. Finally, the effectiveness of the designed controller is verified by simulation experiments which show that the ALM has a very accurate computational value.
\end{abstract}

\begin{keyword}
Model predictive control (MPC), Augmented Lagrange Function (ALM), Trajectory tracking, Measurement uncertainty, Dynamic quantization.
\end{keyword}

\end{frontmatter}
%===============================================================================

\section{Introduction}%\label{sec1}

Autonomous driving technology for traffic accidents prediction is being of great significance to alleviate the safety pressure of regional traffic management \cite{ye2023data,xie2020cooperative,liu2022optimal}. The continuous development of computer technology and artificial intelligence has made the emergence of autonomous driving technology possible \cite{ning2021survey,gao2021autonomous,fernandez2023trustworthy,grigorescu2020survey,kim2022autonomous}. The trajectory tracking is playing an important role in autonomous driving \cite{gw2021autonomous}, and various tracking control algorithms have been proposed for autonomous vehicles, such as PID control, fuzzy control, reinforcement learning, and model predictive control (MPC) \cite{zou2018dynamic,chen2022speed,liu2021reinforcement,choi2021horizonwise}.

MPC is a popular control technique which solves an optimization problem at each time step to determine the optimal control sequence and allows explicit incorporation of the plant uncertainty in the problem formulation \cite{kothare1996robust}. Hence, many works have begun to use MPC approach to solve trajectory tracking problem. In order to improve the performance of trajectory tracking, some practical factors like obstacles, velocity, yaw angle, rollover, and so on are always considered in MPC optimization problems. \cite{guo2018simultaneous} studied obstacle avoidance performance by addressing the MPC optimization problem considering the lateral position, velocity and yaw angle. To improve the path tracking performance of autonomous vehicles, \cite{wang2021path} designed a particle swarm optimization algorithm that considered the lateral offset, the steering frequency and the real-time of the MPC algorithm to realize adaptive optimization for the predictive horizon. In \cite{chu2022trajectory}, MPC with PID feedback control law presented improved performance in tracking accuracy and steering smoothness in vehicle cornering. To predict and avoid rollover, \cite{ghazali2017vehicle} limited rollover as a hard constraint in the MPC problem by using a future error estimation approach. The authors of \cite{ye2019linear} added relaxation factors (soft constraints) to the MPC optimization problem so that the speed and steering wheel can be controlled in phase. Reference \cite{vu2021model} presented a nonlinear model predictive controller considering the dynamic constraints of the vehicle’s physical limitations, the environmental conditions, and the surrounding obstacles to improve the ability of autonomous vehicles. Although these investigations obtained some very nice results, the model predictive controllers were more complex and therefore difficult to solve since many factors in autonomous vehicles have been considered.

In autonomous driving, complex autonomous driving algorithms require communication networks to transmit information \cite{luan2020trajectory,jo2014development}. However, the insertion of communication networks may impact the stability and trajectory tracking accuracy of autonomous vehicles since the adverse effects introduced by network may be involved. In order to address network-induced constraints such as quantization error, transmission delay, and data dropouts, \cite{chang2022quantized} proposed a robust gain-scheduling {$H_\infty$} output feedback controller. In \cite{gu2021path}, a learning-based event-triggered mechanism was introduced to relieve the burden of the autonomous vehicle communication network. Reference \cite{ye2023robust} applied the Bernoulli random distribution approach to deal with data packet dropout and quantized the output signal with a dynamic quantizer. \cite{zhao2022robust} utilized a generalized lumped delay form to unify the time-varying data dropout and network-induced delay, and then a Markovian process is presented to describe the lumped delay as a stochastic distribution. These studies provided us with ideas for applying the approach of MPC in autonomous driving to communication networks. In a networked environment where both data quantization and packet loss may occur, MPC needs to ensure robust stability at each sampling time. Notably, some excellent methods have been applied to solve MPC optimization problems with uncertain data and quantization. In \cite{liu2020robust} a regularized method based on the Lagrange Multiplier Method (LM) was used to obtain the input with an analytical expression for the least-squares problem with uncertain parameters of MPC. Although there have been some good research results in the field of MPC for autonomous vehicle trajectory tracking, the solutions to such problems that simultaneously consider both dynamic quantization and uncertain data are few, especially the calculation methods that have an analytical expression for the solution of this type of MPC optimization problem, which is the motivation for this research.

We focus on improving the computational method for MPC optimization problem of autonomous vehicle trajectory tracking while considering dynamic quantization and measurement deviation. The following summarizes the main contributions of this paper:

$\left. \rm{1}\right)$ According to the vehicle kinematic model, a time-varying trajectory tracking problem based on MPC is constructed, containing dynamic quantization and Bernoulli variables modeled by measurement deviation.

$\left. \rm{2}\right)$ The parameters of the dynamic quantizer for the control input are designed. Then, in order to ensure closed-loop stability of the system, the dynamic parameters of the quantizer and the Bernoulli variables are introduced into the Lyapunov stability condition when designing the terminal weight matrix of the cost function online at each time.

$\left. \rm{3}\right)$ An online iterative model predictive controller is developed for time-varying systems, and the “min-max” optimization problem is solved by the ALM to obtain the analytical solution of input that is convenient for calculation.

The remaining sections are arranged as follows: Section 2 presents the problem and performs modeling. In Section 3, the performance issue involving dynamic quantization is analyzed. Section 4 gives the results of two different trajectory numerical simulation experiments. The final conclusion is drawn in Section 5.

\section{System Description}%\label{sec2}

\subsection{Kinematics Model of Autonomous Vehicle}

\begin{figure}[h]
	\centering
	\includegraphics{model.pdf}
	\caption{Representation of the kinematics model.}
	\label{fig1}
\end{figure}
%
%\begin{figure*}
%	\centerline{\includegraphics[width=1pt,height=1pc,draft]{empty}}
%	\caption{This is the sample figure caption.\label{fig2}}
%\end{figure*}

In a planar Cartesian coordinate system, the kinematic relationship of the autonomous vehicle is shown in Figure \ref{fig1}, where $\varphi$ and $\delta$ are the heading angle and the front wheel deflection angle of the vehicle, respectively. Let $e=[x_e~ y_e~ \varphi_e]^{\rm T}$ describes the actual state of the autonomous vehicle, and $d=[x_d~ y_d~ \varphi_d]^{\rm T}$ describes the desired position. Then, by analyzing the kinematic relationships, we have:
\begin{align}\label{eq1}
	\left[
	\begin{array}{c}
		\dot x_d \\
		\dot y_d \\
		\dot \varphi_d \\
	\end{array}
	\right]=\left[
	\begin{array}{cc}
		\cos \varphi_d \\
		\sin \varphi_d \\
		\dfrac{\tan \delta_d}{l} \\
	\end{array}
	\right]v_d,
\end{align}
where $v_d$ is the desired linear velocity. It is assumed that there is a predetermined path in space. Consider an autonomous vehicle that tracks the path. The parameter $x$ is the tracking error. Then, we establish an equation based on error as follows:
\begin{align}\label{eq2}
	\dot x=\left[
	\begin{array}{c}
		\dot x_e-\dot x_d \\
		\dot y_e-\dot y_d \\
		\dot \varphi_e-\dot \varphi_d \\
	\end{array}
	\right]=f(x_d,y_d)
\end{align}	
The error equation at this point is only a preliminary result obtained after differentiation, and it is a nonlinear expression. We can linearize the system by using first-order Taylor expansion, so that \eqref{eq2} is rewritten in the form of \eqref{eq3}:
\begin{align}\label{eq3}
	\dot {x}=\left[
	\begin{array}{ccc}
		0 & 0 & -v_d\sin \varphi_d \\
		0 & 0 & v_d\cos \varphi_d \\
		0 & 0 & 0 \\
	\end{array}
	\right]x+\left[
	\begin{array}{cc}
		\cos \varphi_d & 0 \\
		\sin \varphi_d & 0 \\
		\dfrac{\tan \delta_d}{l} & \dfrac{v_d}{l\cos^{2} \delta_d} \\
	\end{array}
	\right]u
\end{align}
where $x=[x_e-x_d~ y_e-y_d~ \varphi_e-\varphi_d]^{\rm T},~u=[v_e-v_d~ \delta_e-\delta_d]^{\rm T}$. It is evident that both of the two coefficient matrices in \eqref{eq3} have time-varying properties.

After receiving the status information about the autonomous vehicle, the controller processes it into control signals and sends them to the vehicle. In order to facilitate the transmission and processing of data, a dynamic quantizer has been added to the controller. In addition, inaccurate measurements due to sensor defects and vehicle movements are known as uncertain communications, which may affect system stability if left untreated. Therefore, we process this uncertainty by constructing Bernoulli sequences. After introducing dynamic quantization and data uncertainties, the state equation is discretized by the forward Euler method, and system \eqref{eq3} can be converted to \eqref{eq4}.
\begin{align}\label{eq4}
	x(k+1)=A_{k}x(k)+\gamma (k)B_{k}q_{u}(u(k))
\end{align}
where
\begin{align*}
	A_{k}=\left[
	\begin{array}{ccc}
		1 & 0 & -v_dT\sin \varphi_d \\
		0 & 1 & v_dT\cos \varphi_d \\
		0 & 0 & 1 \\
	\end{array}
	\right],~B_{k}=\left[
	\begin{array}{cc}
		T\cos \varphi_d & 0 \\
		T\sin \varphi_d & 0 \\
		\dfrac{\tan \delta_d}{l} & \dfrac{Tv_d}{l\cos^{2} \delta_d} \\
	\end{array}
	\right],
\end{align*}
$T$ is the sample time, and $\mathit{\gamma} (k)=1$ is a Bernoulli sequence modeled by uncertainty of measurement \cite{ye2023robust}. Define $\gamma (k)=1$, if the pose capture is successful; otherwise, $\gamma (k)=0$. In this paper, the dynamic quantizer $q_{u}(u)$ is used to quantize the control input $U(k)$.

\subsection{Dynamic Quantizer}

Assuming $x\in\mathbb{R}$ is the variable that needs to be quantized. Usually, what we call a quantizer is a description of a piecewise function $q\in\mathbb{R}\rightarrow\mathbb{R}_{i}$, and this piecewise function is a constant function. Therefore, the set $\mathbb{R}$ is divided into a finite quantity of quantization regions of the form $\{x\in\mathbb{R}:q(x)=x_{i}\},x_{i}\in\mathbb{R}_{i}$. For two positive real numbers, $\Delta$ and $M$, the following conditions hold:
\begin{align}\label{eq5}
	\left| q(x)-x \right|\leq \Delta,~~~\rm{if}~\left| \it x\right|\leq\it M;
\end{align}
and
\begin{align}\label{eq6}
	\left| q(x) \right| > M-\Delta,~~~\rm{if}~\left| \it x\right|>\it M
\end{align}
where $M$ and $\Delta$ are the range of the quantizer and the quantization error, respectively. Among these two inequalities, \eqref{eq5} represents the condition where the quantizer isn't saturated, and \eqref{eq6} can be used to detect whether the quantizer is saturated.

The above two conditions demonstrate a static quantization method \cite{liberzon2003hybrid}. In order to make the measurement data more accurate, this article considers a dynamic quantization strategy in the following form:
\begin{align}\label{eq7}
	q_{x}(x)=\mu q(\frac{x}{\mu})
\end{align}
where $q_{x}(x)$ is the dynamic quantizer, and $q$ is one way of the static quantization mentioned earlier. Subsequently, the dynamic quantization range and the error of the quantizer are represented as $M_{x}$ and $\Delta_{x}$, respectively. The variable $\mu>0$ represents the dynamic parameter. Through analysis, it is found that increasing or decreasing the value of $\mu$ can zoom in or zoom out $M_{x}$ and $\Delta_{x}$. So that, $\mu$ can also be considered a zoom factor. By the way, the parameter $\mu$ is at the core of the dynamic quantizer, and we need to design its value based on actual problems. From \eqref{eq5} and \eqref{eq6}, \eqref{eq7} is transformed into a form of error:
\begin{align}\label{eq8}
	&\left| q_{x}(x)-x \right|\leq\mu\Delta_{x},~~~if~\left| \it x\right|\leq\mu\it M_{x}\nonumber\\
	&\left| q_{x}(x)-x \right| > \mu\Delta_{x},~~~if~\left| \it x\right|>\mu\it M_{x}
\end{align}

Similarly, \eqref{eq8} describes the situations of dynamic quantizer unsaturation and saturation. The parameters $\mu\Delta_{x}$ and $\mu\it M_{x}$, respectively, represent the boundary of the quantization error and the quantization range boundary. When the quantizer is not saturated, the signal $x$ is within the quantization range $\mu\it M_{x}$, and when the quantizer is saturated, the quantization error will be larger than $\mu\Delta_{x}$.

\begin{remark}\label{rem1}
	In a static quantization strategy, the quantized level of the quantizer is related to the quantization density. Since a static quantizer only corresponds to one quantization density, its quantization level remains constant. The self-designed dynamic parameter $\mu$ in the dynamic quantizer can be adjusted to optimize the quantized level. Hence, compared to a static quantizer, it can obtain a larger region of attraction. Here, we will adopt a dynamic quantization strategy to handle control input.
\end{remark}

\subsection{Discrete System with Quantization}

Substituting \eqref{eq7} into \eqref{eq4}, and considering \eqref{eq8}, obtain a system equation considering quantization:
\begin{align}\label{eq9}
	x(k+1)=A_{k}x(k)+\gamma (k)B_{k}\left[ u(k)+\hat u(k)\right] 
\end{align}
where $\hat u(k)=\mu \left( q\left( \dfrac{u(k)}{\mu}\right) -\dfrac{u(k)}{\mu}\right)$, representing the quantization error of $u(k)$. The variable $\hat u(k)$ is an add-on to the quantizer.

From the system description, we can see that a tracking error model has been built. The ultimate goal is to enable the autonomous vehicle to travel along the ideal track within the site. This problem can be translated into enabling system \eqref{eq9} to ultimately achieve stability. This requires the proper control input to be designed. To solve the stability problem of the system, we consider robust predictive control and write such a cost function in the form of quadratic as follows:
\begin{align}\label{eq10}
	J(k)=&\sum_{i=0}^{p-1} \left[ x^{\rm T}(k+i|k)Q_{\mathrm{1}}x(k+i|k)+u^{\rm T}(k+i|k) \right.  \\ \nonumber
	&\left. Q_{\mathrm{2}}u(k+i|k)\right]+x^{\rm T}(k+p|k)Wx(k+p|k)
\end{align}
where the matrices $Q_{\mathrm{1}}, Q_{\mathrm{2}}$ and $W$ in quadratic forms are the weighted matrices of the item in which they are located. $Q_{\mathrm{1}}$ and $Q_{\mathrm{2}}$ represent the weights of variable $x$ and control variable $U(k)$, and $W$ corresponds to the weight of the terminal penalty term outside the summation term. The size of $W$ has a significant impact on the stability of the system, so it is necessary to additionally solve the weight matrix $W$. The parameter $p$ represents the predictive time domain of the controller. The $k+i|k$ in brackets indicates the state of the system that predicts time $k+i$ at time $k$.

To construct the optimization problem, the Bernoulli sequence and quantization are further processed here. The probability of successful data transmission is expressed in $\bar{\boldsymbol{\gamma}}$. Then $\bar{\boldsymbol{\gamma}}$ is the mathematical expectation of the Bernoulli sequence $\gamma (k)$. In other words, $E\left\lbrace \gamma (k)\right\rbrace = \bar{\boldsymbol{\gamma}}$. Furthermore, appropriate dynamic parameter $\mu$ needs to be designed. Therefore, a robust predictive controller considering dynamic quantization is used to solve the trajectory tracking problem of autonomous vehicles, as follows:

\begin{problem}\label{Pro1}
	Given the probability of successful data transmission $\bar{\boldsymbol{\gamma}}$, the dynamic quantization range $M_{u}$, the error of the quantizer $\Delta_{u}$ and the reference trajectory, the optimal control input of system \eqref{eq9} can be obtained by solving the optimization problem as follows:
	\begin{align}\label{eq11}
		&u^{*}(k)=\mathop{arg}\mathop{min}\limits_{u(k)}\mathop{max}\limits_{ \Delta_{u},M_{u} }E\left\lbrace J(k)\right\rbrace \nonumber\\
		&s.t.~\eqref{eq9}~and~\eqref{eq10}. 
	\end{align}
	where
	\begin{align}\label{eq12}
		&x(k+1+i|k)=A_{k+i,k}x(k+i|k) \\ \nonumber
		&+\gamma (k+i|k)B_{k+i,k}(u(k+i|k)+\hat u(k+i|k)).
	\end{align}
\end{problem}
In order to simplify the calculation, we make the following assumption:
\begin{align*}
	A_{k+i,k}=A_{k},~B_{k+i,k}=B_{k}.
\end{align*}

\begin{remark}
	We know that classical model predictive control methods require solving optimization problems with LMI constraints at each sampling time. Hence, the computation time for each sampling can be very long, especially in cases where matrix dimensions are large. Since this article is concerned with trajectory tracking problems, in order to ensure the real-time performance of the controller, we will design a predictive controller based on the robust regularized least squares approach \cite{sayed2002regularized}.
\end{remark}

At the end of this section, a definition and a lemma are given, mainly aimed at the design of the terminal weighting matrix $W$ and the dynamic parameter $\mu$, and the computation of the solving process. Design $W$ and $\mu$ based on the following definition:

\begin{definition}\label{Def1}
	The system \eqref{eq9} is mean squared asymptotically stable if $\exists S \in \mathbb{R}^{\rm 2\times3}$, let the terminal weighting matrix $W$ and the dynamic quantizer $q_{u}$ satisfy the following inequality:
	\begin{align}\label{eq13}
		&E\left\lbrace \left\| A_{k}+\gamma (k)B_{k}q_{u}(S) \right\|_{W}^{\rm 2}\right\rbrace+Q_{\mathrm{1}}+\left\| S\right\|_{Q_{\mathrm{2}}}^{\rm 2} - W \\ \nonumber
		&\leq 0, k\in\left[ 1,+\infty\right)
	\end{align}
\end{definition}

\begin{lemma}[\cite{wang1992robust}]
	\label{lem1}
	$\forall a,b\in \mathbb{R}^{n}$, and $X$ is a positive definite matrix that has compatible dimensions with $a$ and $b$, the following inequality holds:
	\begin{align}\label{eq14}
		2a^{\rm T}b \leq a^{\rm T}Xa+b^{\rm T}X^{\mathrm{-1}}b
	\end{align}
\end{lemma}

\section{Main results}%\label{sec3}

For the prediction equation in \eqref{eq12}, the Bernoulli variable should have been $\gamma(k+i|k)$. But because it is an IID variable, we use a zero-order holder strategy for it: $\gamma(k+i|k)=\gamma (k)$. According to \eqref{eq8}, the parameter $\mu$ satisfies the following inequality condition:
\begin{align*}
	\left\| \dfrac{u(k)}{\mu}\right\| \leq M_{u}
\end{align*}
So the dynamic parameter of the quantizer $q_{u}(u(k)$ can be adjusted to:
\begin{align}\label{eq15}
	\mu = \dfrac{\left\| u(k)\right\| }{\sqrt{\alpha_{u}}M_{u}}
\end{align}
where $0<\alpha_{u}<1$ \cite{zheng2021quantized}, its value needs to be further determined. And then, for the quantization error $\hat u(k)$ in system \eqref{eq9}, it has:
\begin{align*}
	\left\| \hat u(k) \right\| = \mu \left\| q\left( \dfrac{u(k)}{\mu}\right) -\dfrac{u(k)}{\mu}\right\| \leq \mu \Delta_{u}
\end{align*}
Bring \eqref{eq15} into the above inequality:
\begin{align}\label{eq16}
	\left\| \hat u(k)\right\| \leq \dfrac{\Delta_{u} }{\sqrt{\alpha_{u}}M_{u}} \left\| u(k)\right\|
\end{align}
\eqref{eq16} is equivalent to:
\begin{align}\label{eq17}
	\hat u^{\rm T}(k)\hat u(k) \leq \dfrac{\Delta_{u}^{\rm 2} }{\alpha_{u}M_{u}^{2}} u^{\rm T}(k)u(k)
\end{align}
Based on the above derivation, introducing the prediction equation \eqref{eq12} into the cost function \eqref{eq10} yields a completely new equality:
\begin{align}\label{eq18}
	J(k) =& \left\| \Phi x(k) + \gamma (k)F\left[ U(k)+\hat{U}(k)\right] \right\|_{Q_{w}}^{2} \\ \nonumber
	&+ \left\| U(k) \right\|_{Q}^{2}
\end{align}
where
\begin{align*}
	&Q_{w} = {\rm diag} \{\underbrace {Q_{\mathrm{1}}, Q_{\mathrm{1}}, \cdots, Q_{\mathrm{1}}}_{p}, W \},Q = {\rm diag} \{\underbrace {Q_{\mathrm{2}}, Q_{\mathrm{2}}, \cdots, Q_{\mathrm{2}}}_{p} \}, \\
	&\Phi = \left[
	\begin{array}{ccccc}
		I \\
		A_{k} \\
		A_{k}^{2} \\
		\vdots \\
		A_{k}^{p}
	\end{array}
	\right], F = \left[
	\begin{array}{cccc}
		0 \\
		B_{k} & 0 \\
		A_{k}B_{k}& B_{k} & \ddots \\
		\vdots & \cdots & \ddots & 0 \\
		A_{k}^{p-1}B_{k} & \cdots & A_{k}B_{k} & B_{k}
	\end{array}
	\right]
\end{align*}
and
\begin{align*}
	U(k) &= \left[ u^{\rm T}(k),u^{\rm T}(k+1|k),\cdots, u^{\rm T}(k+p-1|k) \right]^{\rm T}, \\ \nonumber
	\hat U(k) &= \left[ \hat u^{\rm T}(k),\hat u^{\rm T}(k+1|k),\cdots, \hat u^{\rm T}(k+p-1|k) \right]^{\rm T}, \\ \nonumber
	u^{\rm T}(k)&=u^{\rm T}(k|k),\hat u^{\rm T}(k)=\hat u^{\rm T}(k|k)
\end{align*}

With the help of Definition \ref{Def1}, the matrix $W$ and the parameter $\mu$ will be designed using the following theorem.

\begin{thm}\label{thm1}
	Consider Problem \ref{Pro1}, if $\exists \bar{W}^{\rm T}=\bar{W}>0,~S \in \mathbb{R}^{2\times3},~\bar{S}=S\bar{W}$ and an unknown scalar $\tau$, make the following linear matrix inequality:
	\begin{align}\label{eq19}
		&\left[
		\begin{array}{cccccc}
			-\bar{W} & * & * & * & * & *\\
			\Pi_{1} & -\Pi_{2} & * & * & * & *\\
			B_{k}\bar{S} & \bar{\boldsymbol{\gamma}}B_{k}B_{k}^{\rm T} & -\Pi_{3} & * & * & *\\
			\bar{S} & 0 & 0 & -Q^{-1}_{2} & * & *\\
			\bar{W} & 0 & 0 & 0 & -Q^{-1}_{1} & *\\
			\bar{S} & 0 & 0 & 0 & 0 & - \tau^{-1} I
		\end{array}
		\right] \leq 0
	\end{align}
	hold, in that way the closed-loop system has mean square asymptotic stability, where $\Pi_{1} = A_{k}\bar{W}+\bar{\boldsymbol{\gamma}}B_{k}\bar{S},~\Pi_{2} = \bar{W}-\bar{\boldsymbol{\gamma}}^{2}B_{k}B_{k}^{\rm T},~\Pi_{3} = \hat{\gamma}\bar{W}-B_{k}B_{k}^{\rm T},~\hat{\gamma} = \dfrac{1}{\bar\gamma(1-\bar\gamma)}$. The terminal weighting matrix $W = \bar{W}^{-1}$ in $J(k)$.
\end{thm}

The formation of the LMI in Theorem \ref{thm1} and the stability of the closed-loop system will be explained below.

\begin{pf}
	For uncertain discrete system \eqref{eq9} and performance metric \eqref{eq10}, let the state feedback control law $u(k)=Sx(k)$ and the uncertain term $B_{k}\hat{u}(k)=\Delta B_{k}u(k)$, where $\Delta B_{k}$ represents the system indeterminate parameter. The Lyapunov function is defined as:
	\begin{align*}
		V(x(k)) = x^{\rm T}(k)Wx(k).
	\end{align*}
	The closed-loop system is robust and stable if
	\begin{align*}
		V(x(k+1))-V(x(k)) \leq -\left[ \left\| x(k)\right\|_{Q_{\mathrm{1}}}^{2}+\left\| u(k)\right\|_{Q_{\mathrm{2}}}^{2}\right] 
	\end{align*}
	holds. By deformation, the above equation can be equivalent to:
	\begin{align*}
		\left\| A_{k}+\gamma (k)\left( B_{k} + \Delta B_{k} \right)S \right\|_{W}^{2} - W +Q_{\mathrm{1}}+\left\| S\right\|_{Q_{\mathrm{2}}}^{2} \leq 0
	\end{align*}
	Consider the existence of a random variable, $\gamma (k)$, and take the expectation for inequality. In addition, there is $B_{k}\hat{S}=\Delta B_{k}S$. Thus, the above result can be further rewritten to the form of \eqref{eq13}.
	
	About the dynamic quantizer $q_{u}$ we know: $q_{u}(S)=S+\hat{S}$, and the inequality \eqref{eq13} in Definition \ref{Def1} can be rewritten as:
	\begin{align*}
		E\left\lbrace \left\| A_{k}+\gamma (k)B_{k}\left( S + \hat{S} \right) \right\|_{W}^{2} \right\rbrace+Q_{\mathrm{1}}+\left\| S\right\|_{Q_{\mathrm{2}}}^{2} - W \leq 0
	\end{align*}
	where $\hat S = q_{u}(S)-S$. The following inequality can be obtained by using the same method as:
	\begin{align}\label{eq20}
		0 \geq &\left\| S \right\|_{Q_{\mathrm{2}}}^{2} + \bar{\boldsymbol{\gamma}}\left( 1-\bar{\boldsymbol{\gamma}} \right)\left\|B_{k}\left( S + \hat{S} \right) \right\|_{W}^{2} + Q_{\mathrm{1}} \\ \nonumber
		& +\left\| A_{k} + \bar{\boldsymbol{\gamma}}B_{k}\left( S + \hat{S} \right) \right\|_{W}^{2} - W
	\end{align}
	According to the Schur complement, the above inequality can be transformed into the following form of matrix inequality \cite{xue2010robust}:
	\begin{align*}
		&\left[
		\begin{array}{ccccc}
			-\bar{W}^{-1} & * & * & * & *\\
			A_{k}+\bar{\boldsymbol{\gamma}}B_{k}S+\bar{\boldsymbol{\gamma}}B_{k}\hat{S} & -\bar{W} & * & * & *\\
			B_{k}S+B_{k}\hat{S} & 0 & -\hat{\gamma}\bar{W} & * & *\\
			S & 0 & 0 & -Q^{-1}_{2} & *\\
			I & 0 & 0 & 0 & -Q^{-1}_{1}
		\end{array}
		\right] \leq 0 \\
		&\Longrightarrow
		\left[
		\begin{array}{ccccc}
			-\bar{W}^{-1} & * & * & * & *\\
			A_{k}+\bar{\boldsymbol{\gamma}}B_{k}S & -\bar{W} & * & *\\
			B_{k}S & 0 & -\hat{\gamma}\bar{W} & * & *\\
			S & 0 & 0 & -Q^{-1}_{2} & *\\
			I & 0 & 0 & 0 & -Q^{-1}_{1}
		\end{array}
		\right] + \left[
		\begin{array}{c}
			\hat S^{\rm T} \\
			0 \\
			0 \\
			0 \\
			0
		\end{array}
		\right]
	\end{align*}
	\begin{align*}
		\left[ 
		\begin{array}{c}
			0~~\bar{\boldsymbol{\gamma}}B_{k}^{\rm T}~~B_{k}^{\rm T}~~0~~0
		\end{array}
		\right] + \left[
		\begin{array}{c}
			0 \\
			\bar{\boldsymbol{\gamma}}B_{k} \\
			B_{k} \\
			0 \\
			0
		\end{array}
		\right] \left[ 
		\begin{array}{c}
			\hat S~~0~~0~~0~~0
		\end{array}
		\right] \leq 0
	\end{align*}
	Here we use the deformation of Lemma \ref{lem1} and come to the conclusion that \eqref{eq20} is sure to hold if:
	\begin{align*}
		&\left[
		\begin{array}{ccccc}
			-\bar{W}^{-1} & * & * & * & *\\
			A_{k}+\bar{\boldsymbol{\gamma}}B_{k}S & -\bar{W} & * & *\\
			B_{k}S & 0 & -\hat{\gamma}\bar{W} & * & *\\
			S & 0 & 0 & -Q^{-1}_{2} & *\\
			I & 0 & 0 & 0 & -Q^{-1}_{1}
		\end{array}
		\right] + \left[
		\begin{array}{c}
			0 \\
			\bar{\boldsymbol{\gamma}}B_{k} \\
			B_{k}\\
			0 \\
			0
		\end{array}
		\right] M  \\
		&\left[ 
		\begin{array}{c}
			0~~\bar{\boldsymbol{\gamma}}B_{k}^{\rm T}~~B_{k}^{\rm T}~~0~~0
		\end{array}
		\right] + \left[
		\begin{array}{c}
			\hat S^{\rm T} \\
			0 \\
			0 \\
			0 \\
			0
		\end{array}
		\right] M^{-1} \left[ 
		\begin{array}{c}
			\hat S~~0~~0~~0~~0
		\end{array}
		\right] \leq 0
	\end{align*}
	where the matrix $M$ is a positive definite matrix. According to \eqref{eq17}, there is:
	\begin{align*}
		\hat S^{\rm T}\hat S \leq \dfrac{\Delta_{u}^{\rm 2} }{\alpha_{u}M_{u}^{2}} S^{\rm T}S.
	\end{align*}
	As we all know, the identity matrix $I$ also belongs to the positive definite matrix. For convenience, set $M=I$. At the same time, merge the first and second items, and then substitute the result of the above derivation into the inequality. If the following inequality
	\begin{align*}
		&\left[
		\begin{array}{ccccc}
			-\bar{W}^{-1} & * & * & * & *\\
			A_{k}+\bar{\boldsymbol{\gamma}}B_{k}S & -\Pi_{2} & * & * & *\\
			B_{k}S & \bar{\boldsymbol{\gamma}}B_{k}B_{k}^{\rm T} & -\Pi_{3} & * & *\\
			S & 0 & 0 & -Q^{-1}_{2} & *\\
			I & 0 & 0 & 0 & -Q^{-1}_{1} 
		\end{array}
		\right] \\ \nonumber
		&+ \dfrac{\Delta_{u}^{\rm 2} }{\alpha_{u}M_{u}^{2}} S^{\rm T}S \leq 0
	\end{align*}
	holds, then it can be concluded that \eqref{eq20} must hold. Then use the Shure complement once again and define $\dfrac{\Delta_{u}^{\rm2}}{\alpha_{u}M_{u}^{2}} = \tau$. Finally, multiplying both sides of the above inequality by ${\rm diag} \{ \bar{W}, I, I, I, I, I \}$ on the left and right simultaneously yields \eqref{eq19}. By the way, since the two coefficient matrices $A_{k},~B_{k}$ are updated in real time due to the desired heading angle $\varphi_d$, it is also necessary to solve the above matrix inequalities in real time.
	
	At time $k$, the optimal predictive cost function can be expressed as:
	\begin{align}\label{eq21}
		J^{*}(k)=\sum_{i=0}^{p-1} &\left[ \left\| x(k+i|k)\right\|_{Q_{\mathrm{1}}}^{2} + \left\| u^{*}(k+i|k)\right\|_{Q_{\mathrm{2}}}^{2} \right] \\ \nonumber
		&+\left\| x(k+p|k)\right\|_{W}^{2}
	\end{align}
	The following contractile derivations are made according to \eqref{eq21}:
	\begin{align}\label{eq22}
		&E\big\{ J^{*}(k+1)\big\} - E\left\lbrace J^{*}(k)\right\rbrace= \nonumber\\
		&E\Bigg\{ \sum_{i=0}^{p-1} \left[ \left\| u^{*}(k+1+i|k+1) \right\|_{Q_{\mathrm{2}}}^{2} - \left\| u^{*}(k+i|k)\right\|_{Q_{\mathrm{2}}}^{2}\right.   \nonumber\\
		&+\left. \left\| x(k+1+i|k+1)\right\|_{Q_{\mathrm{1}}}^{2} - \left\| x(k+i|k)\right\|_{Q_{\mathrm{1}}}^{2} \right] \\ \nonumber
		&+\left\| x(k+1+p|k+1) \right\|_{W}^{2} - \left\| x(k+p|k) \right\|_{W}^{2} \Bigg\}
	\end{align}
	Due to the trajectory tracking problem discussed in this article, theoretically, the state error at each moment should be smaller than the previous moment, so the final state error should approach zero infinitely. According to the theories of predictive control, the following equality holds:
	\begin{align}\label{eq23}
		&x(k+p|k+1) = x^{*}(k+p|k), \\ \nonumber
		&u(k+i|k+1) = u^{*}(k+i|k)
	\end{align}
	where $i=0,1\cdots,p-1$. Then, combined with \eqref{eq23}, conduct an in-depth derivation on \eqref{eq22}:
	\begin{align}\label{eq24}
		&E\big\{ J^{*}(k+1)\big\} - E\left\lbrace J^{*}(k)\right\rbrace \\ \nonumber
		\leq& E\Bigg\{ \sum_{i=0}^{p-1} \left[ \left\| u^{*}(k+1+i|k) \right\|_{Q_{\mathrm{2}}}^{2} - \left\| u^{*}(k+i|k)\right\|_{Q_{\mathrm{2}}}^{2}  \right. \Big. \\ \nonumber
		&+\left.\left\| x(k+1+i|k)\right\|_{Q_{\mathrm{1}}}^{2} - \left\| x(k+i|k)\right\|_{Q_{\mathrm{1}}}^{2} \right] \\ \nonumber
		&+\Big.  \left\| x(k+1+p|k) \right\|_{W}^{2} - \left\| x(k+p|k) \right\|_{W}^{2} \Bigg\} \\ \nonumber
		\leq& E\Big\{ \left\| u^{*}(k+p|k) \right\|_{Q_{\mathrm{2}}}^{2} - \left\| u^{*}(k)\right\|_{Q_{\mathrm{2}}}^{2} + \left\| x(k+p|k)\right\|_{Q_{\mathrm{1}}}^{2} \Big.  \\ \nonumber
		&-\Big. \left\| x(k)\right\|_{Q_{\mathrm{1}}}^{2} \Big\} \\ \nonumber
		\leq& 0.
	\end{align}
	At the same time, there is also:
	\begin{align}\label{eq25}
		J^{*}(k)=&\sum_{i=1}^{p-1} \left[ \left\| x(k+i|k)\right\|_{Q_{\mathrm{1}}}^{2} +\left\| u^{*}(k+i|k)\right\|_{Q_{\mathrm{2}}}^{2} \right]  \nonumber\\
		&+\left\| x(k+p|k)\right\|_{W}^{2} + \left\| x(k)\right\|_{Q_{\mathrm{1}}}^{2} \\ \nonumber
		\Longrightarrow& J^{*}(k)\geq \left\| x(k)\right\|_{Q_{\mathrm{1}}}^{2}\geq \left\| x(k)\right\|
	\end{align}
	When $\left\| x(k)\right\|\rightarrow\infty$, one has $J^{*}(k)\rightarrow\infty$. According to Lyapunov's stability theorem, it can be concluded that the closed-loop system has global asymptotic stability in the mean square sense.
\end{pf}

Let $Y(k) = F \hat U(k)$, then the following inequality holds:
\begin{align}\label{eq26}
	\left\| Y(k)\right\| \leq \dfrac{ \Delta_{u}}{\sqrt{\alpha_{u}}M_{u} }\left\| F\right\| \left\| U(k)\right\| .
\end{align}
Assuming $f = \dfrac{ \Delta_{u}}{\sqrt{\alpha_{u}}M_{u} }\left\| F\right\|$, the above derivation can be summarized as follows: $\left\| Y(k)\right\| \leq f\left\| U(k) \right\|$. Then, optimization problem \ref{Pro1} can be transformed into:
\begin{align}\label{eq27}
	&U^{*}(k)= \mathop{arg} \ \mathop{min}\limits_{U(k)} \mathop{max}\limits_{\left\| Y(k)\right\| \leq f\left\| U(k) \right\|} E\left\lbrace J(k)\right\rbrace \nonumber\\
	&s.t.~\eqref{eq12}~and~\eqref{eq18}.
\end{align}
If the above optimization problem has a unique solution, then the global optimal solution $u^{*}(k)$ of Problem \ref{Pro1} is equal to the first term of the optimal control sequence $U^{*}(k)$, i.e., $u^{*}(k)=[I,0,\cdots,0]U^{*}(k)$.

\begin{pf}
	Given the probability of successful data transmission $\bar{\boldsymbol{\gamma}}$, the dynamic quantization range $M_{u}$, and the error of the quantizer $\Delta_{u}$ and the reference trajectory, the optimal control input of system \eqref{eq9} can be obtained by solving the optimization problem as follows:
	\begin{align}
		E\left\lbrace J(k)\right\rbrace
		=&E\left\lbrace \left\| \Phi x(k) + \gamma (k)FU(k) + \gamma (k)Y(k) \right\|_{Q_{w}}^{2}\right\rbrace \nonumber \\
		&+ \left\| U(k) \right\|_{Q}^{2}
	\end{align}\label{eq28}
	
	\begin{remark}
		There is an uncertainty term $\gamma (k)$ in the nominal data of the cost function, while the optimization problem in \eqref{eq27} is a least squares problem. As we all know, if the optimization problem is solved directly at this time, the robustness of the cost function may be weakened, and the final solution will be de-regularization. So we first deal with the random probability by referring to the method mentioned above: $E\left\lbrace \gamma(k) \right\rbrace = \bar{\boldsymbol{\gamma}}$, and then derive a robust regularization least squares solution.
	\end{remark}
	
	Further simplification of equality \eqref{eq28} yields:
	\begin{align}\label{eq29}
		E\left\lbrace J(k)\right\rbrace
		=&\left\| \Phi x(k) + \bar{\boldsymbol{\gamma}}\left( FU(k) + Y(k) \right) \right\|_{Q_{w}}^{2}+  \\ \nonumber
		&\left\| U(k) \right\|_{Q}^{2} + \bar{\boldsymbol{\gamma}}\left( 1-\bar{\boldsymbol{\gamma}} \right) \left\| FU(k) + Y(k) \right\|_{Q_{w}}^{2}
	\end{align}
	To facilitate the determination of whether the optimization problem has solutions, we construct a composite function here: $C\left( U,Y \right) = || \bar{\boldsymbol{\gamma}}\left( FU(k) + Y(k) \right) + \Phi x(k)|| _{Q_{w}}^{\rm2} + \bar{\boldsymbol{\gamma}}\left( 1-\bar{\boldsymbol{\gamma}} \right) \left\| FU(k) + Y(k) \right\|_{Q_{w}}^{\rm2}$. Therefore, based on the calculation results of \eqref{eq29}, the optimization problem in \eqref{eq27} can be rewritten as:
	\begin{align}\label{eq30}
		U^{*}(k)=&\mathop{arg} \ \mathop{min}\limits_{U(k)} \mathop{max}\limits_{\left\| Y(k)\right\| \leq f\left\| U(k) \right\|} \left[ C\left( U,Y \right) \right.  \\ \nonumber
		&+ \left. \left\| U(k) \right\|_{Q}^{2} \right]
	\end{align}
	
	For the optimization problem in \eqref{eq30}, we first analyze its function structure. When considering $Y(k)$ as a fixed constant and observing the maximization problem separately:
	\begin{align}\label{eq31}
		\mathop{max}\limits_{\left\| Y(k)\right\| \leq f\left\| U(k) \right\|} C\left( U,Y \right).
	\end{align}
	It is not difficult to find that $\left( C\left( U,Y \right) + \left\| U(k) \right\|_{Q}^{2}\right) $ about $U(k)$ is convex. Therefore, we can conclude that Problem \ref{Pro1} has a unique global optimal solution $U^{*}(k)$.
	Secondly, considering $U(k)$ as a constant, the maximizing problem in \eqref{eq31} is also convex about $Y(k)$. So the maximization problem in \eqref{eq31} will be solved at the boundary: $\left\| Y(k)\right\| = f\left\| U(k) \right\|$. This allows the inequality constraint to be transformed into an equation constraint.
	
	\begin{remark}
		For optimization problems with constraints, most previous studies have mostly used the LMI technique to deal with inequality constraints \cite{kothare1996robust}. For the equality constraint, the LM is almost always used to derive the optimal solution \cite{sayed2002regularized}. Here we have converted the inequality constraint into an equation constraint, which saves a lot of time for solving. In order to determine the feasible solution more quickly, we use the ALM to solve the optimization problem.
	\end{remark}
	Using the ALM, \eqref{eq31} is constructed as an unconstrained maximum problem:
	\begin{align}\label{eq32}
		C^{*}\left( U \right) =\mathop{max}\limits_{ Y(k),\lambda }&\left( \frac{\sigma}{2} \left(  \left\| Y(k)\right\|^{2} - f^{2}\left\| U(k) \right\|^{2}\right)^{2}+\right.  \\ \nonumber
		&\left. C\left( U,Y \right) - \lambda\left[ \left\| Y(k)\right\|^{2} - f^{2}\left\| U(k) \right\|^{2} \right]\right)
	\end{align}
	We use $C_{\sigma}(Y,\lambda)$ to represent its cost. After calculation, obtain the optimal solutions $Y^{*}$ and $\lambda^{*}$:
	\begin{align}
		&\left( \lambda^{*}I - \bar{\boldsymbol{\gamma}}Q_{w} \right)Y^{*} = \bar{\boldsymbol{\gamma}}Q_{w}\left[ \Phi x(k) + FU(k) \right], \\ \nonumber
		&\left\| Y(k)\right\|^{2} = f^{2}\left\| U(k) \right\|^{2}
	\end{align}\label{eq33}
	During the process of solving $\lambda^{*}$, we found that the Hessian matrix for $Y(k)$ in equation \eqref{eq33} is equal to: $\bar{\boldsymbol{\gamma}}Q_{w} - \lambda I$. In order to ensure the feasibility of the solution, the Hessian matrix must be semi-negative definite. Hence, the prerequisite for $\lambda = \lambda^{*}$ is $\lambda^{*}I \geq \bar{\boldsymbol{\gamma}}Q_{w}$. By using the derivative results in \eqref{eq33}, we have:
	\begin{align}\label{eq34}
		c\left( Y \right) = \left\| Y\right\|^{2} - f^{2}\left\| U \right\|^{2}
	\end{align}
	where the symbol $\left( \lambda^{*}I - \bar{\boldsymbol{\gamma}}Q_{w} \right)^{\dagger}$ represents the pseudo-inverse of $\lambda^{*}I - \bar{\boldsymbol{\gamma}}Q_{w}$. At this point, the maximization problem has been solved, and the next step is to solve the minimization problem. For the convenience of calculation, let $\omega(\lambda) = Q_{w} + \bar{\boldsymbol{\gamma}}Q_{w}\left( \lambda I - \bar{\boldsymbol{\gamma}}Q_{w} \right)^{\dagger} \bar{\boldsymbol{\gamma}}Q_{w}$. So the problem in \eqref{eq30} is equivalent to the following minimization problem:
	\begin{align}\label{eq35}
		\mathop{min}\limits_{U(k)}\left[ \bar{C}(U,\lambda) + \left\| U(k) \right\|_{Q}^{2} \right] 
	\end{align}
	where
	\begin{align*}
		\bar{C}(U,\lambda) =& \left\| \Phi x(k) + FU(k)\right\|_{ \omega(\lambda) }^{2} + \left( 1-\bar{\boldsymbol{\gamma}} \right)\left\| \Phi x(k) \right\|_{ Q_{w} }^{2} + \\
		&\left[ \lambda-\sigma_{k} c(Y_{k})\right] f^{2}\left\| U(k) \right\|^{2}
	\end{align*}
	Based on the sufficient condition: $\lambda \geq \bar{\boldsymbol{\gamma}}\left\| Q_{w} \right\|$. Subsequently,  the results derived from \cite{sayed2002regularized} are used. We ultimately obtain the following theorem.
	
	\begin{thm}\label{thm2}
		Consider a trajectory tracking system \eqref{eq9} with uncertain parameters. With the probability $\bar{\boldsymbol{\gamma}}$, penalty factor $\sigma_{0}$, the dynamic quantization range $M_{u}$, and the error of the quantizer $\Delta_{u}$, the optimal control sequence is obtained by solving the regularization least squares problem in \eqref{eq30}:
		\begin{align}\label{eq36}
			&\left( Q + \left[ \lambda^{*}-\sigma_{k} c(Y_{k})\right] f^{2}I + F^{\rm T} \omega(\lambda^{*})F \right) U^{*}(k) \\ \nonumber
			&=-F^{\rm T}\omega(\lambda^{*}) \Phi x(k)
		\end{align}
		the weight matrices $Q>0,~Q_{w}\geq0$. When substituting \eqref{eq36} into \eqref{eq35}, the parameter $\lambda$ is obtained by the following minimization problem:
		\begin{align}\label{eq37}
			\lambda^{*} = \mathop{arg} \ \mathop{min}\limits_{\lambda \geq \bar{\boldsymbol{\gamma}}\left\| Q_{w} \right\|} L(\lambda)
		\end{align}
		where $L(\lambda) =f^{2}\left[ \lambda-\sigma_{k} c(Y_{k})\right] \left\| U^{*}(k) \right\|^{2} + \left\| U^{*}(k) \right\|_{Q}^{2} +  \left\| \Phi x(k) + FU^{*}(k)\right\|_{ \omega(\lambda) }^{2}$.
	\end{thm}
	
	It is easy to find that $\left( 1-\bar{\boldsymbol{\gamma}} \right)\left\| \Phi x(k) \right\|_{ Q_{w} }^{2}$ is completely independent and not affected by decision variables $U^{*}(k)$ and $\lambda$. Therefore, we construct a new cost function $L(\lambda)$. When $\lambda$ is within $[\bar{\boldsymbol{\gamma}}\left\| Q_{w} \right\|,+\infty)$, the cost $L(\lambda)$ is unimodal with respect to $\lambda$, meaning \eqref{eq37} has a unique global minimum.
\end{pf}

From the above proof results, it can be seen that the optimal solution of the system depend on the unknown terminal weighting matrix $W$ in the cost function $J(k)$ and the dynamic parameter $\mu$.

To guarantee the effectiveness of the optimization problem in \eqref{eq27}, the smallest $L(\lambda)$ must be found to determine $\lambda^{*}$. By solving the LMI in Theorem \ref{thm1}, the dynamic parameter of the dynamic quantizer $\mu$ and the matrix variable $W$ can be calculated. Then, determine the lower limit $\bar{\boldsymbol{\gamma}}\left\| Q_{w}\right\|$ of the variable $\lambda$. Finally, the values of $L(\lambda)$ are compared using an iterative approach. The overall iteration process of the algorithm is shown in Algorithm \ref{alg1}.

\makeatletter
\def\BState{\State\hskip-\ALG@thistlm}
\makeatother
\begin{algorithm}
	\caption{Iteration of robust predictive control with generalized multiplier.}\label{alg1}
	\begin{algorithmic}[1]
		\State At $k$=$0$, set the following parameter values: $v_{d}$, $\varphi_{d}$, $\delta_{d}$, $l$, $r$, $Q_{\mathrm{1}}$, $Q_{\mathrm{2}}$, $T$, $\Delta_{u}$, $M_{u}$, $\bar{\boldsymbol{\gamma}}$, $\sigma_{0}$, $\rho$, $\varepsilon$, $\eta$, and $x(0)$.
		\Procedure{Iteration}{}
		\State Make use of Matlab to solve inequality in \eqref{eq19} under the condition of Theorem \ref{thm1}. When it has feasible solutions, we obtain ${W}$ and $\alpha_{u}$. Use the terminal weight matrix ${W}$ to calculate $\bar{\boldsymbol{\gamma}}\left\| Q_{w}\right\|$ as $\lambda_{1}$. Then select the prediction range $N$, the maximum number of cycles $H$.
		Set $h = 1$, $\lambda_{1}\left( \lambda_{h}\right)$ as initialization $\lambda$.
		\BState \emph{top}:
		\If {$h < H$}
		\State Substitute $\lambda_{h}$ into \eqref{eq26} to calculate $U_{h}$. Next, according to the equation for $L(\lambda)$ in optimization problem in \eqref{eq37}, use the obtained $\lambda_{h}$ and $U_{h}$ to compute $L_{h}$.
		\State \textbf{goto} \emph{loop}.
		\EndIf
		\State \textbf{close};
		\BState \emph{loop}:
		\If {$L_{h-1}-L_{h}<\eta$}
		\State It is proven that the problem in \eqref{eq37} has an optimal solution, which is $\lambda^{*}=\lambda_{h}$.
		\State \textbf{close};
		\EndIf
		\State update the number of iterations: $h$=$h$+$1$. Calculate $Y_{h}$ according to \eqref{eq33}, simultaneously update the numeric value $\lambda_{h}$=$\lambda_{h}$+$\sigma_{h-1}c(Y_{h})$ and the penalty factor $\sigma_{h}=\rho\sigma_{h-1}$.
		\State \textbf{goto} \emph{top}.
		\EndProcedure
		\State We default to $U_{h}$=$U^{*}$, cut off the first component of the control sequence $U_{h}$ as the input $u^{*}$ and substitute it into state equation \eqref{eq4} in the form of feedback to update error $x(k+1)$. Finally, make $k=k+1$ and continue iteration.
	\end{algorithmic}
\end{algorithm}

\begin{figure*}[htbp]
	\centering
	\begin{minipage}[t]{0.48\textwidth}
		\centering
		\subfigure[LM]{
			\includegraphics[width=0.8\linewidth]{lm}}
		\label{fig.LM}
	\end{minipage}
	\begin{minipage}[t]{0.48\textwidth}
		\centering
		\subfigure[ALM]{
			\includegraphics[width=0.8\linewidth]{alm}}
		\label{fig.ALM}
	\end{minipage}
	\caption{When penalty factor $\sigma=2$, the change of Contour line corresponding to the two methods.}
	\label{1}
\end{figure*}

\begin{remark}
	It is worth noting that the trajectory tracking problem needs to ensure real-time performance in the actual process. Using iterative methods to determine the optimal control input $u^{*}$ will inevitably increase the computational complexity. Because during iterative calculation, solving a feasible solution to a matrix inequality will take a lot of time, the size of the prediction range $p$ also affects the amount of calculation. Hence, in order to ensure real-time calculation while improving calculation accuracy, we need to balance calculation time and prediction range. In many motion control systems, the prediction range often adopts the one-step control form of $p=1$ to update the state faster.
\end{remark}

\section{Simulation example}

Before trajectory simulation, we use a simple numerical example to compare the LM with the ALM. Firstly, establish an optimization problem:
\begin{align*}
	\mathop{min}&~~~x+\sqrt{3}y, \\
	s.t.&~~~x^{2}+y^{2}=1
\end{align*}
It is easy to find the optimal solution to this problem as $x^{*}=\left( -\frac{1}{2},-\frac{\sqrt{3}}{2}\right)^{\rm T}$. The analytical expression calculated by LM is: $L(x,y,\lambda)=x+\sqrt{3}y+\lambda\left( x^{2}+y^{2}-1\right)$. Based on ALM, we consider the augmented Lagrange function: $L_{\sigma}(x,y,\lambda)=x+\sqrt{3}y+\lambda\left( x^{2}+y^{2}-1\right) +\frac{\sigma}{2}\left( x^{2}+y^{2}-1\right) ^{2}$.

The Contour line of $L_{\rm 2}(x, y, 0.9)$ is drawn in Figure \ref{1}. The points marked with ``$*$" in the figure are the optimal solutions of the original problem $x^{*}$, and the points marked with ``$\circ$" are the optimal solutions of LM or ALM. It can be seen that the optimal solution of the ALM is closer to the real solution.

\subsubsection{Example 1.}\label{ex1}
	In this simulation, we assume that there is a circular trajectory in virtual space. The radius of this ideal track is $0.8$m, and the orbital center coordinate is $O(0,0,0)$. The ideal linear speed that the autonomous vehicle needs to reach is $0.4~\rm m/s$, and the ideal angular speed is $0.5$ rad/s. For circular trajectories, the desired front wheel deflection angle $\delta_{d}$ is fixed and unchanged. This is because the expected angular velocity is a constant value, and when the ideal state is reached, the front wheels uniformly deflect during movement. And from the geometric relationship, we can also draw this conclusion: $\tan \delta_d = l/r$. The wheelbase of the front and rear wheels $l$ and radius $r$ are predetermined values, in addition, control the sampling time to $T=0.1\rm s$. So we obtained the coefficient matrix of the error model:
	\begin{align*}
		A_{k}&=\left[
		\begin{array}{ccc}
			1 & 0 & -0.04\sin \varphi_d \\
			0 & 1 & 0.04\cos \varphi_d \\
			0 & 0 & 1 \\
		\end{array}
		\right], B_{k}=\left[
		\begin{array}{cc}
			0.1\cos \varphi_d & 0 \\
			0.1\sin \varphi_d & 0 \\
			0.125 & 0.2125 \\
		\end{array}
		\right]
	\end{align*}
	Similarly, the desired heading angle $\varphi_d$ also varies uniformly in real-time. When initializing, we select the zero-time error $x(0)=[-0.1~-0.2~0.2]^{\rm T}$ for the car, initial position $x_{e}(0)=[0.6~0~0]^{\rm T}$, weight matrices $Q_{\mathrm{1}}=2I,~ Q_{\mathrm{2}}=0.05I$, axle length $l=0.2~\rm m$, the quantization error $\Delta_{u}=1$, the range of quantizer $M_{u}=10$, the Bernoulli probability $\bar{\boldsymbol{\gamma}}=0.9$, the prediction range $p=2$, the cycle iteration limit number $H=100$, the iteration step size $\hat{\lambda}=1$, and the coefficient of determination $\eta=0.01$.
	
	In addition, according to \eqref{eq16}, to demonstrate the effect of the dynamic quantizer during simulation, assume $\hat u(k) = \dfrac{\Delta_{u}\cos(k)}{\sqrt{\alpha_{u}}M_{u}} u(k)$. Therefore, the system equation can be rewritten in the form with coefficients:
	\begin{align*}
		x(k+1)=&\left[
		\begin{array}{ccc}
			1 & 0 & -0.04\sin \varphi_d \\
			0 & 1 & 0.04\cos \varphi_d \\
			0 & 0 & 1 \\
		\end{array}
		\right]x(k)+ \\
		&0.9\left( 1+\dfrac{\Delta_{u}\cos(k)}{\sqrt{\alpha_{u}}M_{u}}\right)\left[
		\begin{array}{cc}
			0.1\cos \varphi_d & 0 \\
			0.1\sin \varphi_d & 0 \\
			0.125 & 0.2125 \\
		\end{array}
		\right]u(k)
	\end{align*}
	
	\begin{figure*}[htbp]
		\centering
		\begin{minipage}[t]{0.48\textwidth}
			\centering
			\subfigure[]{
				\label{fig.2}
				\includegraphics[width=0.8\linewidth]{trajectoryO}}
		\end{minipage}
		\begin{minipage}[t]{0.48\textwidth}
			\centering
			\subfigure[]{
				\label{fig.3}
				\includegraphics[width=0.8\linewidth]{xe-xd,ye-yd,phie-phidO}}
		\end{minipage}
		\begin{minipage}[t]{0.48\textwidth}
			\centering
			\subfigure[]{
				\label{fig.4}
				\includegraphics[width=0.8\linewidth]{wvO}}
		\end{minipage}
		\begin{minipage}[t]{0.48\textwidth}
			\centering
			\subfigure[]{
				\label{fig.5}
				\includegraphics[width=0.8\linewidth]{jO}}
		\end{minipage}
		\caption{The circular trajectory simulation results.}
		\label{2}
	\end{figure*}
	
	In the MATLAB simulation environment, based on the designed robust predictive control algorithm, we first substitute the corresponding heading angle $\varphi_d$ at that time into inequality in \eqref{eq31} to calculate $W$ and $\alpha_{u}$. Then use $W$ and $Q_{\mathrm{1}}$ to calculate $\bar{\boldsymbol{\gamma}}\left\| Q_{w}\right\|$, and use the forward method to obtain $\lambda_{1}$. The maximum number of iterations for a prediction horizon of 2 is calculated as $h_{\rm max}=15$ according to the designed algorithm. At each moment, the minimum number of iterations required for the LMI toolbox to calculate a feasible solution is between 21 and 46. The trajectory simulation result is shown in Figure \ref{fig.2}, where the solid line represents the desired trajectory and the dashed line represents the actual trajectory of the autonomous vehicle. The circle represents the starting point of the vehicle. It can be seen that the vehicle successfully tracks the planned trajectory and operates smoothly after passing through approximately one-third of the circle. The state curve of the entire error system is shown in Figure \ref{fig.3}, where the red curve represents the error in the $X$-axis direction, the black curve represents the error in the $Y$-axis direction, and the blue curve represents the deviation of the heading angle. Figure \ref{fig.4} shows the function curve of the actual linear velocity and angular velocity, which are calculated based on the system error $x$, the control input $u$, and the corresponding reference trajectory at each time. It can be observed from Figure \ref{fig.3} and Figure \ref{fig.4} that the autonomous vehicle can reach the ideal state after approximately 3 seconds, the system error tends to zero, and the speed matches the reference speed set. Computing the cost function function from \eqref{eq18}, and the convergence curve of the entire system at each moment is shown in Figure \ref{fig.5}.
	
	Based on this, we can conclude that the robust predictive control algorithm designed for autonomous vehicles can effectively solve the trajectory tracking problem.

\begin{figure*}[htbp]
	\centering
	\begin{minipage}[t]{0.48\textwidth}
		\centering
		\subfigure[]{
			\label{fig.6}
			\includegraphics[width=0.8\linewidth]{track8}}
	\end{minipage}
	\begin{minipage}[t]{0.48\textwidth}
		\centering
		\subfigure[]{
			\label{fig.7}
			\includegraphics[width=0.8\linewidth]{x8}}
	\end{minipage}
	\begin{minipage}[t]{0.48\textwidth}
		\centering
		\subfigure[]{
			\label{fig.8}
			\includegraphics[width=0.8\linewidth]{we,ve8}}
	\end{minipage}
	\begin{minipage}[t]{0.48\textwidth}
		\centering
		\subfigure[]{
			\label{fig.9}
			\includegraphics[width=0.8\linewidth]{J8}}
	\end{minipage}
	\caption{The eight-shaped trajectory simulation results.}
	\label{3}
\end{figure*}

\subsubsection{Example 2.}\label{ex2}
	This time, we have selected an eight-shaped trajectory as the reference path. Its equation with respect to time is $x_d = \sin(0.1t),y_d=\sin(0.2t)$. Compared to a circular trajectory, both the variation of angles and the control of position are more complex. Firstly, we assume that the ideal linear velocity of the autonomous vehicle is $0.2~\rm m/s$. Due to the non-uniformity of angle variation, expressions for both the heading angle of the vehicle and the steering angle of the front wheel need to be calculated in advance. By taking the derivatives of $x_d$ and $y_d$ with respect to $t$ and then taking their ratio, we obtain the equation for the slope of the curve, which corresponds to the heading angle of the vehicle. Finally, we can calculate the steering angle of the front wheel and the expected angular velocity at the corresponding time. The relevant functions calculated are as follows:
	\begin{align*}
		\dfrac{\mathrm{d} y_d}{\mathrm{d} x_d}=&\dfrac{2\cos(0.2t)}{\cos(0.1t)}=\tan\varphi_d \\
		\omega_d=&\dfrac{\mathrm{d} \varphi_d}{\mathrm{d} t}=0.1\tan(0.1t)\sin\varphi_d\cos\varphi_d \\
		&-0.8\sin(0.1t)\cos^{2}\varphi_d \\
		\delta_i=&\arctan(\dfrac{l\omega_i}{v_i}),~i=d,e.
	\end{align*}
	
	Similarly, we set an initial error state and starting point for the vehicle: $x(0)=[0~-0.25~2.2]^{\rm T},~x(0)=[-0.2~-0.8~0]^{\rm T}$. In addition, the remaining parameter conditions are the same as those in Example 1. The final simulation results are shown in Figure \ref{3}.
	
	From Figure \ref{fig.6}, it can be seen that the autonomous vehicle is able to smoothly track the reference trajectory and subsequently travel stably along the path. From Figure \ref{fig.7} and Figure \ref{fig.8}, it can be determined that when the desired linear velocity is set to $0.2~\rm m/s$, it takes approximately 3 seconds to catch up with the reference trajectory. As can be seen from Figure \ref{fig.9}, the stability curve is very smooth. Furthermore, the curve eventually converges to 0. By conducting simulation experiments on these two cases, it can be observed that the proposed regularized robust predictive controller, which takes into account uncertainties and dynamic quantization, is effective in tracking paths.

\section{Conclusion}

This paper investigates an error control problem in trajectory tracking for autonomous vehicles. To begin with, a linear time-varying error model is established based on the vehicle kinematic model. Considering the bias in practical information measurements, a Bernoulli probability is introduced, and a quantization strategy for dynamic parameter updates is adopted. Subsequently, a quadratic cost function is constructed, and a terminal weight matrix is designed in the penalty function to ensure the Lyapunov stability of the system. In this way, the error control problem is transformed into a predictive control problem with LMI constraints. To address the uncertain parameters in the model, a regularized least squares method is used to solve the optimization problem. After calculating the terminal weight matrix that satisfies the system for the asymptotic stability of the mean square at the corresponding time, the optimal control law is obtained through online iteration to update the system state. Finally, the feasibility of the proposed algorithm is verified through two simulation experiments.

The results indicate that the accuracy of the model is guaranteed with the addition of the dynamic quantizer, but it also increases the computation time and reduces the real-time performance of the controller. If the LMI constraints can be solved offline \cite{liu2020robust}, it may be possible to achieve a better balance between these two factors. Additionally, given the many unpredictable uncertainties and external disturbances that exist during actual driving, this adds more challenges to the theoretical establishment of models. In future research, it may be possible to draw on improved MPC algorithms such as nonlinear MPC \cite{vu2021model} and model-reference adaptive MPC \cite{l2020constrained} to address uncertainties and disturbances in practical problems. However, it should be noted that these methods are not directly applicable to trajectory tracking issues. Careful consideration needs to be given to how these ideas can be borrowed. Thus, this challenge remains obstructive and persistent.

%\begin{ack}
%Place acknowledgments here.
%\end{ack}

\bibliography{ref}             % bib file to produce the bibliography
                                                     % with bibtex (preferred)
                                                   
%\begin{thebibliography}{xx}  % you can also add the bibliography by hand

%\bibitem[Able(1956)]{Abl:56}
%B.C. Able.
%\newblock Nucleic acid content of microscope.
%\newblock \emph{Nature}, 135:\penalty0 7--9, 1956.

%\bibitem[Able et~al.(1954)Able, Tagg, and Rush]{AbTaRu:54}
%B.C. Able, R.A. Tagg, and M.~Rush.
%\newblock Enzyme-catalyzed cellular transanimations.
%\newblock In A.F. Round, editor, \emph{Advances in Enzymology}, volume~2, pages
%  125--247. Academic Press, New York, 3rd edition, 1954.

%\bibitem[Keohane(1958)]{Keo:58}
%R.~Keohane.
%\newblock \emph{Power and Interdependence: World Politics in Transitions}.
%\newblock Little, Brown \& Co., Boston, 1958.

%\bibitem[Powers(1985)]{Pow:85}
%T.~Powers.
%\newblock Is there a way out?
%\newblock \emph{Harpers}, pages 35--47, June 1985.

%\bibitem[Soukhanov(1992)]{Heritage:92}
%A.~H. Soukhanov, editor.
%\newblock \emph{{The American Heritage. Dictionary of the American Language}}.
%\newblock Houghton Mifflin Company, 1992.

%\end{thebibliography}

%\appendix
%\section{A summary of Latin grammar}    % Each appendix must have a short title.
%\section{Some Latin vocabulary}              % Sections and subsections are supported  
                                                                         % in the appendices.
\end{document}