Fully Distributed Optimization based CAV Platooning Control under Linear Vehicle Dynamics

Jinglai Shen, Eswar Kumar H. Kammara, Lili Du Jinglai Shen is with Department of Mathematics and Statistics, University of Maryland Baltimore County, MD 21250, USA. Email: [email protected].Eswar Kumar H. K. is with Department of Mathematics and Statistics, University of Maryland Baltimore County, MD 21250, USA. Email: [email protected].Lili Du is with Department of Civil and Coastal Engineering, University of Florida, Gainesville, FL 32608, USA. Email: [email protected].

(Originally November 23, 2020)

Abstract

This paper develops distributed optimization based, platoon centered CAV car following schemes, motivated by the recent interest in CAV platooning technologies. Various distributed optimization or control schemes have been developed for CAV platooning. However, most existing distributed schemes for platoon centered CAV control require either centralized data processing or centralized computation in at least one step of their schemes, referred to as partially distributed schemes. In this paper, we develop fully distributed optimization based, platoon centered CAV platooning control under the linear vehicle dynamics via the model predictive control approach with a general prediction horizon. These fully distributed schemes do not require centralized data processing or centralized computation through the entire schemes. To develop these schemes, we propose a new formulation of the objective function and a decomposition method that decomposes a densely coupled central objective function into the sum of several locally coupled functions whose coupling satisfies the network topology constraint. We then exploit the formulation of locally coupled optimization and operator splitting methods to develop fully distributed schemes. Control design and stability analysis is carried out to achieve desired traffic transient performance and asymptotic stability. Numerical tests demonstrate the effectiveness of the proposed fully distributed schemes and CAV platooning control.

Keywords: Connected and autonomous vehicle, car following control, distributed algorithm, constrained optimization, control and stability

1 Introduction

The recent advancement of connected and autonomous vehicle (CAV) technologies provides a unique opportunity to mitigate urban traffic congestion through innovative traffic flow control and operations. Supported by advanced sensing, vehicle communication, and portable computing technologies, CAVs can sense, share, and process real-time mobility data and conduct cooperative or coordinated driving. This has led to a surging interest in self-driving technologies. Among a number of emerging self-driving technologies, vehicle platooning technology receives substantial attention. Specifically, the vehicle platooning technology links a group of CAVs through cooperative acceleration or speed control. It allows adjacent group members to travel safely at a higher speed with smaller spacing between them and thus has a great potential to increase lane capacity, improve traffic flow efficiency, and reduce congestion, emission, and fuel consumption [2, 7].

There is extensive literature on CAV platooning control. The widely studied approaches include adaptive cruise control (ACC) [8, 10, 11, 18, 24], cooperative adaptive cruise control (CACC) [14, 15, 17, 22], and platoon centered vehicle platooning control [4, 5, 19, 20]. The first two approaches intend to improve an individual vehicle’s safety, mobility, and string stability rather than systematical performance of the entire platoon, even though enhanced system performance is validated by analysis, simulations, or field experiments. In contrast, the platoon centered approach aims to improve the performance of the entire platoon and seeks a control input that optimizes the platoon’s transient traffic dynamics for a smooth traffic flow while achieving stability and other desired long-time dynamical behaviors.

The platoon centered CAV platooning control often gives rise to sophisticated, large-scale optimal control or optimization problems, and requires extensive computation. In order to successfully implement these control schemes, efficient real-time computation is needed [19]. However, due to high computation load and the absence of roadside computing facilities, centralized computation is either inefficient or infeasible [21]. By leveraging portable computing capability of each vehicle, distributed computing is a favorable option because it has potentials to be more adaptive to different platoon network topologies, be more robust to network malfunctions, and accommodate for communication delays effectively [12, 21]. In spite of these advantages, the development of efficient distributed algorithms to solve platoon centered optimization or optimal control problems in real-time is nontrivial, especially under complicated traffic conditions and constraints.

A number of effective distributed control or optimization schemes have been proposed for CAV platooning [19, 21, 22, 25]. In particular, the recent paper [5] develops model predictive control (MPC) based car-following control schemes for CAV platooning by exploiting transportation, control and optimization methodologies. These control schemes take vehicle constraints, transient dynamics and and asymptotic dynamics of the entire platoon into account, and can be computed in a distributed manner. The paper [4] extends these distributed schemes to a mixed traffic flow including both CAVs and human-driven vehicles. However, to the best of our knowledge, the proposed schemes in [4, 5] as well as many other existing distributed or decentralized schemes [9] either require all vehicles to exchange information with a central component for centralized data processing or perform centralized computation in at least one step of these schemes. We refer to these schemes as partially distributed schemes. In contrast, the distributed schemes developed in this paper do not require centralized data processing or carry out centralized computation through the entire schemes and thus are called fully distributed. Distinct advantages of fully distributed schemes include but are not limited to: (i) no data synchronization is needed such that no central computing equipment is required; and (ii) each vehicle only interacts with its nearby vehicles through a vehicle communication network. Hence, these schemes impose less restriction on vehicle communication networks and can be easily implemented on a wide range of vehicle networks. They are also suitable for a large CAV platoon in remote areas where communication network is unreliable or roadside equipments are scarce. Further, they are more robust to network malfunction or cyber attacks.

In this paper, we develop fully distributed optimization based and platoon centered CAV car following control schemes over a general vehicle communication network. We propose a general $p$ -horizon MPC model subject to the linear vehicle dynamics of the CAVs and various physical or safety constraints. Typically, a fully distributed optimization scheme requires the objective function and constraints of the underlying optimization problem to be decoupled [6]. However, the constrained optimization problem arising from the proposed MPC model does not satisfy this requirement since its objective function is densely coupled and its constraints are locally coupled; see Remark 3.2 for more details. Therefore, this paper develops new techniques to overcome this difficulty.

The main contributions of this paper are summarized as follows:

(1)

We propose a new form of the objective function in the MPC model with new sets of weight matrices. This new formulation facilitates the development of fully distributed schemes and closed loop stability analysis whereas it can achieve desired traffic transient performance of the whole platoon. Based on the new formulation, a decomposition method is developed for the strongly convex quadratic objective function. This method decomposes the central objective function into the sum of locally coupled (strongly) convex quadratic functions, where local coupling satisfies the network topology constraint under a mild assumption on network topology. Along with locally coupled constraints in the MPC model, the underlying optimization model is formulated as a locally coupled convex quadratically constrained quadratic program (QCQP).
(2)

Fully distributed schemes are developed for solving the above-mentioned convex QCQP arising from the MPC model using the techniques of locally coupled optimization and operator splitting methods. Specifically, by introducing copies of local coupling variables of each vehicle, an augmented optimization model is formulated with an additional consensus constraint. A generalized Douglas-Rachford splitting method based distributed scheme is developed, where only local information exchange is needed, leading to a fully distributed scheme. Other operator splitting method based distributed scheme are also discussed.
(3)

The new formulation of the weight matrices and objective function leads to different closed loop dynamics in comparison with that in [5]. Besides, since a general $p$ -horizon MPC is considered, it calls for new stability analysis of the closed loop dynamics. We perform detailed stability analysis and choose suitable weight matrices for desired traffic transient performance for a general horizon length $p$ . In particular, we prove that up to a horizon of $p=3$ , the closed loop dynamic matrix is Schur stable. Extensive numerical tests are carried out to test the proposed distributed schemes under different MPC horizon $p$ ’s and to evaluate the closed loop stability and performance.

The rest of the paper is organized as follows. Section 2 introduces the linear vehicle dynamics, state and control constraints, and vehicle communication networks. The model predictive control model with a general prediction horizon $p$ is proposed and formulated as a constrained optimization problem in Section 3; fundamental properties of this optimization problem are established. Section 4 develops fully distributed schemes by exploiting a decomposition method for the central quadratic objective function, locally coupled optimization, and operator splitting methods. Control design and stability analysis for the closed loop dynamics is presented in Section 5 with numerical results given in Section 6. Finally, conclusions are given in Section 7.

2 Vehicle Dynamics, Constraints, and Communication Networks

We consider a platoon consisting of vehicles on a roadway, where the (uncontrolled) leading vehicle is labeled by the index 0 and its $n$ following CAVs are labeled by the indices $i=1,\ldots,n$ , respectively. Let $x_{i},v_{i}$ denote the longitudinal position and speed of the $i$ th vehicle, respectively. Let $\tau>0$ be the sampling time, and each time interval is given by $[k\tau,(k+1)\tau)$ for $k\in\mathbb{Z}_{+}:=\{0,1,2,\ldots\}$ . We consider the following kinematic model for linear vehicle dynamics that is widely adopted in system-level studies with the acceleration $u_{i}(k)$ as the control input:

\displaystyle x_{i}(k+1)\,=\,x_{i}(k)+\tau v_{i}(k)+\frac{\tau^{2}}{2}u_{i}(k),\qquad v_{i}(k+1)\,=\,v_{i}(k)+\tau u_{i}(k).

(1)

State and control constraints. Each vehicle in a platoon is subject to several important state and control constraints. For each $i=1,\ldots,n$ ,

(i)

Control constraints: $a_{\min}\leq u_{i}\leq a_{\max}$ , where $a_{\min}<0$ and $a_{\max}>0$ are pre-specified acceleration and deceleration bounds for a vehicle.
(ii)

Speed constraints: $v_{\min}\leq v_{i}\leq v_{\max}$ , where $0\leq v_{\min}<v_{\max}$ are pre-specified bounds on longitudinal speed for a vehicle;
(iii)

Safety distance constraints: these constraints guarantee sufficient spacing between neighboring vehicles to avoid collision even if the leading vehicle comes to a sudden stop. This gives rise to the safety distance constraint of the following form:

$x_{i-1}-x_{i}\,\geq\,L+r\cdot v_{i}-\frac{(v_{i}-v_{\min})^{2}}{2a_{\min}},$ (2)

where $L>0$ is a constant depending on vehicle length, and $r$ is the reaction time.

Note that constraints (i) and (ii) are decoupled across vehicles, whereas the safety distance constraint (iii) is state-control coupled since such a constraint involves control inputs of two vehicles. This yields challenges to distribution computation. Further, we consider the identical bounds of the acceleration and deceleration and ignore their variations due to road surface condition changes in this paper, although the proposed approach can handle a general case with different acceleration or deceleration bounds.

Communication network topology. We consider a general communication network whose topology is modeled by a graph $\mathcal{G}(\mathcal{V},\mathcal{E})$ , where $\mathcal{V}=\{1,2,\ldots,n\}$ is the set of nodes where the $i$ th node corresponds to the $i$ th CAV, and $\mathcal{E}$ is the set of edges connecting two nodes in $\mathcal{V}$ . Let $\mathcal{N}_{i}$ denote the set of neighbors of node $i$ , i.e., $\mathcal{N}_{i}=\{j\,|\,(i,j)\in\mathcal{E}\}$ . The following assumption on the communication network topology is made throughout the paper:

$\bf A.1$

The graph $\mathcal{G}(\mathcal{V},\mathcal{E})$ is undirected and connected. Further, two neighboring vehicles form a bidirectional edge of the graph, i.e., $(1,2),(2,3),\ldots,(n-1,n)\in\mathcal{E}$ .

Since the graph is undirected, for any $i,j\in\mathcal{V}$ with $i\neq j$ , $(i,j)\in\mathcal{E}$ means that there exists an edge between node $i$ and node $j$ . In other words, vehicle $i$ can receive information from vehicle $j$ and send information to vehicle $j$ , and so does vehicle $j$ . We also assume that the first vehicle can detect and receive $x_{0}$ , $v_{0}$ and $u_{0}$ from the leading vehicle. The setting given by $\bf A.1$ includes many widely used communication networks of CAV platoons, e.g., immediate-preceding, multiple-preceding, and preceding-and-following networks [23].

3 Model Predictive Control for CAV Platooning Control

We exploit the model predictive control (MPC) approach for car following control of a platoon of CAVs. Let $\Delta$ be the desired (constant) spacing between two adjacent vehicles, and $(x_{0},v_{0},u_{0})$ be the position, speed, and control input of the leading vehicle, respectively. Define the vectors: (i) $z(k):=\big{(}x_{0}-x_{1}-\Delta,\ldots,x_{n-1}-x_{n}-\Delta\big{)}(k)\in\mathbb{R}^{n}$ , representing the relative spacing error; (ii) $z^{\prime}(k):=\big{(}v_{0}-v_{1},\ldots,v_{n-1}-v_{n}\big{)}(k)\in\mathbb{R}^{n}$ , representing the relative speed between adjacent vehicles; and (iii) $u(k):=\big{(}u_{1},\ldots,u_{n}\big{)}(k)\in\mathbb{R}^{n}$ , representing the control input. Further, let $w_{i}(k):=u_{i-1}(k)-u_{i}(k)$ for each $i=1,\ldots,n$ , and $w(k):=\big{(}w_{1},\ldots,w_{n}\big{)}(k)\in\mathbb{R}^{n}$ , which stands for the difference of control input between adjacent vehicles. Hence, for any $k\in\mathbb{Z}_{+}$ , $u(k)=-S_{n}w(k)+u_{0}(k)\cdot\mathbf{1}$ , where $\mathbf{1}:=(1,\ldots,1)^{T}$ is the vector of ones, and

S_{n}\,:=\,\begin{bmatrix}1&0&0&\ldots&0\\ 1&1&0&\ldots&0\\ \vdots&\vdots&\ddots&\ddots&\vdots\\ 1&1&\ldots&1&0\\ 1&1&\ldots&1&1\\ \end{bmatrix}\in\mathbb{R}^{n\times n},\qquad S^{-1}_{n}\,=\,\begin{bmatrix}1&&&&\\ -1&1&&&\\ &\ddots&\ddots&&\\ &&-1&1&\\ &&&-1&1\\ \end{bmatrix}\in\mathbb{R}^{n\times n}.

(3)

Given a prediction horizon $p\in\mathbb{N}$ , the $p$ -horizon MPC control is determined by solving the following constrained optimization problem at each $k\in\mathbb{Z}_{+}$ , involving all vehicles’ control inputs for given feasible state $(x_{i}(k),v_{i}(k))^{n}_{i=1}$ and $(x_{0}(k),v_{0}(k),u_{0}(k))$ at time $k$ subject to the vehicle dynamics model (1):

	$\displaystyle\ \mbox{minimize}\ \ J(u(k),\ldots,u(k+p-1)):=$		(4)
	$\displaystyle\ \frac{1}{2}\sum^{p}_{s=1}\Big{(}\underbrace{\tau^{2}u^{T}(k+s-1)S^{-T}_{n}Q_{w,s}S^{-1}_{n}u(k+s-1)}_{\mbox{ride comfort}}+\underbrace{z^{T}(k+s)Q_{z,s}z(k+s)+(z^{\prime}(k+s))^{T}Q_{z^{\prime},s}z^{\prime}(k+s)}_{\mbox{traffic stability and smoothness}}\Big{)}$

subject to: for each $i=1,\ldots,n$ and each $s=1,\ldots,p$ ,

	$\displaystyle a_{\min}$	$\displaystyle\leq u_{i}(k+s-1)\ \leq a_{\max},$	$\displaystyle\qquad v_{\min}\,\leq\,v_{i}(k+s)\ \leq\ v_{\max},$		(5)
		$\displaystyle x_{i-1}(k+s)-x_{i}(k+s)$	$\displaystyle\geq L+r\cdot v_{i}(k+s)-\frac{(v_{i}(k+s)-v_{\min})^{2}}{2a_{\min}},$		(6)

where $Q_{z,s}$ , $Q_{z^{\prime},s}$ and $Q_{w,s}$ are $n\times n$ symmetric positive semidefinite weight matrices to be discussed soon. Note that when $p>1$ , $(x_{0}(k+s+1),v_{0}(k+s+1),u_{0}(k+s))$ are unknown at time $k$ for $s=1,\ldots,p-1$ . In this case, we assume that $u_{0}(k+s)=u_{0}(k)$ for all $s=1,\ldots,p-1$ and use these $u_{0}(k+s)$ ’s and the vehicle dynamics model (1) to predict $(x_{0}(k+s+1),v_{0}(k+s+1))$ for $s=1,\ldots,p-1$ . Here we assume that $u_{0}(k)$ represents the actual acceleration of the leading vehicle at time $k$ .

Remark 3.1.

The three terms in the objective function $J$ intend to minimize traffic flow oscillations via mild control: the first term penalizes the magnitude of control, whereas the second and last terms penalize the variations of the relative spacing and relative speed, respectively. The presence of the matrix $S_{n}$ in the first term is due to the coupled vehicle dynamics through the CAV platoon. To illustrate this, let $\widetilde{w}(k+s-1):=w(k+s-1)-u_{0}(k)\cdot\mathbf{e}_{1}$ for $s=1,\ldots,p$ . Thus $\widetilde{w}(k+s-1)=-S^{-1}_{n}u(k+s-1)$ for each $s=1,\ldots,p$ . Therefore, the first term in $J$ satisfies $\tau^{2}u^{T}(k+s-1)S^{-T}_{n}Q_{w,s}S^{-1}_{n}u(k+s-1)=\tau^{2}\widetilde{w}^{T}(k+s-1)\,Q_{w,s}\,\widetilde{w}(k+s-1)$ for each $s$ .

The weight matrices $Q_{z,s}$ , $Q_{z^{\prime},s}$ , and $Q_{w,s}$ , $s=1,\ldots,p$ determine transient and asymptotic dynamics, and they depend on vehicle network topologies and can be chosen by stability analysis and transient dynamics criteria of the closed loop system. To develop fully distributed schemes for a broad class of vehicle network topologies and to facilitate control design and analysis, we make the following blanket assumption on $Q_{z,s}$ , $Q_{z^{\prime},s}$ , and $Q_{w,s}$ throughout the rest of the paper:

$\bf A.2$

For each $s=1,\ldots,p$ , $Q_{z,s}$ and $Q_{z^{\prime},s}$ are diagonal and positive semidefinite (PSD), and $Q_{w,s}$ is diagonal and positive definite (PD).

The reasons for considering this class of diagonal positive semidefinite or positive definite weight matrices are three folds: (i) Diagonal matrices have a simpler interpretation in transportation engineering so that the selection of such matrices is easier to practitioners. For instance, the diagonal $Q_{z,s}$ and $Q_{z^{\prime},s}$ mean that one imposes penalties on each element of $z(k+s)$ and $z^{\prime}(k+s)$ without considering their coupling. Further, by suitably choosing the weight matrices $Q_{w,s}$ , it can be shown that the ride comfort term in equation (4), which corresponds to acceleration of CAVs, is similar to imposing direct penalties on $u_{i}$ ’s, which simplifies control design. (ii) This class of weight matrices facilitates the development of fully distributed schemes for general vehicle network topologies. (iii) Closed-loop stability and performance analysis is relatively simpler (although still nontrivial) when using this class of weight matrices. The detailed discussions of choosing diagonal, positive semidefinite or positive definite weight matrices for satisfactory closed loop dynamics will be given in Section 5.

The sequential feasibility has been established in [4, 5] for the MPC model (4) when $r\geq\tau$ . Define $\mathcal{P}((x_{i},v_{i})^{n}_{i=0},u_{0})\,:=\,\{u\in\mathbb{R}^{n}\,|\,a_{\min}\leq u_{i}\leq a_{\max},\,v_{\min}\leq v_{i}+\tau u_{i}\leq v_{\max},\,h_{i}(u)\leq 0,\ \forall\,i=1,\ldots,n\},$ where $h_{i}(u):=L+r(v_{i}+\tau u_{i})-\frac{(v_{i}+\tau u_{i}-v_{\min})^{2}}{2a_{\min}}+(x_{i}-x_{i-1})+\tau(v_{i}-v_{i-1})+\frac{\tau^{2}}{2}[u_{i}-u_{i-1}]$ for each $i=1,\ldots,n$ . Specifically, the sequential feasibility implies that for any feasible $x_{i}(k),v_{i}(k),u_{0}(k)$ at time $k$ , the constraint set $\mathcal{P}((x_{i}(k),v_{i}(k))^{n}_{i=0},u_{0}(k))$ is non-empty such that the MPC model (4) has a solution $u_{*}(k)$ such that the constraint set $\mathcal{P}((x_{i}(k+1),v_{i}(k+1))^{n}_{i=0},u_{0}(k+1))$ is non-empty. Using this result, we show below that under a mild assumption, the constraint sets of the MPC model have nonempty interior for any MPC horizon $p\in\mathbb{N}$ . This result is important to the development of distributed algorithms.

Corollary 3.1.

Consider the linear vehicle dynamics (1) and assume $r\geq\tau$ . Suppose the leading vehicle is such that $(v_{0}(k),u_{0}(k))$ is feasible and $v_{0}(k)>v_{\min}$ for all $k\in\mathbb{Z}_{+}$ . Then the constraint set of the $p$ -horizon MPC model (4) has nonempty interior at each $k$ .

Proof.

Fix an arbitrary $k\in\mathbb{Z}_{+}$ . Since $v_{0}(k)>v_{\min}$ , it follows from [5, Proposition 3.1] that there exists a vector denoted by $\widehat{u}(k)$ in the interior of the set $\mathcal{P}((x_{i}(k),v_{i}(k))^{n}_{i=0},u_{0}(k))$ . Let $x_{i}(k+1)$ and $v_{i}(k+1)$ be generated by $\widehat{u}(k)$ (and $(x_{i}(k),v_{i}(k))^{n}_{i=0},u_{0}(k)$ ). Since $v_{0}(k+1)>v_{\min}$ , we deduce via [5, Proposition 3.1] again that there exists a vector denoted by $\widehat{u}(k+1)$ in the interior of the constraint set $\mathcal{P}((x_{i}(k+1),v_{i}(k+1))^{n}_{i=0},u_{0}(k+1))$ . Continuing this process in $p$ -steps, we derive the existence of an interior point in the constraint set of the $p$ -horizon MPC model (4). ∎

3.1 Constrained MPC Optimization Model

Consider the constrained MPC optimization model (4) at an arbitrary but fixed time $k\in\mathbb{Z}_{+}$ subject to the linear vehicle dynamics (1). In view of the following results: for each $s=1,\ldots,p$ ,

	$\displaystyle v_{i}(k+s)$	$\displaystyle=v_{i}(k)+\tau\sum^{s-1}_{j=0}u_{i}(k+j),\quad z^{\prime}(k+s)=z^{\prime}(k)+\tau\sum^{s-1}_{j=0}w(k+j),$
	$\displaystyle z(k+s)$	$\displaystyle=z(k)+s\tau z^{\prime}(k)+\tau^{2}\sum^{s-1}_{j=0}\frac{2(s-j)-1}{2}w(k+j),\quad w(k+s)=S^{-1}_{n}\big{[}-u(k+s)+u_{0}(k)\cdot\mathbf{1}\big{]},$

we formulate (4) as the constrained convex minimization problem (where we omit $k$ since it is fixed):

\begin{array}[]{ll}\mbox{minimize}&J({\mbox{\bf u}}):=\frac{1}{2}{\mbox{\bf u}}^{T}W{\mbox{\bf u}}+c^{T}{\mbox{\bf u}}+\gamma,\\ \mbox{subject to}&{\mbox{\bf u}}_{i}\in\mathcal{X}_{i},\quad(H_{i}({\mbox{\bf u}}))_{s}\leq 0,\quad\forall\,i=1,\ldots,n,\quad\forall\,s=1,\ldots,p,\end{array}

(7)

where ${\mbox{\bf u}}:=({\mbox{\bf u}}_{1},\ldots,{\mbox{\bf u}}_{n})\in\mathbb{R}^{np}$ with ${\mbox{\bf u}}_{i}:=(u_{i}(k),\ldots,u_{i}(k+p-1))\in\mathbb{R}^{p}$ , $W$ is a PD matrix to be shown in Lemma 3.1 below, $c\in\mathbb{R}^{np}$ , $\gamma\in\mathbb{R}$ , each $\mathcal{X}_{i}:=\{z\in\mathbb{R}^{p}\,|\,a_{\min}\cdot\mathbf{1}\leq z\leq a_{\max}\cdot\mathbf{1},\ (v_{\min}-v_{i}(k))\cdot\mathbf{1}\leq\tau S_{p}z\leq(v_{\max}-v_{i}(k))\cdot\mathbf{1}\}$ is a polyhedral set, and each $(H_{i}(\cdot))_{s}$ is a convex quadratic function characterizing the safety distance given by (12). Here $S_{p}$ is the $p\times p$ matrix of the form given by (3). Further, $\mathbf{u}_{0}:=u_{0}(k)\cdot\mathbf{1}_{p}\in\mathbb{R}^{p}$ for the given $u_{0}(k)$ . An important property of the matrix $W$ in (7) is given below.

Lemma 3.1.

Suppose that $Q_{z,s}$ and $Q_{z^{\prime},s}$ are PSD and $Q_{w,s}$ are PD for all $s=1,\ldots,p$ (but not necessarily diagonal). Then the matrix $W$ in (7) is PD.

Proof.

Let u be an arbitrary nonzero vector in $\mathbb{R}^{np}$ . Since $J(\cdot)$ is quadratic, we have $\frac{1}{2}{\mbox{\bf u}}^{T}W{\mbox{\bf u}}=\lim_{\lambda\rightarrow\infty}\frac{J(\lambda{\mbox{\bf u}})}{\lambda^{2}}$ . In view of the equivalent formulation of $J(\cdot)$ given by (4), we deduce that for any $\lambda>0$ , $J(\lambda{\mbox{\bf u}})=J(\lambda u(k),\ldots,\lambda u(k+p-1))\geq\frac{\lambda^{2}}{2}\sum^{p}_{s=1}\tau^{2}u^{T}(k+s-1)S^{-T}_{n}Q_{w,s}S^{-1}_{n}u(k+s-1)>0$ , where the first inequality follows from the fact that $Q_{z,s}$ and $Q_{z^{\prime},s}$ are PSD, and the second inequality holds because $Q_{w,s}$ , and thus $S^{-T}_{n}Q_{w,s}S^{-1}_{n}$ , are PD. Therefore, $\frac{J(\lambda{\mbox{\bf u}})}{\lambda^{2}}\geq\frac{1}{2}\sum^{p}_{s=1}\tau^{2}u^{T}(k+s-1)S^{-1}_{n}Q_{w,s}S^{-1}_{n}u(k+s-1)>0$ , leading to $\frac{1}{2}{\mbox{\bf u}}^{T}W{\mbox{\bf u}}\geq\frac{1}{2}\sum^{p}_{s=1}\tau^{2}u^{T}(k+s-1)S^{-1}_{n}Q_{w,s}S^{-1}_{n}u(k+s-1)>0$ . Hence, $W$ is PD. ∎

To establish the closed form expressions of the matrix $W$ and the vector $c$ in (7), we define the following matrices for any $i,j\in\{1,\ldots,p\}$ :

V_{i,j}\,:=\,S^{-T}_{n}\left[\sum^{p}_{s=\max(i,j)}\Big{(}\frac{\tau^{4}}{4}[2(s-i)+1]\cdot[2(s-j)+1]Q_{z,s}+\tau^{2}Q_{z^{\prime},s}\Big{)}\right]S^{-1}_{n}\in\mathbb{R}^{n\times n}.

Clearly, $V_{i,j}=V_{j,i}$ for any $i,j$ . Moreover, let $\widetilde{Q}_{w,s}:=S^{-T}_{n}Q_{w,s}S^{-1}_{n}$ for $s=1,\ldots,p$ . Hence, the symmetric matrix $W$ is given by $W=E^{T}VE$ , where

V\,=\,\begin{bmatrix}V_{1,1}+\tau^{2}\widetilde{Q}_{w,1}&V_{1,2}&V_{1,3}&\cdots&\cdots&V_{1,p}\\ V_{2,1}&V_{2,2}+\tau^{2}\widetilde{Q}_{w,2}&V_{2,3}&\cdots&\cdots&V_{2,p}\\ \cdots&&\cdots&&\cdots&\\ \cdots&&\cdots&&\cdots&\\ V_{p,1}&V_{p,2}&V_{p,3}&\cdots&\cdots&V_{p,p}+\tau^{2}\widetilde{Q}_{w,p}\end{bmatrix}\in\mathbb{R}^{np\times np},

(8)

and $E\in\mathbb{R}^{np\times np}$ is the permutation matrix satisfying

\begin{bmatrix}u(k)\\ u(k+1)\\ \vdots\\ u(k+p-1)\end{bmatrix}\,=\,E\begin{bmatrix}\mathbf{u}_{1}\\ \mathbf{u}_{2}\\ \vdots\\ \mathbf{u}_{n}\end{bmatrix}.

Specifically, the $(i,j)$ -entry of the matrix $E$ is given by

E_{i,j}=\left\{\begin{array}[]{ll}1&\mbox{if}\ \ i=n\cdot k+s,\ j=p\cdot(s-1)+k+1,\quad\mbox{for}\ \ k=0,\ldots,p-1,\ s=1,\ldots,n;\\ 0,&\mbox{otherwise}.\end{array}\right.

(9)

In particular, when $p=1$ , $E=I_{n}$ .

For a fixed $k\in\mathbb{Z}_{+}$ , we also define for each $s=1,\ldots,p$ ,

\displaystyle d_{s}(k)\,:=\,z(k)+s\tau z^{\prime}(k)+\tau^{2}\sum^{s-1}_{j=0}\frac{2(s-j)-1}{2}S^{-1}_{n}\cdot\mathbf{1}\cdot u_{0}(k),\quad f_{s}(k)

\displaystyle\,:=\,z^{\prime}(k)+\tau\sum^{s-1}_{j=0}S^{-1}_{n}\cdot\mathbf{1}\cdot u_{0}(k).

In light of $S^{-1}_{n}$ given by (3), we have $S^{-1}_{n}\cdot\mathbf{1}=\mathbf{e}_{1}$ . Therefore, we obtain

d_{s}(k)=z(k)+s\tau z^{\prime}(k)+\frac{\tau^{2}}{2}s^{2}\mathbf{e}_{1}u_{0}(k),\qquad f_{s}(k)=z^{\prime}(k)+\tau s\mathbf{e}_{1}u_{0}(k).

(10)

In view of

z(k+s)=d_{s}(k)-\tau^{2}\sum^{s-1}_{j=0}\frac{2(s-j)-1}{2}S^{-1}_{n}u(k+j),\qquad z^{\prime}(k+s)=f_{s}(k)-\tau\sum^{s-1}_{j=0}S^{-1}_{n}u(k+j),

the linear terms in the objective function $J$ are given by

-\sum^{p}_{i=1}\left(\,\sum^{p}_{s=i}\Big{[}\frac{\tau^{2}}{2}[2(s-i)+1]d^{T}_{s}(k)Q_{z,s}+\tau f^{T}_{s}(k)Q_{z^{\prime},s}\Big{]}\,\right)S^{-1}_{n}\cdot u(k+i-1).

(11)

Using the permutation matrix $E$ given in (9), we can write the above linear terms as $c^{T}\mathbf{u}\,=\,\sum^{n}_{i=1}c^{T}_{\mathcal{I}_{i}}\mathbf{u}_{i},$ where $c_{\mathcal{I}_{i}}$ is the subvector of $c$ corresponding to $\mathbf{u}_{i}$ . Since $Q_{z,s}$ and $Q_{z^{\prime},s}$ are diagonal, it is easy to obtain the following lemma via $d_{s}(k),f_{s}(k)$ in (10) and the structure of $S^{-1}_{n}$ given by (3):

Lemma 3.2.

Consider the vector $c=(c_{\mathcal{I}_{1}},\ldots,c_{\mathcal{I}_{n}})$ given above. Then for each $i=1,\ldots,n$ , the subvector $c_{\mathcal{I}_{i}}$ depends only on $z_{i}(k),z^{\prime}_{i}(k),z_{i+1}(k),z^{\prime}_{i+1}(k)$ ’s for $i=1,\ldots,n-1$ , and $c_{\mathcal{I}_{n}}$ depends only on $z_{n}(k),z^{\prime}_{n}(k)$ . Further, only $c_{\mathcal{I}_{1}}$ depends on $u_{0}(k)$ .

The above lemma shows that each $c_{\mathcal{I}_{i}}$ only depends on the information of the adjacent vehicles of vehicle $i$ , and thus can be easily established from any vehicle network. This property is important for developing fully distributed schemes to be shown in Section 4.3.

To find the vector form of the safety constraint, we note that for $s=1,\ldots,p$ ,

x_{i}(k+s)\,=\,x_{i}(k)+s\tau v_{i}(k)+\tau^{2}\sum^{s-1}_{j=0}\frac{2(s-j)-1}{2}u_{i}(k+j),\quad v_{i}(k+s)=v_{i}(k)+\tau\sum^{s-1}_{j=0}u_{i}(k+j).

The safety distance constraint for the $i$ -th vehicle at time $k$ is given by: for $s=1,\ldots,p$ ,

$\displaystyle 0$	$\displaystyle\geq-\Big{[}\,x_{i-1}(k)+s\tau v_{i-1}(k)-(x_{i}(k)+s\tau v_{i}(k))\,\Big{]}-\tau^{2}\sum^{s-1}_{j=0}\frac{2(s-j)-1}{2}[u_{i-1}(k+j)-u_{i}(k+j)]$
	$\displaystyle+\,L+rv_{i}(k)+r\tau\sum^{s-1}_{j=0}u_{i}(k+j)$	(12)
	$\displaystyle-\frac{1}{2a_{\min}}\Big{[}\tau^{2}\big{(}\sum^{s-1}_{j=0}u_{i}(k+j)\big{)}^{2}+2\tau(v_{i}(k)-v_{\min})\sum^{s-1}_{j=0}u_{i}(k+j)+(v_{i}(k)-v_{\min})^{2}\Big{]}:=(H_{i}(\mathbf{u}_{i-1},\mathbf{u}_{i}))_{s},$

where $(H_{i}(\cdot,\cdot))_{s}$ is a convex quadratic function for each $s=1,\ldots,p$ . Hence, the set $\mathcal{Z}_{i}:=\{\mathbf{u}\in\mathbb{R}^{np}\,|\,(H_{i}(\mathbf{u}_{i-1},\mathbf{u}_{i}))_{s}\leq 0,\,\forall\,i=1,\ldots,p\}$ is closed and convex. The problem (7) becomes $\min_{\mathbf{u}}J(\mathbf{u})$ subject to $\mathbf{u}_{i}\in\mathcal{X}_{i}$ and $\mathbf{u}\in\mathcal{Z}_{i}$ for all $i=1,\ldots,n$ , which is a convex quadratically constrained quadratic program (QCQP) and can be solved via a second-order cone program or a semi-definite program in the centralized manner.

Remark 3.2.

The above results show that each $\mathcal{X}_{i}$ ’s are decoupled from the other vehicles, whereas the constraint function $H_{i}$ for vehicle $i$ is locally coupled with its neighboring vehicles. Specifically, $H_{i}$ depends not only on $\mathbf{u}_{i}$ but also on $\mathbf{u}_{i-1}$ of vehicle $i-1$ , which can exchange information with vehicle $i$ . We will explore this local coupling property to develop fully distributed schemes for solving (7).

4 Fully Distributed Algorithms for Constrained Optimization in MPC

We develop fully distributed algorithms for solving the underlying optimization problem given by (7) at each time $k$ using the techniques of locally coupled convex optimization and operator splitting methods.

4.1 Brief Overview of Locally Coupled Convex Optimization

One of major techniques for developing fully distributed schemes for the underlying optimization problem given by (7) is to formulate it as a locally coupled convex optimization problem [6]. To be self-contained, we briefly describe its formulation.

Consider a multi-agent network of $n$ agents whose communication is characterized by a connected and undirect graph $\mathcal{G}(\mathcal{V},\mathcal{E})$ , where $\mathcal{V}=\{1,\ldots,n\}$ is the set of agents, and $\mathcal{E}$ denotes the set of edges. For $i\in\mathcal{V}$ , let $\mathcal{N}_{i}$ be the set of neighbors of agent $i$ , i.e., $\mathcal{N}_{i}=\{j\,|\,(i,j)\in\mathcal{E}\}$ . Let $\{\mathcal{I}_{1},\ldots,\mathcal{I}_{n}\}$ be a disjoint union of the index set $\{1,\ldots,N\}$ . Hence, for any $x\in\mathbb{R}^{N}$ , $(x_{\mathcal{I}_{i}})^{n}_{i=1}$ forms a partition of $x$ . We call $x_{\mathcal{I}_{i}}$ a local variable of each agent $i$ . For each $i$ , define $\widehat{x}_{i}:=\big{(}x_{\mathcal{I}_{i}},(x_{\mathcal{I}_{j}})_{j\in\mathcal{N}_{i}}\big{)}\in\mathbb{R}^{n_{i}}$ . Hence, for each $i$ , $\widehat{x}_{i}$ contains the local variable $x_{\mathcal{I}_{i}}$ and the variables from its neighboring agents (or locally coupled variables). Consider the convex optimization problem

(P):\qquad\min_{x\in\mathbb{R}^{N}}J(x),\qquad\quad\mbox{ where }\quad J(x)\,:=\,\sum^{n}_{i=1}J_{i}(\widehat{x}_{i}),

where $J_{i}:\mathbb{R}^{n_{i}}\rightarrow\mathbb{R}\cup\{+\infty\}$ is an extended-valued, proper, and lower semicontinuous convex function for each $i$ . Clearly, each $J_{i}$ is locally coupled such that the problem $(P)$ bears the name of “locally coupled convex optimization”. Although the problem $(P)$ is seemingly unconstrained, it does include constrained convex optimization since $J_{i}$ may contain the indicator function of a closed convex set. To impose the locally coupled convex constraint explicitly, the problem $(P)$ can be equivalently written as:

(P^{\prime}):\qquad\min_{x\in\mathbb{R}^{N}}\sum^{n}_{i=1}\widehat{J}_{i}(\widehat{x}_{i}),\quad\mbox{subject to }\quad\widehat{x}_{i}\in\mathcal{C}_{i},\quad\forall\ i=1,\ldots,n,

(13)

where for each $i$ , $\widehat{J}_{i}:\mathbb{R}^{n_{i}}\rightarrow\mathbb{R}$ is a real-valued convex function, and $\mathcal{C}_{i}\subseteq\mathbb{R}^{n_{i}}$ is a closed convex set.

By introducing copies of the locally coupled variables for each agent and imposing certain consensus constraints on these copies, the paper [6] formulates the problem $(P^{\prime})$ (or equivalently $(P)$ ) as a separable consensus convex optimization problems. Under suitable assumptions, Douglas-Rachford and other operator splitting based distributed schemes are developed; details can be found in [6].

4.2 Decomposition of a Strongly Convex Quadratic Objective Function

The framework of locally coupled optimization problems requires that both an objective function and constraints are expressed in a locally coupled manner. Especially, the central objective function in (13) is expected to be written as the sum of a few locally coupled functions preserving certain desired properties, e.g., the (strong) convexity if the central objective function is so, where local coupling satisfies network topology constraints. While the constraints of the problem (7) have been shown to be locally coupled (cf. Remark 3.2), the central strongly convex quadratic objective function, particularly its quadratic term $\frac{1}{2}\mathbf{u}^{T}W\mathbf{u}$ , is highly coupled and thus need to be decomposed into the sum of locally coupled (strongly) convex quadratic functions, where the local coupling should satisfy the network topology constraint. In this subsection, we address this decomposition problem under a mild assumption on network topology. Specifically, our results yield a decomposition into convex and strongly convex functions.

We start from a slightly general setting. Let $\boldsymbol{\lambda}:=(\lambda_{1},\ldots,\lambda_{n})\in\mathbb{R}^{n}$ and $\Lambda={\rm diag}(\boldsymbol{\lambda})=\mbox{diag}(\lambda_{1},\ldots,\lambda_{n})$ be a diagonal matrix, i.e., $\boldsymbol{\lambda}$ is the vector representation of the diagonal entries of $\Lambda$ . Therefore, the following matrix is tridiagonal:

S^{-T}_{n}\Lambda S^{-1}_{n}=\begin{bmatrix}\lambda_{1}+\lambda_{2}&-\lambda_{2}&&&\\ -\lambda_{2}&\lambda_{2}+\lambda_{3}&-\lambda_{3}&&\\ &\ddots&\ddots&\ddots&\\ &&-\lambda_{n-1}&\lambda_{n-1}+\lambda_{n}&-\lambda_{n}\\ &&&-\lambda_{n}&\lambda_{n}\end{bmatrix}.

(14)

Consider a general $p\in\mathbb{N}$ . Let $\Theta$ be a symmetric block diagonal matrix given by

\Theta=\begin{bmatrix}\Theta_{1,1}&\Theta_{1,2}&\cdots&\cdots&\Theta_{1,p}\\ \Theta_{2,1}&\Theta_{2,2}&\cdots&\cdots&\Theta_{2,p}\\ \cdots&&\cdots&&\cdots&\\ \cdots&&\cdots&&\cdots&\\ \Theta_{p,1}&\Theta_{p,2}&\cdots&\cdots&\Theta_{p,p}\end{bmatrix}\in\mathbb{R}^{np\times np},

where $\Theta_{i,j}={\rm diag}(\boldsymbol{\theta}_{i,j})\in\mathbb{R}^{n\times n}$ is diagonal for some $\boldsymbol{\theta}_{i,j}\in\mathbb{R}^{n}$ , and $\boldsymbol{\theta}_{i,j}=\boldsymbol{\theta}_{j,i}$ for all $i,j=1,\ldots,p$ . Let $(\boldsymbol{\theta}_{i,j})_{k}$ denote the $k$ th entry of the vector $\boldsymbol{\theta}_{i,j}$ . For each $i=1,\ldots,n$ , define the matrix

U_{i}\,:=\,\begin{bmatrix}(\boldsymbol{\theta}_{1,1})_{i}&(\boldsymbol{\theta}_{1,2})_{i}&\cdots&(\boldsymbol{\theta}_{1,p})_{i}\\ (\boldsymbol{\theta}_{2,1})_{i}&(\boldsymbol{\theta}_{2,2})_{i}&\cdots&(\boldsymbol{\theta}_{2,p})_{i}\\ \cdots&\cdots&\cdots&\\ (\boldsymbol{\theta}_{p,1})_{i}&(\boldsymbol{\theta}_{p,2})_{i}&\cdots&(\boldsymbol{\theta}_{p,p})_{i}\end{bmatrix}\in\mathbb{R}^{p\times p}.

(15)

It can be shown that $\Theta=E^{T}{\rm diag}(U_{1},\ldots,U_{n})E$ , where $E$ is the permutation matrix given by (9). Hence, $\Theta$ is PD (resp. PSD) if and only if each $U_{i}$ is PD (resp. PSD).

Let

\displaystyle V

\displaystyle=

\displaystyle\underbrace{\begin{bmatrix}S^{-T}_{n}&&&&\\ &S^{-T}_{n}&&&\\ &&\ddots&&\\ &&&S^{-T}_{n}\end{bmatrix}}_{:=\mathbf{S}^{-T}}\Theta\underbrace{\begin{bmatrix}S^{-1}_{n}&&&&\\ &S^{-1}_{n}&&&\\ &&\ddots&&\\ &&&S^{-1}_{n}\end{bmatrix}}_{:=\mathbf{S}^{-1}}=\begin{bmatrix}V_{1,1}&V_{1,2}&\cdots&V_{1,p}\\ V_{2,1}&V_{2,2}&\cdots&V_{2,p}\\ \cdots&&\cdots&\cdots\\ V_{p,1}&V_{p,2}&\cdots&V_{p,p}\\ \end{bmatrix},

where $V_{i,j}:=S^{-T}_{n}\Theta_{i,j}S^{-1}_{n}$ is symmetric and tridiagonal. Letting $E$ be the permutation matrix given by (9), a straightforward computation shows that $E^{T}VE$ is a symmetric block tridiagonal matrix given by

W=E^{T}VE=\begin{bmatrix}W_{1,1}&W_{1,2}&&&\\ W_{2,1}&W_{2,2}&W_{2,3}&&\\ &\ddots&\ddots&\ddots&\\ &&W_{n-1,n-2}&W_{n-1,n-1}&W_{n-1,n}\\ &&&W_{n,n-1}&W_{n,n}\end{bmatrix}\in\mathbb{R}^{np\times np},

where each $W_{i,j}\in\mathbb{R}^{p\times p}$ is symmetric and $W_{i,j}=W_{j,i}$ . Furthermore, for each $i=1,\ldots,n$ and $j\in\{i,i+1\}$ ,

W_{i,j}=\begin{bmatrix}(V_{1,1})_{i,j}&(V_{1,2})_{i,j}&\cdots&(V_{1,p})_{i,j}\\ (V_{2,1})_{i,j}&(V_{2,2})_{i,j}&\cdots&(V_{2,p})_{i,j}\\ \cdots&\cdots&\cdots&\\ (V_{p,1})_{i,j}&(V_{p,2})_{i,j}&\cdots&(V_{p,p})_{i,j}\end{bmatrix}\in\mathbb{R}^{p\times p},

where $(V_{r,s})_{i,j}$ denotes the $(i,j)$ -entry of the block $V_{r,s}$ . In view of $V_{i,j}=S^{-T}_{n}\Theta_{i,j}S^{-1}_{n}$ and (14), we have that $W_{i,i}=U_{i}+U_{i+1}$ and $W_{i,i+1}=-U_{i+1}$ for $i=1,\ldots,n-1$ , and $W_{n,n}=U_{n}$ . Moreover, since $W=E^{T}VE=E^{T}\mathbf{S}^{-T}\Theta\,\mathbf{S}^{-1}E$ , $W$ is PD (resp. PSD) if and only if $\Theta$ is PD (resp. PSD), which is also equivalent to that each $U_{i}$ is PD (resp. PSD); see the comment after (15).

In what follows, we consider PSD (resp. PD) matrix decomposition for a PSD (resp. PD) $W$ generated by $\boldsymbol{\theta}_{i,j}\in\mathbb{R}^{n}_{+}$ for $i=1,\ldots,n$ and $j\geq i$ . The goal of this decomposition is to construct PSD matrices $\widetilde{W}^{s}\in\mathbb{R}^{np\times np}$ for $s=1,\ldots,n$ such that the following conditions hold:

(i)

\widetilde{W}^{1}\,=\,\begin{bmatrix}(\widetilde{W}^{1})_{1,1}&(\widetilde{W}^{1})_{1,2}&&&\\ (\widetilde{W}^{1})_{2,1}&(\widetilde{W}^{1})_{2,2}&&&\\ &&0&\cdots&0&\\ &&\vdots&\cdots&\vdots&\\ &&0&\cdots&0\end{bmatrix},\quad\widetilde{W}^{n}\,=\,\begin{bmatrix}0&\cdots&0&&&\\ \vdots&\cdots&\vdots&&&\\ 0&\cdots&0&&\\ &&&(\widetilde{W}^{n})_{n-1,n-1}&(\widetilde{W}^{n})_{n-1,n}\\ &&&(\widetilde{W}^{n})_{n,n-1}&(\widetilde{W}^{n})_{n,n}\\ \end{bmatrix};

(ii)

for each $s=2,\ldots,n-1$ ,

\widetilde{W}^{s}\,=\,\begin{bmatrix}{\mathbf{0}}_{(i-2)p\times(i-2)p}&&&&\\ &(\widetilde{W}^{s})_{s-1,s-1}&(\widetilde{W}^{s})_{s-1,s}&0&&\\ &(\widetilde{W}^{s})_{s,s-1}&(\widetilde{W}^{s})_{s,s}&(\widetilde{W}^{s})_{s,s+1}&&\\ &0&(\widetilde{W}^{s})_{s+1,s}&(\widetilde{W}^{s})_{s+1,s+1}&&&\\ &&&&0&\cdots&0&\\ &&&&\vdots&\cdots&\vdots&\\ &&&&0&\cdots&0\end{bmatrix};\ \ \ \mbox{and}

(iii)

$W=\sum^{n}_{s=1}\widetilde{W}^{s}$ .

For notational simplicity, let $\widehat{W}^{s}$ denote the possibly nonzero block in each $\widetilde{W}^{s}$ . Specifically,

\widehat{W}^{1}:=\begin{bmatrix}(\widetilde{W}^{1})_{1,1}&(\widetilde{W}^{1})_{1,2}\\ (\widetilde{W}^{1})_{2,1}&(\widetilde{W}^{1})_{2,2}\end{bmatrix}\in\mathbb{R}^{2p\times 2p},\qquad\widehat{W}^{n}:=\begin{bmatrix}(\widetilde{W}^{n})_{n-1,n-1}&(\widetilde{W}^{n})_{n-1,n}\\ (\widetilde{W}^{n})_{n,n-1}&(\widetilde{W}^{n})_{n,n}\end{bmatrix}\in\mathbb{R}^{2p\times 2p},

and for each $s=2,\ldots,n-1$ ,

\widehat{W}^{s}:=\begin{bmatrix}(\widetilde{W}^{s})_{s-1,s-1}&(\widetilde{W}^{s})_{s-1,s}&0\\ (\widetilde{W}^{s})_{s,s-1}&(\widetilde{W}^{s})_{s,s}&(\widetilde{W}^{s})_{s,s+1}\\ 0&(\widetilde{W}^{s})_{s+1,s}&(\widetilde{W}^{s})_{s+1,s+1}\end{bmatrix}\in\mathbb{R}^{3p\times 3p}.

When $W$ is PD, we also want each $\widehat{W}^{s}$ in the above decomposition to be PD.

Proposition 4.1.

Let $W$ be a PSD matrix generated by $\boldsymbol{\theta}_{i,j}\in\mathbb{R}^{n}_{+}$ for $i=1,\ldots,n$ . Then there exist PSD matrices $\widetilde{W}^{s},\ s=1,\ldots,n$ satisfying the above conditions. Moreover, suppose $W$ is PD. Then there exist PD matrices $\widehat{W}^{s},\ s=1,\ldots,n$ such that their corresponding $\widetilde{W}^{s}$ ’s satisfy the above conditions.

Proof.

Let $W$ be generated by $\boldsymbol{\theta}_{i,j}$ ’s such that $W$ is PSD, and let $U_{i}$ ’s be defined in (15) corresponding to $\boldsymbol{\theta}_{i,j}$ ’s. Note that each $U_{i}$ is PSD as $W$ is PSD. Let

\widetilde{W}^{1}\,=\,\begin{bmatrix}U_{1}&0&&&\\ 0&0&&&\\ &&0&\cdots&0\\ &&\vdots&\cdots&\vdots\\ &&0&\cdots&0\end{bmatrix},\qquad\widetilde{W}^{n}\,=\,\begin{bmatrix}0&\cdots&0&&\\ \vdots&\cdots&\vdots&&\\ 0&\cdots&0&&\\ &&&U_{n}&-U_{n}\\ &&&-U_{n}&U_{n}\\ \end{bmatrix},

and for each $s=2,\ldots,n-1$ ,

\widetilde{W}^{s}\,=\,\begin{bmatrix}{\mathbf{0}}_{(i-2)p\times(i-2)p}&&&&\\ &U_{s}&-U_{s}&0&&\\ &-U_{s}&U_{s}&0&&\\ &0&0&0&&&\\ &&&&0&\cdots&0&\\ &&&&\vdots&\cdots&\vdots&\\ &&&&0&\cdots&0\end{bmatrix}.

Since each $U_{i}$ is PSD, so is $\widetilde{W}^{s}$ for each $s=1,\ldots,n$ . Clearly, $W=\sum^{n}_{s=1}\widetilde{W}^{s}$ .

Now suppose $W$ is PD. Hence, each $U_{i}$ given by (15) is PD. Define

\breve{W}^{1}:=\frac{1}{2}\begin{bmatrix}U_{1}+U_{2}&-U_{2}\\ -U_{2}&U_{2}\end{bmatrix},\qquad\breve{W}^{2}:=\frac{1}{2}\begin{bmatrix}U_{1}+U_{2}&-U_{2}&0\\ -U_{2}&U_{2}+U_{3}&-U_{3}\\ 0&-U_{3}&U_{3}\end{bmatrix},\qquad\breve{W}^{n}:=\frac{1}{2}\begin{bmatrix}U_{n}&-U_{n}\\ -U_{n}&U_{n}\end{bmatrix},

and for each $s=3,\ldots,n-1$ ,

\breve{W}^{s}:=\frac{1}{2}\begin{bmatrix}U_{s}&-U_{s}&0\\ -U_{s}&U_{s}+U_{s+1}&-U_{s+1}\\ 0&-U_{s+1}&U_{s+1}\end{bmatrix}\in\mathbb{R}^{3p\times 3p}.

Note that $\breve{W}^{1}=\frac{1}{2}\left\{\begin{bmatrix}U_{1}&0\\ 0&0\end{bmatrix}+\begin{bmatrix}U_{2}&-U_{2}\\ -U_{2}&U_{2}\end{bmatrix}\right\}$ and the two matrices on the right hand side are both PSD and the intersection of their null spaces is the zero subspace. Hence, $\breve{W}^{1}$ is PD. Similarly, $\breve{W}^{2}$ is PD, and the other $\breve{W}^{s}$ ’s are PSD. Since $\breve{W}^{1}$ is PD, we see that for an arbitrary $\delta_{1}\in(0,\lambda_{\min}(\breve{W}^{1}))$ , $\widehat{W}^{1}:=\breve{W}^{1}-\delta_{1}\cdot I_{2p}$ is PD. Hence,

\grave{W}^{2}:=\breve{W}^{2}+\delta_{1}\cdot\begin{bmatrix}I_{p}&&\\ &I_{p}&\\ &&0\end{bmatrix}=\frac{1}{2}\begin{bmatrix}U_{1}+U_{2}+2\delta_{1}\cdot I_{p}&-U_{2}&0\\ -U_{2}&U_{2}+U_{3}+2\delta_{1}\cdot I_{p}&-U_{3}\\ 0&-U_{3}&U_{3}\end{bmatrix}

is also PD. Therefore, for an arbitrary $\delta_{2}\in(0,\lambda_{\min}(\grave{W}^{2}))$ , the matrix $\widehat{W}^{2}:=\grave{W}^{2}-\delta_{2}\cdot\begin{bmatrix}0&&\\ &I_{p}&\\ &&I_{p}\end{bmatrix}$ is PD. Further, it is easy to show that the matrix $\grave{W}^{3}:=\breve{W}^{3}+\delta_{2}\cdot\begin{bmatrix}I_{p}&&\\ &I_{p}&\\ &&0\end{bmatrix}$ is PD such that for any $\delta_{3}\in(0,\lambda_{\min}(\grave{W}^{3}))$ , the matrix $\widehat{W}^{3}:=\grave{W}^{3}-\delta_{3}\cdot\begin{bmatrix}0&&\\ &I_{p}&\\ &&I_{p}\end{bmatrix}$ is PD. Continuing this process by induction, we see that $\widehat{W}^{s}$ is PD for all $s=4,\ldots,n-1$ and $\widehat{W}^{n-1}:=\grave{W}^{n-1}-\delta_{n-1}\cdot\begin{bmatrix}0&&\\ &I_{p}&\\ &&I_{p}\end{bmatrix}$ is PD for an arbitrary $\delta_{n-1}\in(0,\lambda_{\min}(\grave{W}^{n-1}))$ , where $\grave{W}^{n-1}$ is PD. Finally, define $\widehat{W}^{n}:=\breve{W}^{n}+\delta_{n-1}\cdot I_{2p},$ which is clearly PD. Using these PD $\widehat{W}^{s},\,s=1,\ldots,n$ , we construct $\widetilde{W}^{s}$ by setting the possibly nonzero block in each $\widetilde{W}^{s}$ as $\widehat{W}^{s}$ . Specifically,

\begin{bmatrix}(\widetilde{W}^{1})_{1,1}&(\widetilde{W}^{1})_{1,2}\\ (\widetilde{W}^{1})_{2,1}&(\widetilde{W}^{1})_{2,2}\end{bmatrix}=\widehat{W}^{1}\in\mathbb{R}^{2p\times 2p},\qquad\begin{bmatrix}(\widetilde{W}^{n})_{n-1,n-1}&(\widetilde{W}^{n})_{n-1,n}\\ (\widetilde{W}^{n})_{n,n-1}&(\widetilde{W}^{n})_{n,n}\end{bmatrix}=\widehat{W}^{n}\in\mathbb{R}^{2p\times 2p},

and for each $s=2,\ldots,n-1$ ,

\begin{bmatrix}(\widetilde{W}^{s})_{s-1,s-1}&(\widetilde{W}^{s})_{s-1,s}&0\\ (\widetilde{W}^{s})_{s,s-1}&(\widetilde{W}^{s})_{s,s}&(\widetilde{W}^{s})_{s,s+1}\\ 0&(\widetilde{W}^{s})_{s+1,s}&(\widetilde{W}^{s})_{s+1,s+1}\end{bmatrix}=\widehat{W}^{s}\in\mathbb{R}^{3p\times 3p}.

A straightforward calculation shows that $W=\sum^{n}_{s=1}\widetilde{W}^{s}$ , yielding the desired result. ∎

To obtain a desired decomposition using the above proposition, we observe that the matrix $V$ in (8) is given by $\mathbf{S}^{-T}\Theta\mathbf{S}^{-1}$ for some matrix $\Theta$ of the form given below (14) whose blocks are positive combinations of $Q_{z,s},Q_{z^{\prime},s}$ and $Q_{w,s}$ . Since $Q_{z,s}$ and $Q_{z^{\prime},s}$ are diagonal and PSD and $Q_{w,s}$ are diagonal and PD, each block of $\Theta$ is diagonal and PD or PSD. Moreover, by Lemma 3.1, $W$ is PD. Hence, there are uncountably many ways to construct positive $\delta_{s}$ , and thus PD $\widehat{W}^{s}$ , as shown in the above proposition. Therefore, we obtain the following strongly convex decomposition for the objective function $J$ in (7), where we set the constant $\gamma=0$ without loss of generality:

J(\mathbf{u})=\frac{1}{2}{\mbox{\bf u}}^{T}W{\mbox{\bf u}}+c^{T}{\mbox{\bf u}}=\sum^{n}_{i=1}\frac{1}{2}{\mbox{\bf u}}^{T}\widetilde{W}^{i}{\mbox{\bf u}}+\sum^{n}_{i=1}c^{T}_{\mathcal{I}_{i}}{\mbox{\bf u}}_{i}=J_{1}({\mbox{\bf u}}_{1},{\mbox{\bf u}}_{2})+\sum^{n-1}_{i=2}J_{i}({\mbox{\bf u}}_{i-1},{\mbox{\bf u}}_{i},{\mbox{\bf u}}_{i+1})+J_{n}({\mbox{\bf u}}_{n-1},{\mbox{\bf u}}_{n})

where the strongly convex functions $J_{i}$ are given by

$\displaystyle J_{1}({\mbox{\bf u}}_{1},{\mbox{\bf u}}_{2})$	$\displaystyle:=$	$\displaystyle\frac{1}{2}\begin{bmatrix}{\mbox{\bf u}}^{T}_{1}&{\mbox{\bf u}}^{T}_{2}\end{bmatrix}\widehat{W}_{1}\begin{bmatrix}{\mbox{\bf u}}_{1}\\ {\mbox{\bf u}}_{2}\end{bmatrix}+c^{T}_{\mathcal{I}_{1}}{\mbox{\bf u}}_{1},$
$\displaystyle J_{i}({\mbox{\bf u}}_{i-1},{\mbox{\bf u}}_{i},{\mbox{\bf u}}_{i+1})$	$\displaystyle:=$	$\displaystyle\frac{1}{2}\begin{bmatrix}{\mbox{\bf u}}^{T}_{i-1}&{\mbox{\bf u}}^{T}_{i}&{\mbox{\bf u}}^{T}_{i+1}\end{bmatrix}\widehat{W}_{i}\begin{bmatrix}{\mbox{\bf u}}_{i-1}\\ {\mbox{\bf u}}_{i}\\ {\mbox{\bf u}}_{i+1}\end{bmatrix}+c^{T}_{\mathcal{I}_{i}}{\mbox{\bf u}}_{i},\qquad\forall\ i=2,\ldots,n-1,$	(16)
$\displaystyle J_{n}({\mbox{\bf u}}_{n-1},{\mbox{\bf u}}_{n})$	$\displaystyle:=$	$\displaystyle\frac{1}{2}\begin{bmatrix}{\mbox{\bf u}}^{T}_{n-1}&{\mbox{\bf u}}^{T}_{n}\end{bmatrix}\widehat{W}_{n}\begin{bmatrix}{\mbox{\bf u}}_{n-1}\\ {\mbox{\bf u}}_{n}\end{bmatrix}+c^{T}_{\mathcal{I}_{n}}{\mbox{\bf u}}_{n}.$

Remark 4.1.

Clearly, the above decomposition method is applicable to any vehicle communication network satisfying the assumption $\bf A.1$ in Section 2, i.e., $(i,i+1)\in\mathcal{E}$ for all $i=1,\ldots,n-1$ . Besides, various alternative approaches can be developed to construct PD matrices $\widehat{W}^{s}$ using the similar idea given in the above proposition. Further, a similar decomposition method can be developed for other vehicle communication networks different from the the cyclic-like graph. For instance, if such a graph contains edges other than $(i,i+1)\in\mathcal{E}$ , one can add or subtract certain small terms pertaining to these extra edges in relevant matrices, which will preserve the desired PD property.

In what follows, we write each $J_{i}$ as $J_{i}({\mbox{\bf u}}_{i},(\mathbf{u}_{j})_{j\in\mathcal{N}_{i}})$ for notational convenience, where $\mathcal{N}_{i}$ denotes the set of neighbors of vehicle $i$ in a vehicle network such that $i-1,i+1\in\mathcal{N}_{i}$ for $i=2,\ldots,n-1$ and $2\in\mathcal{N}_{1},n-1\in\mathcal{N}_{n}$ .

4.3 Operator Splitting Method based Fully Distributed Algorithms

For illustration simplicity, we consider the cyclic like network topology through this subsection, although the proposed schemes can be easily extended to other network topologies under a suitable assumption ((cf. Remark 4.1). In this case, $\mathcal{N}_{1}=\{2\}$ , $\mathcal{N}_{n}=\{n-1\}$ , and $\mathcal{N}_{i}=\{i-1,i+1\}$ for $i=2,\ldots,n-1$ .

Define the constraint set

\displaystyle\mathcal{P}

\displaystyle\,:=\,\{\mathbf{u}=(\mathbf{u}_{1},\ldots,\mathbf{u}_{n})\in\mathbb{R}^{np}\,|\,\mathbf{u}_{i}\in\mathcal{X}_{i},\ \mathbf{u}\in\mathcal{Z}_{i},\ i=1,\ldots,n\}.

Recall that $\mathcal{P}$ is defined by convex quadratic functions. The underlying optimization problem (7) at time $k$ becomes $\min_{\mathbf{u}}J(\mathbf{u})$ subject to $\mathbf{u}\in\mathcal{P}$ .

We formulate this problem as a locally coupled convex optimization problem [6] and solve it using fully distributed algorithms. Specifically, in view of the decompositions given by (16), the objective function in (7) can be written as

J(\mathbf{u})=\sum^{n}_{i=1}J_{i}(\mathbf{u}_{i},(\mathbf{u}_{j})_{j\in\mathcal{N}_{i}}),

In view of Remark 3.2, the safety constraints are also locally coupled. Let $\boldsymbol{\delta}_{S}$ denote the indicator function of a (closed convex) set $S$ . Define, for each $i=1,\ldots,n$ ,

\displaystyle\widehat{J}_{i}(\mathbf{u}_{i},(\mathbf{u}_{j})_{j\in\mathcal{N}_{i}})

\displaystyle\,:=\,J_{i}(\mathbf{u}_{i},(\mathbf{u}_{j})_{j\in\mathcal{N}_{i}})+\boldsymbol{\delta}_{\mathcal{X}_{i}}({\mbox{\bf u}}_{i})+\boldsymbol{\delta}_{\mathcal{Z}_{i}}(\mathbf{u}_{i-1},\mathbf{u}_{i}).

As in [6], define $\widehat{\mathbf{u}}_{i}:=\big{(}\mathbf{u}_{i},(\mathbf{u}_{i,j})_{j\in\mathcal{N}_{i}}\big{)}$ , where the new variables $\mathbf{u}_{i,j}$ represent the predicted values of $\mathbf{u}_{j}$ of vehicle $j$ in the neighbor of vehicle $i$ , and let $\widehat{\mathbf{u}}:=(\widehat{\mathbf{u}}_{i})_{i=1,\ldots,n}\in\mathbb{R}^{\ell}$ . Define the consensus subspace

\mathcal{A}\,:=\,\{\widehat{\mathbf{u}}\,|\,\mathbf{u}_{i,j}=\mathbf{u}_{j},\ \forall\,(i,j)\in\mathcal{E}\}.

Then the underlying optimization problem (7) can be equivalently written as

\min_{\widehat{\mathbf{u}}}\ \ \sum^{n}_{i=1}\widehat{J}_{i}(\widehat{\mathbf{u}}_{i}),\quad\mbox{ subject to }\quad\widehat{\mathbf{u}}\in\mathcal{A}.

Let $\mathcal{P}_{i}:=\{\widehat{\mathbf{u}}_{i}\,|\,\mathbf{u}_{i}\in\mathcal{X}_{i},\ (H_{i}(\mathbf{u}_{i,i-1},\mathbf{u}_{i}))_{s}\leq 0,\,\forall\,s=1,\ldots,p\}$ for notational simplicity. Then the underlying optimization problem becomes

\min_{\widehat{\mathbf{u}}}\,F(\widehat{\mbox{\bf u}}):=\sum^{n}_{i=1}J_{i}(\widehat{\mathbf{u}}_{i})+\sum^{n}_{i=1}\boldsymbol{\delta}_{\mathcal{P}_{i}}(\widehat{\mathbf{u}}_{i})+\boldsymbol{\delta}_{\mathcal{A}}(\widehat{\mathbf{u}}),

(17)

where $F:\mathbb{R}^{\ell}\rightarrow\mathbb{R}\cup\{+\infty\}$ denotes the extended-valued objective function. Thus $F$ is the sum of two indictor functions of closed convex sets and the convex quadratic function given by $J(\widehat{\mathbf{u}}):=\sum^{n}_{i=1}J_{i}(\widehat{\mathbf{u}}_{i})$ , by slightly abusing the notation. Note that $\mathcal{A}$ is polyhedral. It is easy to show via Corollary 3.1 that the Slater’s condition holds under the mild assumptions given in Corollary 3.1, e.g., $v_{0}(k)>v_{\min}$ for all $k\in\mathbb{Z}_{+}$ . Hence, by [13, Corollary 23.8.1], $\partial F(\widehat{\mbox{\bf u}})=\sum^{n}_{i=1}\Big{(}\nabla J_{i}(\widehat{\mathbf{u}}_{i})+\mathcal{N}_{\mathcal{P}_{i}}(\widehat{\mathbf{u}}_{i})\Big{)}+\mathcal{N}_{\mathcal{A}}(\widehat{\mathbf{u}})$ in light of $\partial\boldsymbol{\delta}_{C}(x)=\mathcal{N}_{C}(x)$ , where $\mathcal{N}_{C}(x)$ denotes the normal cone of a closed convex set $C$ at $x\in C$ . Finally, the formulation given by (17) is a locally coupled convex optimization problem; see Section 4.1. This formulation allows one to develop fully distributed schemes. Particularly, in fully distributed computation, each vehicle $i$ only knows $\widehat{\mathbf{u}}_{i}$ and $\widehat{J}_{i}$ (i.e., $J_{i}$ and $\mathcal{P}_{i}$ ) but does not know $\widehat{\mathbf{u}}_{j}$ and $\widehat{J}_{j}$ with $j\neq i$ . Each vehicle $i$ will exchange information with its neighboring vehicles to update $\widehat{\mathbf{u}}_{i}$ via a distributed scheme.

Algorithm 1 Generalized Douglas-Rachford Splitting Method based Fully Distributed Algorithm

1: Choose constants

0<\alpha<1

and

\rho>0

2: Initialize

k=0

, and choose an initial point

z^{0}

3: while the stopping criteria is not met do

4: for

i=1,\ldots,n

5: Compute

\overline{z}^{k}_{i}

using equation (18), and let

w^{k+1}_{i}\leftarrow\overline{z}^{k}_{i}

6: end for

7: for

i=1,\ldots,n

z^{k+1}_{i}\leftarrow z^{k}_{i}+2\alpha\cdot\Big{[}\mbox{Prox}_{\rho\widehat{J}_{i}}\big{(}2w^{k+1}_{i}-z^{k}_{i}\big{)}-w^{k+1}_{i}\Big{]}

9: end for

10:

k\leftarrow k+1

11: end while

12: return

\widehat{\mathbf{u}}^{*}=w^{k}

We introduce more notation first. For a proper, lower semicontinuous convex function $f:\mathbb{R}^{n}\rightarrow\mathbb{R}\cup\{+\infty\}$ , let $\mbox{Prox}_{f}(\cdot)$ denote the proximal operator, i.e., for any given $x\in\mathbb{R}^{n}$ ,

\mbox{Prox}_{f}(x)\,:=\,\operatornamewithlimits{\arg\min}_{z\in\mathbb{R}^{n}}\ f(z)+\frac{1}{2}\|z-x\|^{2}_{2}.

Further, $\Pi_{C}$ denotes the Euclidean projection onto a closed convex set $C$ . Using this notation, we present a specific operator splitting method based distributed scheme for solving (17). By grouping the first two sums (with separable variables) in the objective function of (17), we apply the generalized Douglas-Rachford splitting algorithm [6]. Recall that $\widehat{J}_{i}(\widehat{\mathbf{u}}_{i}):=J_{i}(\widehat{\mathbf{u}}_{i})+\boldsymbol{\delta}_{\mathcal{P}_{i}}(\widehat{\mathbf{u}}_{i})$ for each $i=1,\ldots,n$ . For any constants $\alpha$ and $\rho$ satisfying $0<\alpha<1$ and $\rho>0$ , this algorithm is given by:

\displaystyle w^{k+1}\,=\,\Pi_{\mathcal{A}}(z^{k}),\qquad z^{k+1}\,=\,z^{k}+2\alpha\cdot\Big{[}\mbox{Prox}_{\rho\widehat{J}_{1}+\cdots+\rho\widehat{J}_{n}}\big{(}2w^{k+1}-z^{k}\big{)}-w^{k+1}\Big{]}.

It is shown in [3, 6] that the sequence $(w^{k})$ converges to the unique minimizer $\mathbf{\widehat{u}}^{*}$ of the underlying optimization problem (17). In the above scheme, $\Pi_{\mathcal{A}}$ is the orthogonal projection onto the consensus subspace $\mathcal{A}$ such that the following holds: for any $\widehat{\mathbf{u}}:=(\widehat{\mathbf{u}}_{1},\ldots,\widehat{\mathbf{u}}_{n})$ where $\widehat{\mathbf{u}}_{i}:=\big{(}\mathbf{u}_{i},(\mathbf{u}_{ij})_{j\in\mathcal{N}_{i}}\big{)}$ , $\overline{\mathbf{u}}:=\Pi_{\mathcal{A}}(\widehat{\mathbf{u}})$ is given by [6, Section IV]:

\overline{\mathbf{u}}_{j}=\overline{\mathbf{u}}_{ij}=\frac{1}{1+|\mathcal{N}_{j}|}\Big{(}\widehat{\mathbf{u}}_{j}+\sum_{k\in\mathcal{N}_{j}}\widehat{\mathbf{u}}_{kj}\Big{)},\qquad\forall\,(i,j)\in\mathcal{E}.

(18)

Furthermore, due to the decoupled structure of $\widehat{J}_{i}$ ’s, we obtain the distributed version of the above algorithm, which is also summarized in Algorithm 1:


$\displaystyle w^{k+1}_{i}$	$\displaystyle\,=\,\overline{z}^{k}_{i},\qquad i=1,\ldots,n;$	(19a)
$\displaystyle z^{k+1}_{i}$	$\displaystyle\,=\,z^{k}_{i}+2\alpha\cdot\Big{[}\mbox{Prox}_{\rho\widehat{J}_{i}}\big{(}2w^{k+1}_{i}-z^{k}_{i}\big{)}-w^{k+1}_{i}\Big{]},\quad i=1,\ldots,n.$	(19b)

In the distributed scheme (19), the first step (or Line 5 of Algorithms 1) is a consensus step, and the consensus computation is carried out in a fully distributed and synchronous manner as indicated in [6, Section IV]. The second step in (19) (or Line 8 of Algorithms 1) does not need inter-agent communication [6] and is performed using local computation only. For effective computation in the second step, recall that $\widehat{J}_{i}(\widehat{\mathbf{u}}_{i}):=J_{i}(\widehat{\mathbf{u}}_{i})+\boldsymbol{\delta}_{\mathcal{P}_{i}}(\widehat{\mathbf{u}}_{i})$ for each $i=1,\ldots,n$ such that the proximal operator $\mbox{Prox}_{\rho\widehat{J}_{i}}(\widehat{\mathbf{u}}_{i})$ becomes

\mbox{Prox}_{\rho\widehat{J}_{i}}(\widehat{\mathbf{u}}_{i})\,=\,\operatornamewithlimits{\arg\min}_{z\in\mathcal{P}_{i}}\ J_{i}(z)+\frac{1}{2\rho}\|z-\widehat{\mathbf{u}}_{i}\|^{2}_{2}.

Since $\mathcal{P}_{i}$ is the intersection of a polyhedral set and a quadratically constrained convex set and $J_{i}$ is a quadratic convex function, $\mbox{Prox}_{\rho\widehat{J}_{i}}(\widehat{\mathbf{u}}_{i})$ is formulated as a QCQP and can be solved via a second order cone program [1] or a semidefinite program. Efficient numerical packages, e.g., SeDuMi [16], can be used for solving the QCQP. Lastly, a typical (global) stopping criterion in the scheme (19) (or Algorithm 1) is defined by the error tolerance of two neighboring $z^{k}$ ’s, i.e., $\|z^{k+1}-z^{k}\|_{2}\leq\varepsilon$ , where $\varepsilon>0$ is an error tolerance. For distributed computation, one can use its local version, i.e., $\|z^{k+1}_{i}-z^{k}_{i}\|_{2}\leq\varepsilon/n$ , as a stopping criterion for each vehicle.

Remark 4.2.

Other distributed algorithms can be used to solve the underlying optimization problem (17). For example, the three operator splitting schemes developed in [3] can be applied. To describe such schemes, let $\widehat{L}:=\max_{i=1,\ldots,n}\|\widehat{W}^{i}\|_{2}>0$ . the Hessian $HJ(\widehat{\mathbf{u}})=\mbox{diag}(\widehat{W}^{i})_{i=1,\ldots,n}$ . Hence, $\nabla J$ is $\widehat{L}$ -Lipschitz continuous and thus $1/L$ -cocoercive. Further, the two indicator functions are proper, closed, and convex functions. Choose the constants $\gamma,\lambda$ such that $0<\gamma<2/\widehat{L}$ and $0<\lambda<2-\frac{\gamma\widehat{L}}{2}$ . Then for any initial condition $z^{0}$ , the iterative scheme is given by [3, Algorithm 1]:

\displaystyle w^{k+1}\,=\,\Pi_{\mathcal{A}}(z^{k}),\qquad z^{k+1}\,=\,z^{k}+\lambda\cdot\Big{[}\Pi_{\mathcal{P}_{1}\times\cdots\times\mathcal{P}_{n}}\big{(}2w^{k+1}-z^{k}-\gamma\nabla J(w^{k+1})\big{)}-w^{k+1}\Big{]}.

In view of the similar discussions for consensus computation and decoupled structure of the projection $\Pi_{\mathcal{P}_{1}\times\cdots\times\mathcal{P}_{n}}$ , we obtain the distributed version of the above algorithm:

	$\displaystyle w^{k+1}_{i}$	$\displaystyle\,=\,\overline{z}^{k}_{i},\qquad i=1,\ldots,n;$
	$\displaystyle z^{k+1}_{i}$	$\displaystyle\,=\,z^{k}_{i}+\lambda\cdot\Big{[}\Pi_{\mathcal{P}_{i}}\big{(}2w^{k+1}_{i}-z^{k}_{i}-\gamma\cdot[\widehat{W}^{i}w^{k+1}_{i}+c_{\mathcal{I}_{i}}]\big{)}-w^{k+1}_{i}\Big{]},\quad i=1,\ldots,n.$

In this scheme, the Euclidean projection $\Pi_{\mathcal{P}_{q,i}}$ can be formulated as a QCQP and be solved via a second order cone program.

When each $\widehat{W}^{i}$ is PD, each $J_{i}$ is strongly convex. Thus $\nabla J$ is $\mu$ -strongly monotone with $\mu=\min_{i=1,\ldots,n}\lambda_{\min}(\widehat{W}^{i})$ , i.e., $(x-y)^{T}(\nabla J(x)-\nabla J(y))\geq\mu\|x-y\|^{2}_{2},\forall\,x,y$ . Since the subdifferential of the indicator function of a closed convex set is monotone, an accelerated scheme developed in [3, Algorithm 2] can be exploited. In particular, let $\eta$ be a constant with $0<\eta<1$ , and $\gamma_{0}\in(0,2/(\widehat{L}\cdot(1-\eta))$ . Set the initial points for an arbitrary $z^{0}$ , $w^{0}=\Pi_{\mathcal{A}}(z^{0})$ and $v^{0}=(z^{0}-w^{0})/\gamma_{0}$ . The distributed version of this scheme is given by:

	$\displaystyle w^{k+1}_{i}$	$\displaystyle\,=\,\overline{z}^{k}_{i}+\gamma_{k}\overline{v}^{k}_{i},\qquad i=1,\ldots,n;$
	$\displaystyle v^{k+1}_{i}$	$\displaystyle\,=\,\frac{1}{\gamma_{k}}\big{(}z^{k}_{i}+\gamma_{k}v^{k}_{i}-w^{k+1}_{i}\big{)},\qquad i=1,\ldots,n;$
	$\displaystyle\gamma_{k+1}$	$\displaystyle\,=\,-\widetilde{\mu}\gamma^{2}_{k}+\sqrt{(\widetilde{\mu}\gamma^{2}_{k})^{2}+\gamma^{2}_{k}};$
	$\displaystyle z^{k+1}_{i}$	$\displaystyle\,=\Pi_{\mathcal{P}_{i}}\Big{(}w^{k+1}_{i}-\gamma_{k+1}v^{k+1}_{i}-\gamma_{k+1}[\widehat{W}^{i}w^{k+1}_{i}+c_{\mathcal{I}_{i}}]\Big{)},\quad i=1,\ldots,n,$

where $\widetilde{\mu}:=\eta\cdot\mu$ . It is shown in [3, Theorem 1.2] that $(w^{k})$ converges to the unique minimizer $\mathbf{\widehat{u}}^{*}$ with $\|w^{k}-\widehat{\mathbf{u}}^{*}\|_{2}=O(1/(k+1))$ .

5 Control Design and Stability Analysis of the Closed Loop Dynamics

In this section, we discuss how to choose the weight matrices $Q_{z,s}$ , $Q_{z^{\prime},s}$ and $Q_{w,s}$ to achieve the desired closed loop performance, including stability and traffic transient dynamics. For the similar reasons given in [5, Section 5], we focus on the constraint free case.

Under the linear vehicle dynamics, the closed-loop system is also a linear system. Specifically, the linear closed-loop dynamics are given by

z(k+1)\,=\,z(k)+\tau z^{\prime}(k)+\frac{\tau^{2}}{2}w(k),\qquad z^{\prime}(k+1)\,=\,z^{\prime}(k)+\tau w(k),

(20)

where $w(k)$ is a unique solution to an unconstrained optimization problem arising from the MPC and is a linear function of $z(k)$ and $z^{\prime}(k)$ to be determined as follows.

Case (i): $p=1$ . In this case, we write $Q_{z,1},Q_{z^{\prime},1},Q_{w,1}$ as $Q_{z},Q_{z^{\prime}},Q_{w}$ respectively. Then the objective function becomes

J(w(k))\,=\,\frac{1}{2}\Big{[}z^{T}(k+1)Q_{z},z(k+1)+(z^{\prime}(k+1))^{T}Q_{z^{\prime}}z^{\prime}(k+1)\Big{]}+\frac{\tau^{2}}{2}\widetilde{w}^{T}(k)Q_{w}\widetilde{w}(k),

where we recall that $\widetilde{w}(k)=w(k)-u_{0}(k)\mathbf{e}_{1}$ . It follows from the similar argument in [5, Section 5] that the the closed-loop system is given by the following linear system:

\displaystyle\begin{bmatrix}z(k+1)\\ z^{\prime}(k+1)\end{bmatrix}\,=\,\underbrace{\begin{bmatrix}I_{n}-\frac{\tau^{2}}{4}\widehat{W}Q_{z}&\tau I_{n}-\widehat{W}\Big{(}\frac{\tau^{3}}{4}Q_{z}+\frac{\tau}{2}Q_{z^{\prime}}\Big{)}\\ -\frac{\tau}{2}\widehat{W}Q_{z}&I_{n}-\widehat{W}\Big{(}\frac{\tau^{2}}{2}Q_{z}+Q_{z^{\prime}}\Big{)}\end{bmatrix}}_{A_{\mbox{c}}}\begin{bmatrix}z(k)\\ z^{\prime}(k)\end{bmatrix}+\begin{bmatrix}\frac{\tau^{2}}{2}I_{n}\\ \tau I_{n}\end{bmatrix}\widehat{W}Q_{w}\mathbf{e}_{1}\cdot u_{0}(k),

(21)

where $A_{\mbox{c}}$ is the closed loop dynamics matrix, and

\widehat{W}\,:=\,\left[\frac{\tau^{2}Q_{z}}{4}+Q_{z^{\prime}}+Q_{w}\right]^{-1}.

(22)

The matrix $A_{\mbox{c}}$ in (21) plays an important role in the closed loop stability and desired transient dynamical performance. Since $Q_{z},Q_{z^{\prime}}$ and $Q_{w}$ are all diagonal and PSD (resp. PD), we have $Q_{z}=\mbox{diag}(\boldsymbol{\alpha})$ , $Q_{z^{\prime}}=\mbox{diag}(\boldsymbol{\beta})$ , and $Q_{w}=\mbox{diag}(\boldsymbol{\zeta})$ , where $\boldsymbol{\alpha},\boldsymbol{\beta}\in\mathbb{R}^{n}_{+}$ and $\boldsymbol{\zeta}\in\mathbb{R}^{n}_{++}$ with $\boldsymbol{\alpha}=(\alpha_{i})^{n}_{i=1}$ , $\boldsymbol{\beta}=(\beta_{i})^{n}_{i=1}$ , and $\boldsymbol{\zeta}=(\zeta_{i})^{n}_{i=1}$ . Hence, we write the matrix $A_{\mbox{c}}$ as $A_{\mbox{c}}(\boldsymbol{\alpha},\boldsymbol{\beta},\boldsymbol{\zeta},\tau)$ to emphasize its dependence on these parameters. The following result asserts the asymptotic stability of the linear closed-loop dynamics; its proof resembles that for [5, Proposition 5.1] and is thus omitted.

Proposition 5.1.

Given any $\tau\in\mathbb{R}_{++}$ and any $\boldsymbol{\alpha},\boldsymbol{\beta},\boldsymbol{\zeta}\in\mathbb{R}^{n}_{++}$ , the matrix $A_{\mbox{c}}(\boldsymbol{\alpha},\boldsymbol{\beta},\boldsymbol{\zeta},\tau)$ is Schur stable, i.e., each eigenvalue $\mu\in\mathbb{C}$ of $A_{\mbox{c}}(\boldsymbol{\alpha},\boldsymbol{\beta},\boldsymbol{\zeta},\tau)$ satisfies $|\mu|<1$ . Moreover, for any eigenvalue $\mu_{i}$ of $A_{\mbox{c}}$ , the following hold:

(1)

if $\mu_{i}$ is non-real, then $|\mu_{i}|^{2}=\frac{\zeta_{i}}{d_{i}}$ , where $d_{i}:=\frac{\alpha_{i}\tau^{2}}{4}+\beta_{i}+\zeta_{i}$ ;
(2)

if $\mu_{i}$ is real, then $1-(\frac{\alpha_{i}\tau^{2}}{2}+\beta_{i})\frac{1}{d_{i}}<\mu_{i}<1-\frac{\alpha_{i}\tau^{2}}{4d_{i}}$ .

Case (ii): $p>1$ . Fix $k$ . For a general $p\in\mathbb{N}$ , let $\mathbf{w}:=(w(k),\ldots,w(k+p-1))\in\mathbb{R}^{np}$ . Recall that for each $s=1,\ldots,p$ ,

z(k+s)\,=\,z(k)+s\tau z^{\prime}(k)+\tau^{2}\sum^{s-1}_{j=0}\frac{2(s-j)-1}{2}w(k+j),\quad z^{\prime}(k+s)=z^{\prime}(k)+\tau\sum^{s-1}_{j=0}w(k+j),

and with slightly abusing notation, the objective function is

	$\displaystyle J(\underbrace{w(k),\ldots,w(k+p-1)}_{\mathbf{w}})$	$\displaystyle=$	$\displaystyle\frac{1}{2}\sum^{p}_{s=1}\Big{(}\tau^{2}\widetilde{w}^{T}(k+s-1)Q_{w,s}\widetilde{w}(k+s-1)$
			$\displaystyle\qquad\quad+\ z^{T}(k+s)Q_{z,s}z(k+s)\,+\,(z^{\prime}(k+s))^{T}Q_{z^{\prime},s}z^{\prime}(k+s)\Big{)},$

where $\widetilde{w}(k+s):=w(k+s)-u_{0}(k)\cdot\mathbf{e}_{1}$ introduced in Remark 3.1. It follows from the similar development in Section 3.1 that

J(\mathbf{w})=\frac{1}{2}\mathbf{w}^{T}\mathbf{H}\mathbf{w}+\mathbf{w}^{T}\left(\mathbf{G}\begin{bmatrix}z(k)\\ z^{\prime}(k)\end{bmatrix}-u_{0}(k)\mathbf{g}\right)+\widetilde{\gamma},

where $\widetilde{\gamma}$ is a constant. By a similar argument as in Lemma 3.1, it can be shown that $\mathbf{H}\in\mathbb{R}^{pn\times pn}$ is a symmetric PD matrix. Further, it resembles the matrix $V$ in (8) (by replacing $S^{-1}_{n}$ with $I_{n}$ ), i.e.,

\mathbf{H}\,=\,\begin{bmatrix}\breve{\mbox{\bf H}}_{1,1}+\tau^{2}Q_{w,1}&\breve{\mbox{\bf H}}_{1,2}&\breve{\mbox{\bf H}}_{1,3}&\cdots&\cdots&\breve{\mbox{\bf H}}_{1,p}\\ \breve{\mbox{\bf H}}_{2,1}&\breve{\mbox{\bf H}}_{2,2}+\tau^{2}Q_{w,2}&\breve{\mbox{\bf H}}_{2,3}&\cdots&\cdots&\breve{\mbox{\bf H}}_{2,p}\\ \cdots&&\cdots&&\cdots&\\ \cdots&&\cdots&&\cdots&\\ \breve{\mbox{\bf H}}_{p,1}&\breve{\mbox{\bf H}}_{p,2}&\breve{\mbox{\bf H}}_{p,3}&\cdots&\cdots&\breve{\mbox{\bf H}}_{p,p}+\tau^{2}Q_{w,p}\end{bmatrix}\in\mathbb{R}^{np\times np},

(23)

where $\breve{\mbox{\bf H}}_{i,j}$ ’s are diagonal PD matrices given by

\breve{\mbox{\bf H}}_{i,j}\,:=\,\sum^{p}_{s=\max(i,j)}\Big{(}\frac{\tau^{4}}{4}[2(s-i)+1]\cdot[2(s-j)+1]Q_{z,s}+\tau^{2}Q_{z^{\prime},s}\Big{)}\in\mathbb{R}^{n\times n}.

Moreover, it follows from (10) and (11) that the matrix $\mathbf{G}$ and constant vector $\mathbf{g}$ are

\mathbf{G}\,:=\,\begin{bmatrix}\mathbf{G}_{1,1}&\mathbf{G}_{1,2}\\ \vdots&\vdots\\ \mathbf{G}_{p,1}&\mathbf{G}_{p,2}\end{bmatrix}\in\mathbb{R}^{pn\times 2n},\qquad\mathbf{g}\,:=\,\tau^{2}\begin{bmatrix}Q_{w,1}\mathbf{e}_{1}\\ \vdots\\ Q_{w,p}\mathbf{e}_{1}\end{bmatrix}\in\mathbb{R}^{pn}.

(24)

where $\mathbf{G}_{i,1},\mathbf{G}_{i,2}\in\mathbb{R}^{n\times n}$ are given by: for each $i=1,\ldots,p$ ,

\displaystyle\mathbf{G}_{i,1}

\displaystyle=

\displaystyle\tau^{2}\sum^{p}_{s=i}\frac{2(s-i)+1}{2}Q_{z,s},\qquad\quad\mathbf{G}_{i,2}\,=\,\tau^{3}\sum^{p}_{s=i}s\frac{2(s-i)+1}{2}Q_{z,s}+\tau\sum^{p}_{s=i}Q_{z^{\prime},s}.

Hence, the optimal solution is $\mathbf{w}_{*}=(w_{*}(k),w_{*}(k+1),\ldots,w_{*}(k+p-1))=-\mathbf{H}^{-1}\Big{(}\mathbf{G}\begin{bmatrix}z(k)\\ z^{\prime}(k)\end{bmatrix}-u_{0}(k)\mathbf{g}\Big{)}$ , and $w_{*}(k)=-\begin{bmatrix}I_{n}&0&\cdots&0\end{bmatrix}\mathbf{H}^{-1}\Big{(}\mathbf{G}\begin{bmatrix}z(k)\\ z^{\prime}(k)\end{bmatrix}-u_{0}(k)\mathbf{g}\Big{)}$ . Define the matrix $\mathbf{K}$ and the vector $\mathbf{d}$ as

\mathbf{K}:=-\begin{bmatrix}I_{n}&0&\cdots&0\end{bmatrix}\mathbf{H}^{-1}\mathbf{G}\in\mathbb{R}^{n\times 2n},\qquad\mathbf{d}:=\begin{bmatrix}I_{n}&0&\cdots&0\end{bmatrix}\mathbf{H}^{-1}\mathbf{g}\in\mathbb{R}^{n}.

(25)

The closed loop system becomes

\begin{bmatrix}z(k+1)\\ z^{\prime}(k+1)\end{bmatrix}=\underbrace{\left\{\begin{bmatrix}I_{n}&\tau I_{n}\\ 0&I_{n}\end{bmatrix}+\begin{bmatrix}\frac{\tau^{2}}{2}I_{n}\\ \tau I_{n}\end{bmatrix}\mathbf{K}\right\}}_{A_{\mbox{c}}}\begin{bmatrix}z(k)\\ z^{\prime}(k)\end{bmatrix}+\begin{bmatrix}\frac{\tau^{2}}{2}I_{n}\\ \tau I_{n}\end{bmatrix}u_{0}(k)\cdot\mathbf{d},

where $A_{\mbox{c}}$ is the closed loop dynamics matrix, and the subscript of $A_{\mbox{c}}$ represents the closed loop.

Since $Q_{z,s},Q_{z^{\prime},s}$ are diagonal PSD and $Q_{w,s}$ are diagonal PD for all $s=1,\ldots,p$ , we write them as $Q_{z,s}=\mbox{diag}(\boldsymbol{\alpha}^{s})$ , $Q_{z^{\prime},s}=\mbox{diag}(\boldsymbol{\beta}^{s})$ , and $Q_{w,s}=\mbox{diag}(\boldsymbol{\zeta}^{s})$ , where $\boldsymbol{\alpha}^{s},\boldsymbol{\beta}^{s}\in\mathbb{R}^{n}_{+}$ and $\boldsymbol{\zeta}^{s}\in\mathbb{R}^{n}_{++}$ for all $s=1,\ldots,p$ with $\boldsymbol{\alpha}^{s}=(\alpha^{s}_{i})^{n}_{i=1}$ , $\boldsymbol{\beta}^{s}=(\beta^{s}_{i})^{n}_{i=1}$ , and $\boldsymbol{\zeta}^{s}=(\zeta^{s}_{i})^{n}_{i=1}$ for each $s$ . Let $\boldsymbol{\alpha}:=(\boldsymbol{\alpha}^{1},\ldots,\boldsymbol{\alpha}^{p})$ , $\boldsymbol{\beta}:=(\boldsymbol{\beta}^{1},\ldots,\boldsymbol{\beta}^{p})$ , and $\boldsymbol{\zeta}:=(\boldsymbol{\zeta}^{1},\ldots,\boldsymbol{\zeta}^{p})$ . We write the matrix $A_{\mbox{c}}$ as $A_{\mbox{c}}(\boldsymbol{\alpha},\boldsymbol{\beta},\boldsymbol{\zeta},\tau)$ to emphasize its dependence on these parameters.

It can be shown that there exists a permutation matrix $\widehat{E}\in\mathbb{R}^{2n\times 2n}$ such that $\widetilde{A}:=\widehat{E}^{T}A_{\mbox{c}}\widehat{E}$ is a block diagonal matrix, i.e., $\widetilde{A}=\mbox{diag}(\widetilde{A}_{1},\widetilde{A}_{2},\ldots,\widetilde{A}_{n})$ whose each block $\widetilde{A}_{i}\in\mathbb{R}^{2\times 2}$ is given by

\widetilde{A}_{i}\,=\,\begin{bmatrix}1&\tau\\ 0&1\end{bmatrix}+\begin{bmatrix}\frac{\tau^{2}}{2}\\ \tau\end{bmatrix}\widetilde{\mathbf{K}}_{i},\qquad\forall\ i=1,\ldots,n.

Here for each $i=1,\ldots,n$ , $\widetilde{\mathbf{K}}_{i}:=-\mathbf{e}^{T}_{1}\,\widetilde{\mathbf{H}}^{-1}_{i}\,\widetilde{\mathbf{G}}_{i}\in\mathbb{R}^{1\times 2}$ , where

\widetilde{\mbox{\bf H}}_{i}\,:=\,\begin{bmatrix}(\breve{\mbox{\bf H}}_{1,1}+\tau^{2}Q_{w,1})_{i,i}&(\breve{\mbox{\bf H}}_{1,2})_{i,i}&(\breve{\mbox{\bf H}}_{1,3})_{i,i}&\cdots&\cdots&(\breve{\mbox{\bf H}}_{1,p})_{i,i}\\ (\breve{\mbox{\bf H}}_{2,1})_{i,i}&(\breve{\mbox{\bf H}}_{2,2}+\tau^{2}Q_{w,2})_{i,i}&(\breve{\mbox{\bf H}}_{2,3})_{i,i}&\cdots&\cdots&(\breve{\mbox{\bf H}}_{2,p})_{i,i}\\ \cdots&&\cdots&&\cdots&\\ \cdots&&\cdots&&\cdots&\\ (\breve{\mbox{\bf H}}_{p,1})_{i,i}&(\breve{\mbox{\bf H}}_{p,2})_{i,i}&(\breve{\mbox{\bf H}}_{p,3})_{i,i}&\cdots&\cdots&(\breve{\mbox{\bf H}}_{p,p}+\tau^{2}Q_{w,p})_{i,i}\end{bmatrix}\in\mathbb{R}^{p\times p},

and

\widetilde{\mbox{\bf G}}_{i}\,:=\,\begin{bmatrix}(\mathbf{G}_{1,1})_{i,i}&(\mathbf{G}_{1,2})_{i,i}\\ \vdots&\vdots\\ (\mathbf{G}_{p,1})_{i,i}&(\mathbf{G}_{p,2})_{i,i}\end{bmatrix}\in\mathbb{R}^{p\times 2}.

Note that $\widetilde{\mathbf{H}}=\mbox{diag}(\widetilde{\mathbf{H}}_{1},\widetilde{\mathbf{H}}_{2},\ldots,\widetilde{\mathbf{H}}_{n})=E^{T}\mathbf{H}E$ for the permutation matrix $E\in\mathbb{R}^{pn\times pn}$ given by (9). Since $\mathbf{H}$ is PD, so are all the $\widetilde{\mbox{\bf H}}_{i}$ ’s.

As examples, we give the closed form expressions of $\widetilde{\mbox{\bf H}}_{i}$ and $\widetilde{\mbox{\bf G}}_{i}$ for some small $p$ ’s. When $p=1$ , $\widetilde{\mathbf{H}}_{i}=\tau^{2}\big{(}\frac{\tau^{2}}{4}\alpha^{1}_{i}+\beta^{1}_{i}+\zeta^{1}_{i}\big{)}$ and $\widetilde{\mathbf{G}}_{i}=\begin{bmatrix}\frac{\tau^{2}}{2}\alpha^{1}_{i}&\frac{\tau^{3}}{2}\alpha^{1}_{i}+\tau\beta^{1}_{i}\end{bmatrix}$ for each $i=1,\ldots,s$ . When $p=2$ , we have, for each $i=1,\ldots,n$ ,

\widetilde{\mathbf{H}}_{i}=\tau^{2}\begin{bmatrix}\frac{\tau^{2}}{4}\alpha^{1}_{i}+\frac{9\tau^{2}}{4}\alpha^{2}_{i}+\beta^{1}_{i}+\beta^{2}_{i}+\zeta^{1}_{i}&\frac{3\tau^{2}}{4}\alpha^{2}_{i}+\beta^{2}_{i}\\ \frac{3\tau^{2}}{4}\alpha^{2}_{i}+\beta^{2}_{i}&\frac{\tau^{2}}{4}\alpha^{2}_{i}+\beta^{2}_{i}+\zeta^{2}_{i}\end{bmatrix}\in\mathbb{R}^{2\times 2},

and

\widetilde{\mathbf{G}}_{i}=\begin{bmatrix}\tau^{2}\frac{\alpha^{1}_{i}+3\alpha^{2}_{i}}{2}&\frac{\tau^{3}}{2}\alpha^{1}_{i}+3\tau^{3}\alpha^{2}_{i}+\tau(\beta^{1}_{i}+\beta^{2}_{i})\\ \frac{\tau^{2}}{2}\alpha^{2}_{i}&\tau^{3}\alpha^{2}_{i}+\tau\beta^{2}_{i}\end{bmatrix}\in\mathbb{R}^{2\times 2}.

Lemma 5.1.

Let $p=2$ . For any $\tau>0$ , $(\alpha^{1}_{i},\beta^{1}_{i},\zeta^{1}_{i})>0$ and $0\neq(\alpha^{2}_{i},\beta^{2}_{i},\zeta^{2}_{i})\geq 0$ for each $i=1,\ldots,n$ , the matrix $A_{\mbox{c}}(\boldsymbol{\alpha},\boldsymbol{\beta},\boldsymbol{\zeta},\tau)$ is Schur stable, i.e., its spectral radius is strictly less than 1.

Proof.

By the previous argument, it suffices to show that each $\widetilde{A}_{i}$ is Schur stable for $i=1,\ldots,n$ . Fix an arbitrary $i$ . Letting $\widetilde{\mathbf{K}}_{i}=\begin{bmatrix}c_{1}&c_{2}\end{bmatrix}$ , we have

\widetilde{A}_{i}=\begin{bmatrix}1+\frac{\tau^{2}}{2}c_{1}&\tau+\frac{\tau^{2}}{2}c_{2}\\ \tau c_{1}&1+\tau c_{2}\end{bmatrix},

where

c_{1}=-\frac{d_{2}\alpha^{1}_{i}+\alpha^{2}_{i}(2\beta^{2}_{i}+3\zeta^{2}_{i})}{2d^{\prime}},\qquad c_{2}=-\frac{d_{2}(\frac{\tau^{2}}{2}\alpha^{1}_{i}+\beta^{1}_{i})+\frac{3\tau^{2}}{2}\alpha^{2}_{i}\beta^{2}_{i}+\zeta^{2}_{i}(3\tau^{2}\alpha^{2}_{i}+\beta^{2}_{i})}{\tau d^{\prime}},

and $d^{\prime}:=\det(\widetilde{\mathbf{H}}_{i})/\tau^{4}$ , and $d_{s}=\frac{\tau^{2}}{4}\alpha^{s}_{i}+\beta^{s}_{i}+\zeta^{s}_{i}$ for $s=1,2$ . Hence, $d^{\prime}=d_{1}d_{2}+\tau^{2}\alpha^{2}_{i}(\beta^{2}_{i}+\frac{9}{4}\zeta^{2}_{i})+\beta^{2}_{i}\zeta^{2}_{i}$ . Define

\alpha^{\prime}:=d_{2}\alpha^{1}_{i}+\alpha^{2}_{i}(2\beta^{2}_{i}+3\zeta^{2}_{i}),\ \quad\ \beta^{\prime}:=d_{2}\beta^{1}_{i}+\frac{\tau^{2}}{2}\alpha^{2}_{i}\beta^{2}_{i}+\zeta^{2}_{i}\Big{(}\frac{3}{2}\tau^{2}\alpha^{2}_{i}+\beta^{2}_{i}\Big{)},\ \quad\ \gamma^{\prime}:=d_{2}\zeta^{1}_{i}.

Clearly, $\alpha^{\prime},\beta^{\prime},\gamma^{\prime}$ are all positive for any $\tau>0$ , $(\alpha^{1}_{i},\beta^{1}_{i},\zeta^{1}_{i})>0$ and $0\neq(\alpha^{2}_{i},\beta^{2}_{i},\zeta^{2}_{i})\geq 0$ . Moreover, we deduce from a somewhat lengthy but straightforward calculation that $d^{\prime}=\frac{\tau^{2}}{4}\alpha^{\prime}+\beta^{\prime}+\gamma^{\prime}>0$ . Hence, $c_{1}=-\frac{\alpha^{\prime}}{2d^{\prime}}$ , $c_{2}=-\frac{\tau^{2}\alpha^{\prime}/2+\beta^{\prime}}{\tau d^{\prime}}$ , and

\widetilde{A}_{i}=\begin{bmatrix}1-\frac{\alpha^{\prime}\tau^{2}}{4d^{\prime}}&\tau\Big{(}1-\big{(}\frac{\alpha^{\prime}\tau^{2}}{4}+\frac{\beta^{\prime}}{2}\big{)}\frac{1}{d^{\prime}}\Big{)}\\ -\frac{\alpha^{\prime}\tau}{2d^{\prime}}&1-\Big{(}\frac{\alpha^{\prime}\tau^{2}}{2}+\beta^{\prime}\Big{)}\frac{1}{d^{\prime}}\end{bmatrix}\in\mathbb{R}^{2\times 2}.

It follows from [5, Proposition 5.1] that $\widetilde{A}_{i}$ is Schur stable, and so is $A_{\mbox{c}}(\boldsymbol{\alpha},\boldsymbol{\beta},\boldsymbol{\zeta},\tau)$ . ∎

Using a similar technique but more lengthy calculations, it can be shown that when $p=3$ , the matrix $A(\boldsymbol{\alpha},\boldsymbol{\beta},\boldsymbol{\zeta},\tau)$ is Schur stable for $\tau>0$ , $(\alpha^{1}_{i},\beta^{1}_{i},\zeta^{1}_{i})>0$ , $0\neq(\alpha^{2}_{i},\beta^{2}_{i},\zeta^{2}_{i})\geq 0$ and $0\neq(\alpha^{3}_{i},\beta^{3}_{i},\zeta^{3}_{i})\geq 0$ for each $i=1,\ldots,n$ . For $p>4$ , we expect that the same result holds (supported by numerical experience) although its proof becomes much more complicated. Nevertheless, it is observed that in the $p$ -horizon MPC, when the parameters $\alpha^{s}_{i}$ , $\beta^{s}_{i}$ (and possibly including $\zeta^{s}_{i}$ ) with $s\geq 3$ are medium or large, large control inputs are generated, which causes control or speed saturation and may lead to undesired close-loop dynamics. Motivated by this observation, we obtain the following stability result for small $(\alpha^{s}_{i},\beta^{s}_{i})\geq 0$ for $s=3,\ldots,p$ .

Proposition 5.2.

Let $p\geq 3$ . For any $\tau>0$ , $(\alpha^{1}_{i},\beta^{1}_{i},\zeta^{1}_{i})>0$ and $0\neq(\alpha^{2}_{i},\beta^{2}_{i},\zeta^{2}_{i})\geq 0$ for each $i=1,\ldots,n$ , and $\zeta^{s}_{i}>0$ for $s=3,\ldots,p$ and $i=1,\ldots,n$ , there exists a positive constant $\overline{\varepsilon}$ such that for any $\alpha^{s}_{i},\beta^{s}_{i}\in[0,\overline{\varepsilon})$ for $s=3,\ldots,p$ and $i=1,\ldots,n$ , the matrix $A_{\mbox{c}}(\boldsymbol{\alpha},\boldsymbol{\beta},\boldsymbol{\zeta},\tau)$ is Schur stable.

Proof.

Consider $p\geq 3$ . Fix arbitrary $\tau>0$ , $(\alpha^{1}_{i},\beta^{1}_{i},\zeta^{1}_{i})>0$ , and $0\neq(\alpha^{2}_{i},\beta^{2}_{i},\zeta^{2}_{i})\geq 0$ for each $i=1,\ldots,n$ and $\zeta^{s}_{i}>0$ for $s=3,\ldots,p$ and $i=1,\ldots,n$ . Suppose $\alpha^{s}_{i}=\beta^{s}_{i}=0$ for all $s=3,\ldots,p$ and $i=1,\ldots,n$ . Then $Q_{z,s}=Q_{z^{\prime},s}=0$ for all $s\geq 3$ . Hence, ${\mbox{\bf H}}_{i,j}=0$ for all $i\geq 3$ and any $j$ . Thus it is easy to show that for each $i=1,\ldots,n$ ,

\widetilde{\mathbf{H}}_{i}=\begin{bmatrix}\widetilde{\mathbf{H}}^{2}_{i}&&&\\ &\tau^{2}\zeta^{3}_{i}&&\\ &&\ddots&\\ &&&\tau^{2}\zeta^{p}_{i}\end{bmatrix}\in\mathbb{R}^{p\times p},\qquad\widetilde{\mathbf{G}}_{i}=\begin{bmatrix}\widetilde{\mathbf{G}}^{2}_{i}\\ 0\\ \vdots\\ 0\end{bmatrix}\in\mathbb{R}^{p\times 2},

where $\widetilde{\mathbf{H}}^{2}_{i}\in\mathbb{R}^{2\times 2}$ and $\widetilde{\mathbf{G}}^{2}_{i}\in\mathbb{R}^{2\times 2}$ correspond to $p=2$ given before. Hence, $\widetilde{\mathbf{K}}_{i}:=-\mathbf{e}^{T}_{1}\widetilde{\mathbf{H}}^{-1}_{i}\widetilde{\mathbf{G}}_{i}=-\mathbf{e}^{T}_{1}(\widetilde{\mathbf{H}}^{2}_{i})^{-1}\widetilde{\mathbf{G}}^{2}_{i}$ . It follows from Lemma 5.1 that $A_{\mbox{c}}(\boldsymbol{\alpha},\boldsymbol{\beta},\boldsymbol{\zeta},\tau)$ is Schur stable, i.e., its spectral radius is strictly less than 1. Since the spectral radius of $A_{\mbox{c}}(\boldsymbol{\alpha},\boldsymbol{\beta},\boldsymbol{\zeta},\tau)$ is continuous in $\alpha^{s}_{i},\beta^{s}_{i}$ for all $s=3,\ldots,p$ and $i=1,\ldots,n$ , a small perturbation to $\alpha^{s}_{i},\beta^{s}_{i}$ for all $s=3,\ldots,p$ and $i=1,\ldots,n$ still leads to the Schur stable matrix $A_{\mbox{c}}$ . This yields the desired result. ∎

Based on the above results, one may choose $Q_{z,s},Q_{z^{\prime},s},Q_{w,s}$ in the following way. Let $u_{s},v_{s}\in\mathbb{R}^{n}_{+}$ and $w_{s}\in\mathbb{R}^{n}_{++}$ be positive or nonnegative vectors of the same order. Let $\eta>1$ (e.g., $\eta=5$ or higher) be a constant and let $\kappa_{z},\kappa_{z^{\prime}}$ and $\kappa_{w}$ be some positive constants. Then let

Q_{z,s}=\frac{\kappa_{z}}{(s-1)^{\eta}}\mbox{diag}(u_{s}),\quad Q_{z^{\prime},s}=\frac{\kappa_{z^{\prime}}}{(s-1)^{\eta}}\mbox{diag}(v_{s}),\quad Q_{w,s}=\frac{\kappa_{w}}{(s-1)^{\eta}}\mbox{diag}(w_{s}),\qquad s=2,\ldots,p.

6 Numerical Results

6.1 Numerical Experiments and Weight Matrices Design

We conduct numerical tests to evaluate the performance of the proposed fully distributed schemes and the CAV platooning control. In these tests, we consider a platoon of an uncontrolled leading vehicle labeled by the index 0 and ten (i.e., $n=10$ ) CAVs following the leading vehicle. The following physical parameters are used for the CAVs and their constraints throughout this section unless otherwise stated: the desired spacing $\Delta=50m$ , the vehicle length $L=5m$ , the sample time $\tau=1s$ , the reaction time $r=\tau=1s$ , the acceleration and deceleration limits $a_{\max}=1.35m/s^{2}$ and $a_{\min}=-8m/s^{2}$ , and the speed limits $v_{\max}=27.78m/s$ and $v_{\min}=10m/s$ . The initial state of the platoon is $z(0)=z^{\prime}(0)=0$ and $v_{i}(0)=25m/s$ for all $i=0,1,\ldots,n$ . Further, the vehicle communication network is given by the cyclic-like graph, i.e., the bidirectional edges of the graph are $(1,2),(2,3),\ldots,(n-1,n)\in\mathcal{E}$ .

When $n=10$ , a particular choice of these weight matrices is given as follows: for $p=1$ ,

$\displaystyle\boldsymbol{\alpha}^{1}$	$\displaystyle=$	$\displaystyle\big{(}38.85,40.2,41.55,42.90,44.25,45.60,46.95,48.30,49.65,51.00\big{)}:=\widetilde{\boldsymbol{\alpha}},$
$\displaystyle\boldsymbol{\beta}^{1}$	$\displaystyle=$	$\displaystyle\big{(}130.61,136.21,141.82,147.42,153.03,158.64,164.24,169.85,175.46,181.06\big{)}:=\widetilde{\boldsymbol{\beta}},$
$\displaystyle\boldsymbol{\zeta}^{1}$	$\displaystyle=$	$\displaystyle\big{(}62,74,90,92,106,194,298,402,454,480\big{)}:=\widetilde{\boldsymbol{\zeta}}.$

For $p\geq 2$ , we choose $\boldsymbol{\alpha}^{1}=\widetilde{\boldsymbol{\alpha}}-\mathbf{1}$ , $\boldsymbol{\beta}^{1}=\widetilde{\boldsymbol{\beta}}-\mathbf{1}$ , $\boldsymbol{\zeta}^{1}=\widetilde{\boldsymbol{\zeta}}-\mathbf{1}$ , and

\boldsymbol{\alpha}^{s}=\frac{0.0228}{(s-1)^{4}}\times\widetilde{\boldsymbol{\alpha}},\quad\boldsymbol{\beta}^{s}=\frac{0.044}{(s-1)^{4}}\times\widetilde{\boldsymbol{\beta}},\quad\boldsymbol{\zeta}^{s}=\frac{0.0026}{(s-1)^{4}}\times\widetilde{\boldsymbol{\zeta}},\quad s=2,\ldots,p.

The above vectors $\boldsymbol{\alpha}^{s},\boldsymbol{\beta}^{s},\boldsymbol{\zeta}^{s}$ define the weight matrices $Q_{z,s},Q_{z^{\prime},s},Q_{w,s}$ for $s=1,\ldots,5$ , which further yield the closed loop dynamics matrix $A_{\mbox{c}}$ ; see the discussions below (25). It is shown that when these weights are used, the spectral radius of $A_{\mbox{c}}$ is $0.8498$ for $p=1$ , and $0.8376$ for $p=2,\ldots,5$ , respectively.

Discussion on the selection of MPC horizon. We discuss the choice of the MPC prediction horizon $p$ based on numerical tests as follows. Our numerical experience shows that for $p>1$ , the weight matrices $Q_{z,1},Q_{z^{\prime},1}$ and $Q_{w,1}$ play a more important role for the closed loop dynamics. For fixed $Q_{z,1},Q_{z^{\prime},1}$ and $Q_{w,1}$ with the large penalties in $Q_{z,s},Q_{z^{\prime},s}$ and $Q_{w,s}$ for $s>1$ , the closed loop dynamics may be mildly improved but at the expense of undesired large control. Hence, we choose smaller penalties in $Q_{z,s},Q_{z^{\prime},s}$ and $Q_{w,s}$ for $s>1$ , which only lead to slightly better closed loop performance compared with the case of $p=1$ . Further, when a large $p$ is used, the underlying optimization problem has a larger size, resulting in longer computation time and slow convergence of the proposed distributed scheme. Besides, the current MPC model assumes that the future $u_{0}(k+s)=u_{0}(k)$ for all $s=1,\ldots,p-1$ at each $k$ . This assumption is invalid when the true $u_{0}(k+s)$ is substantially different from $u_{0}(k)$ , which implies that the prediction performance is poor for a large $p$ . For these reasons, it is recommended that a smaller $p$ be used, for example, $p\leq 5$ .

The following scenarios are used to evaluate the proposed CAV platooning control.

$\bullet$ Scenario 1: The leading vehicle performs instantaneous deceleration/acceleration and then keeps a constant speed for a while. The goal of this scenario is to test if the platoon can maintain stable spacing and speed when the leading vehicle is subject to acceleration or deceleration disturbances. The motion profile of the leading vehicle is as follows: the leading vehicle decelerates from $k=51s$ to $k=54s$ with the deceleration $-2m/s^{2}$ , and maintains a constant speed till $k=100s$ . After $k=100s$ , it restores to its original speed $25m/s$ with the acceleration $1m/s^{2}$ .

$\bullet$ Scenario 2: The leading vehicle performs periodical acceleration/deceleration. The goal of this scenario is to test whether the proposed control scheme can reduce periodical spacing and speed fluctuation. The motion profile of the leading vehicle in this scenario is as follows: the leading vehicle periodically changes its acceleration and deceleration from $k=51s$ to $k=100s$ with the period $T=4s$ and acceleration/deceleration $\pm 1m/s^{2}$ . Then it maintains its original constant speed $25m/s$ after $k=100s$ .

$\bullet$ Scenario 3: In this scenario, we aim to test the performance of the proposed control scheme in a real traffic environment, particularly when the leading vehicle undergoes traffic oscillations. We use real world trajectory data from an oscillating traffic flow to generate the leading vehicle’s motion profile. Specifically, we consider NGSIM data on eastbound I-80 in San Francisco Bay area in California. We use the data of position and speed of a real vehicle to generate its control input at each second and treat this vehicle as a leading vehicle. Since the maximum of acceleration of this vehicle is close to $2m/s^{2}$ , we choose $a_{\max}=2m/s^{2}$ . All the other parameters or physical limits remain the same. The experiment setup of this scenario is: $z_{i}(0)=0m$ , $v_{i}(0)=25m/s$ for each $i$ , and the time length is $45s$ . To further test the proposed CAV platooning control in a more realistic traffic setting in Scenario 3, random noise is added to each CAV to simulate dynamical disturbances, model mismatch, signal noise, communication delay, and road condition perturbations. In particular, at each $k$ , the random noise with the normal distribution $\mathcal{N}(0.04,0)$ is added to the first CAV, and the noise with the normal distribution $\mathcal{N}(0.02,0)$ is added to each of the rest of the CAVs. Here a larger noise is imposed to the first CAV since there are more noises and disturbances between the leading vehicle and the first CAV.

6.2 Performance of Fully Distributed Schemes and CAV Platooning Control

The generalized Douglas-Rachford splitting method based distributed algorithm (i.e., Algorithm 1) is tested. For each MPC horizon $p$ , the parameters $\alpha$ , $\rho$ , and the error tolerance for the stopping criteria in this algorithm are chosen to achieve desired numerical accuracy and efficiency; see the discussions below (19) for error tolerances and Table 1 for a list of these parameters and error tolerances. In particular, we choose a larger error tolerance for a larger $p$ to meet the desired computation time requirement of one second per vehicle. For comparison, we also test the three operator splitting based distributed scheme and its accelerated version given in Remark 4.2, where we choose $\delta_{i}=\lambda_{\min}(\grave{W}^{i})/2$ , $\gamma=1.9/{\widehat{L}}$ and $\lambda=1.05$ . Here $\widehat{L}$ is the Lipschitz constant defined in Remark 4.2. For the accelerated scheme, we let $\eta=0.2$ and $\gamma_{0}=1.9/(0.8\times\widehat{L})$ .

Table 1: Parameters in Algorithm 1 for different MPC horizon

p

’s

MPC horizon	$p=1$	$p=2$	$p=3$	$p=4$	$p=5$
$\alpha$	0.95	0.95	0.95	0.8	0.8
$\rho$	0.3	0.3	0.3	0.1	0.1
Error tolerance	$10^{-3}$	$2\times 10^{-3}$	$5\times 10^{-3}$	$7\times 10^{-3}$	$1.25\times 10^{-2}$

Initial guess warm-up. For a given $p$ , the augmented locally coupled optimization problem (17) has nearly $3np$ scalar variables and $3np$ scalar constraints when the cyclic-like network topology is considered. These sizes can be even larger for other network topologies satisfying the assumption A1. Hence, when $p$ is large, the underlying optimization problem is of large size, which may affect the numerical performance of the distributed schemes. Several techniques are developed to improve the efficiency of the proposed Douglas-Rachford distributed schemes for real-time computation, particularly for a large $p$ . For illustration, we discuss the initial guess warm-up technique as follows. When implementing the proposed scheme, we often choose a numerical solution obtained from the last step as an initial guess for the current step and run the proposed Douglas-Rachford scheme. Such the choice of an initial guess usually works well when two neighboring control solutions are relatively close. However, it is observed that the convergence of the proposed distributed scheme is sensitive to an initial guess, especially when the CAV platoon is subject to significant traffic oscillations, which results in highly different control solutions between two neighboring instants. In this case, using a neighboring control solution as an initial guess leads to rather slow convergence. To solve this problem, we propose an initial guess warm-up technique, motivated by the fact that control solutions are usually unconstrained for most of $k$ ’s. Specifically, we first compute an unconstrained solution in a fully distributed manner, which can be realized by setting $\mathcal{P}_{i}$ as the Euclidean space in Algorithm 1. This step can be efficiently computed since the proximal operator is formulated by an unconstrained quadratic program and has a closed form solution. In fact, letting $J_{i}(\widehat{\mathbf{u}}_{i})=\frac{1}{2}\widehat{\mathbf{u}}^{T}_{i}\widehat{W}_{i}\widehat{\mathbf{u}}_{i}+c^{T}_{\mathcal{I}_{i}}\widehat{\mathbf{u}}_{i}$ , the closed form solution to the proximal operator is given by $\mbox{Prox}_{\rho J_{i}}(\widehat{\mathbf{u}}^{\prime}_{i})=-(\rho\widehat{W}_{i}+I)^{-1}(\rho c_{\mathcal{I}_{i}}-\widehat{\mathbf{u}}^{\prime}_{i})$ , where $\widehat{W}_{i}$ is PD. We then project this unconstrained solution onto the constrained set in one step. Due to the decoupled structure of the problem (17), this one-step projection can be computed in a fully distributed manner. We thus use this projected solution as an initial guess for the Douglas-Rachford scheme. Numerical experience shows that this new initial guess significantly improves computation time and solution quality when $p$ is large.

Table 2: Scenario 1: computation time and numerical accuracy

MPC horizon	Computation time per CAV (s)		Relative numer. error
MPC horizon	Mean	Variance	Mean	Variance
$p=1$	0.0248	0.0017	$3.4\times 10^{-4}$	$1.9\times 10^{-7}$
$p=2$	0.0603	0.0034	$1.5\times 10^{-3}$	$2.6\times 10^{-6}$
$p=3$	0.1596	0.0764	$3.2\times 10^{-3}$	$1.1\times 10^{-5}$
$p=4$	0.1528	0.1500	$4.0\times 10^{-3}$	$1.7\times 10^{-5}$
$p=5$	0.2365	0.2830	$6.6\times 10^{-3}$	$5.7\times 10^{-5}$

Table 3: Scenario 2: computation time and numerical accuracy

MPC horizon	Computation time per CAV (s)		Relative numer. error
MPC horizon	Mean	Variance	Mean	Variance
$p=1$	0.0464	0.0039	$4.0\times 10^{-4}$	$1.9\times 10^{-7}$
$p=2$	0.1086	0.0153	$1.1\times 10^{-3}$	$1.4\times 10^{-6}$
$p=3$	0.3296	$0.2593$	$3.2\times 10^{-3}$	$1.13\times 10^{-5}$
$p=4$	0.5049	0.6257	$5.9\times 10^{-3}$	$4.6\times 10^{-5}$
$p=5$	0.5784	0.7981	$1.13\times 10^{-2}$	$1.3\times 10^{-5}$

Table 4: Scenario 3: computation time and numerical accuracy

MPC horizon	Computation time per CAV (s)		Relative numer. error
MPC horizon	Mean	Variance	Mean	Variance
$p=1$	0.0825	0.0023	$1.30\times 10^{-3}$	$3.5\times 10^{-6}$
$p=2$	0.2011	0.0051	$7.5\times 10^{-3}$	$1.6\times 10^{-4}$
$p=3$	0.5830	0.3462	$1.20\times 10^{-2}$	$4.2\times 10^{-4}$
$p=4$	0.8904	0.4685	$1.69\times 10^{-2}$	$3.3\times 10^{-4}$
$p=5$	0.9967	0.7467	$3.25\times 10^{-2}$	$1.3\times 10^{-4}$

Table 5: Scenario 3: computation time and numerical accuracy with initial guess warm-up

MPC horizon	Computation time per CAV (s)		Relative numer. error
MPC horizon	Mean	Variance	Mean	Variance
$p=1$	0.0243	0.0023	$5.0\times 10^{-4}$	$7.0\times 10^{-7}$
$p=2$	0.0097	$0.0017$	$2.6\times 10^{-3}$	$1.6\times 10^{-5}$
$p=3$	0.0579	$0.0253$	$2.2\times 10^{-3}$	$1.1\times 10^{-5}$
$p=4$	0.1063	0.1103	$3.7\times 10^{-3}$	$2.4\times 10^{-5}$
$p=5$	0.1258	0.1155	$8.5\times 10^{-3}$	$1.5\times 10^{-5}$

Performance of distributed schemes. Distributed algorithms are implemented on MATLAB and run on a computer of the following processor with 4 cores: Intel(R) Core(TM) i7-8550U CPU @ $1.80GHz$ and RAM: $16.0GB$ . We test the fully distributed Algorithm 1 for Scenarios 1-3. At each $k\in\mathbb{N}$ , we use the optimal solution obtained from the last step as an initial guess unless otherwise stated. To evaluate the numerical accuracy of the proposed distributed scheme, we compute the relative error between the numerical solution from the distributed scheme and that from a high precision centralized scheme when the latter solution, labeled as the true solution, is nonzero. The mean and variance of computation time per vehicle and relative errors for different MPC horizon $p$ ’s in noise-free Scenarios 1-3 are displayed in Table 2-4, respectively. The numerical performance for Scenario 3 under noises is similar to that without noise; it is thus omitted due to space limit.

It is observed from the numerical results that when the MPC horizon $p$ increases, more computation time is needed with mildly deteriorating numerical accuracy. This observation agrees with the discussion on the choice of $p$ given in Section 6.1, which suggests a relatively small $p$ for practical computation. Besides, we have tested the proposed initial guess warm-up technique on Scenario 3 for different $p$ ’s using the same parameters and error tolerances for Algorithm 1; see Table 1. To compute a warm-up initial guess using an iterative distributed scheme, we use the same $\alpha$ and $\rho$ for each $p$ with error tolerance $5\times 10^{-4}$ for $p=1$ and $10^{-3}$ for the other $p$ ’s. A summary of the numerical results is shown in Table 5. Compared with the results given in Table 4 without initial guess warm-up, the averaging computation time is reduced by at least $80\%$ and the relative numerical error is reduced by at least two thirds for $p\geq 2$ when the initial guess warm-up is used. This shows that the initial guess warm-up technique considerably improves the numerical efficiency and accuracy, and it is especially suitable for real-time computation when a large $p$ is used. Hence, we conclude that Algorithm 1, together with the initial guess warm-up technique, is suitable for real-time computation with satisfactory numerical precision.

We have also tested the three-operator splitting based distributed scheme and its accelerated version given in Remark 4.2. These schemes provide satisfactory computation time and numerical accuracy when $p$ is small. For example, when $p=1$ , the mean of computation time per CAV is 0.0553 seconds with the variance 0.0284 for Scenario 1 and 0.219 seconds with the variance 0.138 for Scenario 2, respectively. However, for a slightly large $p$ , e.g., $p\geq 3$ , it takes much longer than 1 second for an individual CAV to complete computation. These schemes are thus not suitable for real-time computation for $p\geq 3$ .

Refer to caption — (a) Time history of spacing changes.

Performance of the CAV platooning control. We discuss the closed-loop performance of the proposed CAV platooning control for the three aforementioned scenarios with different MPC horizon $p$ ’s. In each scenario, we evaluate the performance of the spacing between two neighboring vehicles (i.e., $S_{i-1,i}(k):=x_{i-1}(k)-x_{i}(k)=z_{i}(k)+\Delta$ ), the vehicle speed $v_{i}(k)$ , and the control input $u_{i}(k),\,i=1,\ldots,n$ for $p=1,2,3,4,5$ . Due to the paper length limit, we present the closed-loop performance for $p=1$ and $p=5$ only for each scenario; see Figures 1-3 for (noise free) Scenarios 1-3 respectively, and Figure 4 for Scenario 3 with noises. In fact, it is observed from these figures (and the other tests) that there is little difference in control performance between $p=1$ and a higher $p$ , e.g., $p=5$ . We comment more on the closed-loop performance of each scenario as follows:

(i)

Scenario 1. Figure 1 shows that the spacing between the uncontrolled leading vehicle and the first CAV, i.e., $S_{0,1}$ , has mild deviation from the desired $\Delta$ when the leading vehicle performs instantaneous acceleration or deceleration, while the spacings between the other CAVs remain the desired constant $\Delta$ . Further, it takes about $35s$ for $S_{0,1}$ to converge to the steady state with the maximum spacing deviation $2.66m$ . The similar performance can be observed for the vehicle speed and control input. In particular, it can be seen that all the CAVs show the exactly same speed change and control, implying that the CAV platoon performs a “coordinated” motion with “consensus” under the proposed platooning control.
(ii)

Scenario 2. Figure 2 displays that under the periodic acceleration/deceleration of the leading vehicle, the CAV platoon also demonstrates a “coordinated” motion with “consensus”. For example, only $S_{0,1}$ demonstrates mild fluctuation, whereas the spacings between the other CAVs remain the desired constant, and all the CAVs show the exactly same speed change and control. Moreover, under the proposed platooning control, the oscillations of $S_{0,1}$ are relatively small with the maximal magnitude less than $0.22m$ . Such oscillations quickly decay to zero within $30s$ when the leading vehicle stops its periodical acceleration/deceleration.
(iii)

Scenario 3. In this scenario, the leading vehicle undergoes various traffic oscillations through the time window of $45s$ . In spite of such oscillations, it is seen from Figure 3 that only $S_{0,1}$ demonstrates small spacing variations with the maximum magnitude less than $1m$ , but the spacings between the other CAVs remain almost constant $\Delta$ through the entire time window. This shows that the CAV platoon also demonstrates a “coordinated” motion with “consensus” as in Scenarios 1-2.
(iv)

Scenario 3 subject to noises. Figure 4 shows the control performance of the CAV platoon in Scenario 3 under noises. It can be seen that there are more noticeable spacing deviations from the desired constant $\Delta$ for all CAVs due to the noises. However, the variation of $S_{0,1}$ remains to be within $1m$ , and the maximum deviation of each $S_{i-1,i}$ with $i\geq 2$ is less than $0.5m$ . Further, the profiles of the CAV speed and control still demonstrate a nearly “coordinated” motion in spite of the noises.

In summary, the proposed platooning control effectively mitigates traffic oscillations of the spacing and vehicle speed of the platoon; it actually achieves a (nearly) consensus motion of the entire CAV platoon even under small random noises and perturbations. Compared with other platoon centered approaches, e.g., [5], the proposed control scheme performs better since it uses different weight matrices that lead to decoupled closed loop dynamics; this choice of the weight matrices also facilitates the development of fully distributed computation.

7 Conclusion

The present paper develops fully distributed optimization based MPC schemes for CAV platooning control under the linear vehicle dynamics. Such schemes do not require centralized data processing or computation and are thus applicable to a wide range of vehicle communication networks. New techniques are exploited to develop these schemes, including a new formulation of the MPC model, a decomposition method for a strongly convex quadratic objective function, formulating the underlying optimization problem as locally coupled optimization, and Douglas-Rachford method based distributed schemes. Control design and stability analysis of the closed loop dynamics is carried out for the new formulation of the MPC model. Numerical tests are conducted to illustrate the effectiveness of the proposed fully distributed schemes and CAV platooning control. Our future research will address nonlinear vehicle dynamics and robust issues under uncertainties, e.g., model mismatch, sensing errors, and communication delay.

8 Acknowledgements

This research work is supported by the NSF grants CMMI-1901994 and CMMI-1902006. We thank Benjamin Hyatt of University of Maryland Baltimore County for his contribution to the proof of the closed loop stability when $p=3$ .

References

[1] F. Alizadeh and D. Goldfarb. Second-order cone programming. Mathematical Programming, 95(1):3–51, 2003.
[2] C. Bergenhem, S. Shladover, E. Coelingh, C. Englund, and S. Tsugawa. Overview of platooning systems. In Proceedings of the 19th ITS World Congress, Oct 22-26, Vienna, Austria (2012), 2012.
[3] D. Davis and W. Yin. A three-operator splitting scheme and its optimization applications. Set-valued and Variational Analysis, 25(4):829–858, 2017.
[4] S. Gong and L. Du. Cooperative platoon control for a mixed traffic flow including human drive vehicles and connected and autonomous vehicles. Transportation Research Part B: Methodological, 116:25–61, 2018.
[5] S. Gong, J. Shen, and L. Du. Constrained optimization and distributed computation based car following control of a connected and autonomous vehicle platoon. Transportation Research Part B: Methodological, 94:314–334, 2016.
[6] J. Hu, Y. Xiao, and J. Liu. Distributed algorithms for solving locally coupled optimization problems on agent networks. In 2018 IEEE Conference on Decision and Control, pages 2420–2425. IEEE, 2018.
[7] P. Kavathekar and Y. Chen. Vehicle platooning: A brief survey and categorization. In ASME 2011 International Design Engineering Technical Conferences and Computers and Information in Engineering Conference, pages 829–845. American Society of Mechanical Engineers, 2011.
[8] A. Kesting, M. Treiber, M. Schönhof, and D. Helbing. Adaptive cruise control design for active congestion avoidance. Transportation Research Part C: Emerging Technologies, 16(6):668–683, 2008.
[9] J. Koshal, A. Nedic, and U. V. Shanbhag. Multiuser optimization: distributed algorithms and error analysis. SIAM Journal on Optimization, 21(3):1046–1081, 2011.
[10] S. Li, K. Li, R. Rajamani, and J. Wang. Model predictive multi-objective vehicular adaptive cruise control. IEEE Transactions on Control Systems Technology, 19(3):556–566, 2011.
[11] G. Marsden, M. McDonald, and M. Brackstone. Towards an understanding of adaptive cruise control. Transportation Research Part C: Emerging Technologies, 9(1):33–51, 2001.
[12] M. Mesbahi and M. Egerstedt. Graph Theoretic Methods in Multiagent Networks. Princeton University Press, 2010.
[13] R. T. Rockafellar. Convex Analysis. Princeton University Press, 1970.
[14] S. E. Shladover, C. Nowakowski, X.-Y. Lu, and R. Ferlis. Cooperative adaptive cruise control: definitions and operating concepts. Transportation Research Record: Journal of the Transportation Research Board, (2489):145–152, 2015.
[15] S. E. Shladover, D. Su, and X.-Y. Lu. Impacts of cooperative adaptive cruise control on freeway traffic flow. Transportation Research Record: Journal of the Transportation Research Board, 2324:63–70, 2012.
[16] J. F. Sturm. Using SeDuMi 1.02, a MATLAB toolbox for optimization over symmetric cones. Optimization Methods and Software, 11(1-4):625–653, 1999.
[17] B. Van Arem, C. J. Van Driel, and R. Visser. The impact of cooperative adaptive cruise control on traffic-flow characteristics. IEEE Transactions on Intelligent Transportation Systems, 7(4):429–436, 2006.
[18] J. Vander Werf, S. E. Shladover, M. A. Miller, and N. Kourjanskaia. Effects of adaptive cruise control systems on highway traffic flow capacity. Transportation Research Record, 1800(1):78–84, 2002.
[19] J. Wang, S. Gong, S. Peeta, and L. Lu. A real-time deployable model predictive control-based cooperative platooning approach for connected and autonomous vehicles. Transportation Research Part B: Methodological, 128:271–301, 2019.
[20] M. Wang, W. Daamen, S. P. Hoogendoorn, and B. van Arem. Rolling horizon control framework for driver assistance systems. Part II: Cooperative sensing and cooperative control. Transportation Research Part C: Emerging Technologies, 40:290–311, 2014.
[21] M. Wang, W. Daamen, S. P. Hoogendoorn, and B. van Arem. Cooperative car-following control: Distributed algorithm and impact on moving jam features. IEEE Transactions on Intelligent Transportation Systems, 17(5):1459–1471, 2016.
[22] S. Zhao and K. Zhang. A distributionally robust stochastic optimization-based model predictive control with distributionally robust chance constraints for cooperative adaptive cruise control under uncertain traffic conditions. Transportation Research Part B: Methodological, 138:144–178, 2020.
[23] Y. Zheng, S. E. Li, K. Li, F. Borrelli, and J. K. Hedrick. Distributed model predictive control for heterogeneous vehicle platoons under unidirectional topologies. IEEE Transactions on Control Systems Technology, 25(3):899–910, 2016.
[24] Y. Zhou, S. Ahn, M. Chitturi, and D. A. Noyce. Rolling horizon stochastic optimal control strategy for ACC and CACC under uncertainty. Transportation Research Part C: Emerging Technologies, 83:61–76, 2017.
[25] Y. Zhou, M. Wang, and S. Ahn. Distributed model predictive control approach for cooperative car-following with guaranteed local and string stability. Transportation Research part B: Methodological, 128:69–86, 2019.