Distributed Model Predictive Control Under Inexact Primal-Dual Gradient Optimization Based on Contraction Analysis

TL;DR

This paper proposes a distributed model predictive control approach for linear systems with coupled constraints, utilizing primal-dual gradient optimization and contraction analysis to ensure convergence, stability, and reduced computation.

Contribution

It introduces a novel DMPC method based on inexact primal-dual gradient optimization with contraction analysis for convergence and stability guarantees.

Findings

Ensures recursive feasibility and stability under inexact solutions.

Uses contraction theory for convergence analysis of primal-dual methods.

Validates the approach with numerical simulation results.

Abstract

This paper develops a distributed model predictive control (DMPC) strategy for a class of discrete-time linear systems with consideration of globally coupled constraints. The DMPC under study is based on the dual problem concerning all subsystems, which is solved by means of the primal-dual gradient optimization in a distributed manner using Laplacian consensus. To reduce the computational burden, the constraint tightening method is utilized to provide a capability of premature termination with guaranteeing the convergence of the DMPC optimization. The contraction theory is first adopted in the convergence analysis of the primal-dual gradient optimization under discrete-time updating dynamics towards a nonlinear objective function. Under some reasonable assumptions, the recursive feasibility and stability of the closed-loop system can be established under the inexact solution. A numerical simulation is given to verify the performance of the proposed strategy.

Tables

Algorithm 1 The Distributed Primal-Dual Gradient Algorithm

Input: $x_i\left(i=1,\cdots ,l\right)$

Output: ${\widetilde{𝒖}}_i\left(i=1,\cdots ,l\right)$

Required: $ϵ$

Initialization: ${\widetilde{𝒖}}_i\left(0\right)=0$, ${𝝀}_i^0=0$ and proper ${𝜸}_i^0$ for $i=1,\cdots ,l$.

Obtain stopping criterion $\overline{k}$ based on $\left(\text{<a><span>54</span></a>}\right)$;

1: while $k\le \overline{k}$

2:  for $i=1,\cdots ,l$ (in parallel) do

3:   Exchange ${𝝀}_i^k$ with its neighbor $j\in {𝒩}_i$;

4:   Obtain ${𝜸}_i^{k+1}$ from $\left(\text{<a><span>32b</span></a>}\right)$, ${𝝀}_i^{k+1}$ according to $\left(\text{<a><span>32a</span></a>}\right)$,

${\widetilde{𝒖}}_i$ following $\left(\text{<a><span>31</span></a>}\right)$, respectively;

5:  end for

6:  Move $k$ to $k+1$.

7: end while

Algorithm 2 The DMPC Algorithm
1: At time instant $t$, each subsystem $i$ samples its state $x_i\left(t\right)$;
2: Each subsystem $i$ obtains $u_i\left(t\right)$ by following Algorithm 1
with $x_i\left(t\right)$;
3: Each subsystem $i$ apples its control input;
4: Let $t=t+1$ and go to step 1 until next sampling instant.

Equations

$$\xi =f\left(\xi ,k\right),$$

$$d\left(x\left(k\right),x^{\ast }\left(k\right)\right)\le C{\tau }^{k}d\left(x\left(0\right),x^{\ast }\left(0\right)\right),\forall k>0,$$

$$\frac{\partial f^{\mathrm{T}}}{\partial x}M\frac{\partial f}{\partial x}-M\le \left({\tau }^2-1\right)M$$

$$x_i\left(t+1\right)=A{x}_i\left(t\right)+B{u}_i\left(t\right),i=1,2,\cdots ,l,$$

$$\sum _{i=1}^l{\Phi }_i^x{x}_i\left(t\right)+{\Phi }_i^u{u}_i\left(t\right)\le 1_p,\forall t\ge 0,$$

$$J_i\left(x_i,u_i\right)=\sum _{s=0}^{N}F\left(x_i\left(s|t\right),u_i\left(s|t\right)\right)+V_f\left(x_i\left(N|t\right)\right)$$

$$J_i\left(x_i,u_i\right)$$

$$+\frac{m_J}{2}{\left({\nabla }_{u_i}{J}_i\left(x_i,u_i^{\prime }\right)\right)}^{\mathrm{T}}\left(u_i^{\prime }-u_i\right),$$

$$J_i\left(x_i,u_i\right)$$

$$+{\left({\nabla }_{u_i}{J}_i\left(x_i,u_i^{\prime }\right)\right)}^{\mathrm{T}}\left(u_i-u_i^{\prime }\right).$$

$$m_J\mathbf{I}\le {\nabla }_{u_i}^2{J}_i\left(x_i,u_i\right)\le L_J\mathbf{I}$$

$$\nabla J_i\left(x_i,u_i\right)-\nabla J_i\left(x_i,u_i^{\ast }\right)=H_i\left(u_i\right)\left(u_i-u_i^{\ast }\right).$$

$${\tilde{\mathcal{U}}}_i\left(x_i\left(t\right)\right):=$$

$$+B{\tilde{\mathbf{u}}}_i\left(s|t\right),{\tilde{\mathbf{x}}}_i\left(0|t\right)=x_i\left(t\right),{\tilde{\mathbf{x}}}_i\left(s|t\right)\in {\mathcal{X}}_i,$$

$${\tilde{\mathbf{u}}}_i\left(s|t\right)\in {\mathcal{U}}_i,{\tilde{\mathbf{x}}}_i\left(N|t\right)\in {\Omega }_i^f,s\in {\mathbb{Z}}_{\left[0,N-1\right]}\},$$

$$V_f\left(A_K{x}_i\left(t\right)\right)-V_f\left(x_i\left(t\right)\right)\le -F\left(x_i\left(t\right),K{x}_i\left(t\right)\right)$$

$$\sum _{i=1}^l\left({\Phi }_i^x+{\Phi }_i^u{K}_i\right)x_i\left(s|t\right)\le 1_p,\forall s\in {\mathbb{Z}}_{>0}.$$

$$u_i^{\ast }\left(s|t\right)=\arg \underset{{\tilde{\mathbf{u}}}_i\left(s|t\right)}{\min }\sum _{i=1}^l{J}_i\left({\tilde{\mathbf{x}}}_i\left(s|t\right),{\tilde{\mathbf{u}}}_i\left(s|t\right)\right)$$

$$s.t.$$

$$\sum _{i=1}^l{\Phi }_i^x{\tilde{\mathbf{x}}}_i\left(s|t\right)+{\Phi }_i^u{\tilde{\mathbf{u}}}_i\left(s|t\right)\le 1_p,$$

$$\sum _{i=1}^l{\Phi }_i^x{\tilde{\mathbf{x}}}_i\left(s|t\right)+{\Phi }_i^u{\tilde{\mathbf{u}}}_i\left(s|t\right)\le \left(1-ϵl\left(s+1\right)\right)1_p,$$

$$\sum _{i=1}^l{f}_i\left(x_i\left(t\right),{\tilde{\mathbf{u}}}_i\left(s|t\right)\right)\le b\left(ϵ\right),$$

$$f_i\left(x_i\left(s|t\right),{\tilde{\mathbf{u}}}_i\left(s|t\right)\right):=F_i{x}_i\left(s|t\right)+G_i{\tilde{\mathbf{u}}}_i\left(s|t\right)$$

$$b\left(ϵ\right):={\left[\left(1-ϵl\right)1_p^{\mathrm{T}},\left(1-2ϵl\right)1_p^{\mathrm{T}},\cdots ,\left(1-ϵNl\right)1_p^{\mathrm{T}}\right]}^{\mathrm{T}}.$$

$$\underline{\zeta }\mathbf{I}\le G_i^{\mathrm{T}}{G}_i\le \overline{\zeta }\mathbf{I}$$

$$\sum _{i=1}^l\left({\Phi }_i^x+{\Phi }_i^u{K}_i\right)x_i\left(s|t\right)\le \left(1-ϵNl\right)1_p,$$

$${\tilde{\mathbf{u}}}_i^{\ast }\left(s|t\right)=\arg \underset{{\tilde{\mathbf{u}}}_i\left(t\right)}{\min }\sum _{i=1}^l{J}_i\left(x_i\left(t\right),{\tilde{\mathbf{u}}}_i\left(s|t\right)\right)$$

$$s.t.$$

$$\sum _{i=1}^l{f}_i\left(x_i\left(t\right),{\tilde{\mathbf{u}}}_i\left(s|t\right)\right)\le b\left(ϵ\right),$$

$$\mathcal{L}\left(x_i,{\tilde{\mathbf{u}}}_i,\mathbf{\lambda }\right)=\sum _{i=1}^l{J}_i\left(x_i,{\tilde{\mathbf{u}}}_i\right)+{\mathbf{\lambda }}^{\mathrm{T}}\left(\sum _{i=1}^l{f}_i\left(x_i,{\tilde{\mathbf{u}}}_i\right)-b\left(ϵ\right)\right)$$

$$\underset{\mathbf{\lambda }\ge 0}{\max }\underset{{\tilde{\mathbf{u}}}_i\in {\tilde{\mathcal{U}}}_i\left(x_i\right)}{\min }\mathcal{L}\left(x_i,{\tilde{\mathbf{u}}}_i,\mathbf{\lambda }\right)$$

$$\underset{\mathbf{\lambda }\ge 0}{\min }\underset{{\tilde{\mathbf{u}}}_i\in {\tilde{\mathcal{U}}}_i\left(x_i\right)}{\max }-\mathcal{L}\left(x_i,{\tilde{\mathbf{u}}}_i,\mathbf{\lambda }\right).$$

$${\nabla }_{{\tilde{\mathbf{u}}}_i^{\ast }}\mathcal{L}\left({\tilde{\mathbf{u}}}_i^{\ast },{\mathbf{\lambda }}^{\ast }\right)$$

$${\nabla }_{{\mathbf{\lambda }}^{\ast }}\mathcal{L}\left({\tilde{\mathbf{u}}}_i^{\ast },{\mathbf{\lambda }}^{\ast }\right)$$

$$\underset{\mathbf{\lambda }\ge 0}{\min }{\Psi }_i\left(\mathbf{\lambda }\right)$$

Full text

I Introduction

Model predictive control (MPC) is of tremendous interest in recent years for its broad applications ranging from industrial control systems [1], aerospace [2], sensor network control [3], etc. Owing to the increasing capability of computing power and the accelerated optimization algorithms, MPC has been utilized from the traditional process control to fields like complex dynamical systems [4, 5], networked control systems [6], multi-agent systems [7] and so on. It is worth to point out that for reducing the computational burden, the early termination is adopted in the optimization process to fulfill the real-time requirements, leading the solution to be inexact [8, 9].

For large-scale systems, however, distributed MPC (DMPC) should be investigated to further reduce the computation resources compared with the centralized MPC. The most commonly used DMPC formulation is to describe the DMPC as a convex optimization problem that should be solved in a distributed fashion [10, 11, 12], which can be referred to distributed consensus-based optimization [13, 14]. For given multi-agent systems, the goal of all the subsystems is to minimize a global objective function collaboratively without sharing the private objective function throughout the optimization process. Generally speaking, the distributed algorithms for convex optimization fall into three categories known as primal methods, dual methods and primal-dual methods. For primal methods, the convex optimization problem is solved in the primal domain with guaranteeing the consensus by introducing the penalty functions reflecting the disagreement [15]. For dual methods, each subsystem solves the dual problem in a distributed manner to seek consensus [16]. For the primal-dual methods, to obtain the saddle-point, the primal and dual problems are updated simultaneously associated with Lagrangian multipliers [17]. It is worth to note that each of the methods has its advantages determining from the global objective function and the constraints formulations in the optimization problem.

The coupled systems concerning DMPC scheme have received increasing attention recently as dynamical couplings are ubiquitous in practical applications. In [18], a distributed receding horizon control strategy was studied for a class of dynamically coupled nonlinear systems with consideration of decoupled constraints. In [19], the coupled probabilistic constraints were investigated in the context of distributed stochastic MPC. The compromise for satisfying the coupled constraints in a distributed way is to allow only one subsystem to optimize at each time instant, which is a widely utilized method. However, the distributed optimization problems have been rarely studied under the globally coupled constraints in the framework of DMPC. The limitation of the globally coupled constraints render the existing results on coupled systems cannot be generalized directly. The authors in [20] presented a DMPC scheme for a group of discrete-time systems taking the local and global constraints into account. A dual problem was established to solve the DMPC optimization problem based on the Alternating Direction Multiplier Method (ADMM). It, however, is worth mentioning that the objective function considered in [20] was with a quadratic form, which implies it cannot be directly extended to the scenario with a general nonlinear objective function.

For nonlinear systems analysis, a well-known method is the contraction theory. By virtue of the fluid mechanics and differential geometry, the contraction theory has been first introduced in [21]. The traditional approach established by using Riemannian manifolds has been extended to many applications such as distributed nonlinear systems [22], stochastic incremental systems [23], etc. Furthermore, some recent results were developed inspired by Finsler manifolds [24, 25]. Nevertheless, it is worth noting that the existing results were obtained in the context of continuous-time dynamical systems. To the best of our knowledge, the contraction theory has been rarely investigated for discrete-time dynamical systems. Moreover, only a few results have been addressed using the contraction theory to analyze the convergence of the optimization algorithm [26].

This paper formulates the DMPC optimization problem as a distributed consensus optimization problem (DCOP). To provide the capability of early termination with guaranteeing the convergence, a tightening constraint is constructed in the optimization. In addition, the primal-dual gradient optimization is adopted to solve the DCOP. The contraction theory based on Riemannian manifolds is utilized to analyze the convergence of the primal-dual gradient dynamics (PDGD). The contributions of this paper are mainly in three-fold as follows.

1. (1)

This paper investigates the DMPC subject to globally coupled constraints which can be formulated as a distributed consensus-based optimization problem. Furthermore, the inexact solver is taken into account to reduce unnecessary computations. The constraint tightening method is performed to allow premature termination with guaranteeing the convergence of the optimization process. 2. (2)

The objective function in the DMPC optimization problem is considered as a nonlinear function in this paper. The dual problem is utilized to solve the DMPC optimization problem. Moreover, the local copies of the Lagrangian multiplier in the dual problem are introduced to achieve fully distributed. Thereafter, the primal-dual gradient dynamics are established to solve the consensus-based optimization problem. Owing to the tightening constraints, the local copies without being required to reach consensus but need to fulfill some specified bounds. 3. (3)

Inspired by the contraction theory based on Riemannian manifolds, some sufficient conditions are given to guarantee the convergence of the PDGD. It is worth noting that the contraction theory is first used to analyze the optimization convergence in the context of discrete-time updating dynamics. In addition, the recursive feasibility and stability of the DMPC algorithm are rigorously analyzed.

The remainder of this paper is organized as follows. Section II presents some necessary preliminaries adopted throughout this paper. In Section III, the considered optimization problem is formulated. The theoretical results are demonstrated in Section IV. Thereafter, Section V describes the proposed algorithm. We analyze the recursive feasibility and stability in Section VI. A numerical example is given in Section VII to verify the effectiveness of the proposed algorithm. Section VIII summarizes the paper.

The notations used in this paper are stated in the following. Denote real and natural number set as $\mathbb{R}$ and $\mathbb{N}$, respectively. The superscripts $\mathrm{T}$ and $+$ stand for transposition of a given matrix and the successor states. The subscripts of $\mathbb{N}_{\geq a}$ and $\mathbb{N}_{\left[b,\,c\right]}$ represent the integers in the intervals $\left[a,\,\infty\right)$ and $\left[b,c\right]$, respectively. Given two matrices $X$ and $Y$, the Kronecker product is denoted by $X\otimes Y$. Define the $P$-weighted norm of a given vector $x$ as $\left\|x\right\|_{P}=\sqrt{x^{\mathrm{T}}Px}$ with respect to the positive definite matrix $P$.

II Preliminaries

II-A Graph Theory

In this paper, the multi-agent system containing $l$ subsystems communicate with each other according to a weighted undirected graph denoted as $\mathcal{G}=\left(\mathcal{V},\mathcal{E},\mathcal{A}\right)$, where the set of agents is represented as $\mathcal{V}=\left\{1,2,\cdots,l\right\}$, the set $\mathcal{E}=\mathcal{V}\times\mathcal{V}$ collects the undirected edges indicating the interconnections among subsystems, and the adjacent matrix $\mathcal{A}=[a_{ij}]_{l\times l}$ is symmetric with $a_{ii}=0$ implying that there is no self-edge in the graph, $a_{ij}=a_{ji}>0$ denoting the undirected weight if $\left\{i,j\right\}\in\mathcal{E}$ and $a_{ij}=0$ otherwise. Collect the neighbors of vertex $i$ in a set $\mathcal{N}_{i}=\left\{j\in\mathcal{E}:a_{ij}>0\right\}$. The degree matrix is $\mathcal{D}=\text{diag}\{d_{i}\}\in\mathbb{R}^{l\times l}$, where $d_{i}=\sum_{j\in\mathcal{N}_{i}}a_{ij}$. The Laplacian matrix is defined as $L=\mathcal{D}-\mathcal{A}$, and its eigenvalue decomposition is $L=Q\text{diag}\left(\nu_{1},\nu_{2},\cdots,\nu_{M}\right)Q^{\mathrm{T}}$, where $Q$ is an orthogonal matrix and $0=\nu_{1}<\nu_{2}\leq\cdots\leq\nu_{l}$. We can obtain $L=\sqrt{L}\sqrt{L}$ by introducing $\sqrt{L}=Q\text{diag}\left(\sqrt{\nu_{1}},\sqrt{\nu_{2}},\cdots,\sqrt{\nu_{M}}\right)Q^{\mathrm{T}}$.

II-B Contraction Theory

Some crucial definitions on Riemannian geometry are recapped for convergence analysis of the proposed algorithm by using contraction theory. For more details, please refer to [21, 27] and references therein. The Riemannian metric of two vectors $\delta_{1},\delta_{2}$ on the tangent space of a given state manifold is a smoothly varying inner product $\left\langle\delta_{1},\delta_{2}\right\rangle_{x}=\delta_{1}^{\mathrm{T}}M\left(x\right)\delta_{2}$ with respect to a positive matrix function $M\left(x\right)$. Thus, we have $\|\delta\|_{x}=\sqrt{\left\langle\delta_{1},\delta_{2}\right\rangle_{x}}$. In this paper, we assume the matrix $M\left(x\right)=M$ to be constant. Given a pair of points $x\in\delta_{1}$ and $y\in\delta_{2}$, let the set of smooth curves connecting $x$ and $y$ be $\Gamma\left(x,y\right)$. There exists a piecewise smooth mapping $\gamma\in\Gamma\left(x,y\right)$ satisfying $\gamma\left(0\right)=x$ and $\gamma\left(1\right)=y$. Define the Riemannian length as $L\left(\gamma\right):=\int_{0}^{1}\|\gamma_{s}\|_{\gamma}ds$, the Riemannian energy $E\left(\gamma\right):=\int_{0}^{1}\|\gamma_{s}\|_{\gamma}^{2}ds$, where $\gamma_{s}:=\partial\gamma_{s}/\partial s$. Denote the Riemannian distance as $d\left(x,y\right)=\inf_{\gamma\in\Gamma\left(x,y\right)}L\left(\gamma\right)$. In this paper, we define $E\left(x,y\right):=d\left(x,y\right)^{2}$.

Consider an autonomous nonlinear discrete-time system in the following form

[TABLE]

where $f$ is a smooth and differentiable function with respect to $\xi$. The differential dynamics can be denoted as $\delta\xi\left(k+1\right)=\frac{\partial f}{\partial x}\delta\xi\left(k\right)$. The optimal state trajectory is defined as a forward-complete solution of $\left(\ref{autonomousSys}\right)$. The Euler discretization of an exponentially controllable system possesses geometric convergence speed [28]. Thus, the optimal state trajectory $x^{*}$ is said to be global exponentially controllable if

[TABLE]

where the positive scalar $C$ and the convergence rate $\tau\in\left(0,1\right)$ are independent of the initial states.

In accordance with [21], the contraction region for the discrete-time system is similarly defined as follows.

Definition 1: Given a discrete-time system $x\left(k+1\right)=f\left(x\left(k\right),k\right)$, a region of state space is called a contraction region, if there exists a uniformly positive definite constant metric $M$, such that

[TABLE]

with the convergence rate $\tau\in\left(0,1\right)$.

III Problem Formulation

Consider the multi-agent system containing $l$ subsystems under a weighted undirected graph. Each subsystem can be formulated as the following linear discrete-time dynamics:

[TABLE]

with $x_{i}\left(t\right)\in\mathcal{X}_{i}$ and $u_{i}\left(t\right)\in\mathcal{U}_{i}$, where $x_{i}$ and $u_{i}$ are the state and input vectors of subsystem $i$, $\mathcal{X}_{i}\subset\mathbb{R}^{n}$ and $\mathcal{U}_{i}\subset\mathbb{R}^{m}$ are compact sets of the local constraints on state and input containing the origin as their inner point, respectively. Moreover, the following globally coupled constraint is taken into account:

[TABLE]

where $\Phi_{i}^{x}\in\mathbb{R}^{p\times n}$, $\Phi_{i}^{u}\in\mathbb{R}^{p\times m}$ are the matrices defining the coupled constraint, and $\boldsymbol{1}_{p}$ represents the $p$-vector with all ones.

For subsystem $i$ at time instant $t$, define the objective function $J_{i}\left(x_{i},u_{i}\right):\mathbb{R}^{n}\times\mathbb{R}^{m}\rightarrow\mathbb{R}$ as

[TABLE]

where $N$ is the prediction horizon. The following assumptions are crucial for theoretical analysis in this paper.

Assumption 1: [14] The local objective function $J_{i}\left(x_{i},u_{i}\right)$ is twice-continuously differentiable and $m_{J}$-strongly convex with Lipschitz gradient $L_{J}$ with respect to $u_{i}$, ie., for all $u_{i},u_{i}^{\prime}\in\mathcal{U}_{i}$,

[TABLE]

Remark 1: Assumption 1 is standard in convergence analysis of convex optimization. We can rewrite Eq. $\left(\ref{assumption1}\right)$ into another form as

[TABLE]

with $0<m_{J}\leq L_{J}<\infty$.

Assumption 2: [29] Under Assumption 1, for any $\boldsymbol{u}_{i}$ and the optimal $u^{*}$, there exists an invertible symmetric matrix $H_{i}\left(u_{i}\right)$ satisfying $m_{J}\boldsymbol{I}\leq H_{i}\left(u_{i}\right)\leq L_{J}\boldsymbol{I}$ such that

[TABLE]

Remark 2: Notice that $H_{i}\left(u_{i}\right)$ is a time-varying matrix with respect to $u_{i}$. For convenience, we use $H_{i}$ for short in this paper.

Inspired by [30, 31], we define

[TABLE]

where $\tilde{\boldsymbol{x}}_{i}\left(s|t\right)$ and $\tilde{\boldsymbol{u}}_{i}\left(s|t\right)$ depict the predicted state and input sequences for $s\in\mathbb{N}_{\left[0,N-1\right]}$, respectively, $N\in\mathbb{Z}_{>0}$ represents the prediction horizon and $\Omega_{i}^{f}:=\left\{x_{i}|V_{f}\left(x_{i}\right)\leq\eta_{f}\right\}$ denotes the terminal set. In addition, the following assumption is given for stability analysis.

Assumption 3: [32] Given a positive scalar $\eta>\eta_{f}$, there exists a local state-feedback control gain $K$ such that $u_{i}=Kx_{i}\in\mathcal{U}_{i}$. Moreover, it holds that

[TABLE]

with $A_{K}=A+BK$ for all $x_{i}\in\Omega_{i}:=\left\{x_{i}|V_{f}\left(x_{i}\right)\leq\eta\right\}$.

It is worth to clarify that the terminal set used in $\left(\ref{U_x}\right)$ is smaller than which is defined in Assumption 3. In addition, the following condition should be satisfied for all $x_{i}\left(t\right)\in\Omega_{i}^{f}$ to fulfill the globally coupled constraint in $\left(\ref{coupledConstraint}\right)$.

[TABLE]

In what follows, we formulate the standard MPC optimization problem as

[TABLE]

with $s\in\mathbb{Z}_{\left[0,N-1\right]}$, where $\tilde{\mathcal{U}}_{i}\left(x_{i}\left(t\right)\right)$ is defined in $\left(\ref{U_x}\right)$ and $\left(\ref{standardMPC:coupled}\right)$ implying the satisfaction of the globally coupled constraint in $\left(\ref{coupledConstraint}\right)$.

To accelerate the optimization, one efficient way is to allow premature termination of the optimization problem in $\left(\ref{standardMPC}\right)$ with guaranteeing the convergence, which can avoid unnecessary computations. It is worth mentioning that the early termination may result in errors even infeasibility of the optimization problem. Taking the inexactness of optimization solver into consideration, the $\epsilon$-strict feasibility is investigated by introducing the tightening constraint in this paper.

Definition 2: [9] Given a polytopic constraint as $f\left(x\right)\leq b$, the $\epsilon$-strictly feasible solution is the vector $x$ which satisfies $f\left(x\right)\leq b+\epsilon\boldsymbol{1}$ with $\boldsymbol{1}$ in proper dimension.

Inspired by the tightening constraint method, we can transform the global coupled constraint in $\left(\ref{standardMPC:coupled}\right)$ as follows

[TABLE]

with $s\in\mathbb{Z}_{\left[0,N-1\right]}$ where $\epsilon\in\mathbb{R}_{\left(0,\frac{1}{Nl}\right)}$ is a user-defined parameter reflecting the tolerance to the violation of the coupled constraint. For convenience, we rewrite $\left(\ref{tighteningConstraint}\right)$ as a standard form in optimization.

[TABLE]

with $s\in\mathbb{Z}_{\left[0,N-1\right]}$, where

[TABLE]

with $F_{i}\in\mathbb{R}^{Np\times n}$ and $G_{i}\in\mathbb{R}^{Np\times Nm}$ being appropriate matrices obtained from $\Phi_{i}^{x}$ and $\Phi_{i}^{u}$ by expressing $\tilde{\boldsymbol{x}}_{i}\left(s|t\right)$ in terms of $x_{i}\left(t\right)$ and $\tilde{\boldsymbol{u}}_{i}\left(s|t\right)$, and

[TABLE]

The following assumption is established based on the linear independence property of the coupled constraint.

Assumption 4: There exist two scalars $\underline{\zeta}>0$ and $\overline{\zeta}>0$ such that the full row rank matrix $F_{i}$ satisfying

[TABLE]

where $\boldsymbol{I}$ is with appropriate dimension.

Moreover, the coupled constraint in the terminal set can be rewritten as the following form, correspondingly.

[TABLE]

with $s\in\mathbb{Z}_{>0}$ for $x_{i}\in\Omega_{i}^{f}$. By utilizing the constraint in $\left(\ref{tighteningConstraintOpt}\right)$, the DMPC optimization problem in $\left(\ref{standardMPC}\right)$ can be rewritten as follows

[TABLE]

with $s\in\mathbb{Z}_{\left[0,N-1\right]}$, where $\tilde{\mathcal{U}}_{i}\left(x_{i}\right)$ is defined in $\left(\ref{U_x}\right)$.

IV Main Results

IV-A The Dual Form

To solve the optimization problem in $\left(\ref{erMPC}\right)$, the dual form is considered in the following. The Lagrangian of optimization problem $\left(\ref{erMPC}\right)$ is given as

[TABLE]

for $\tilde{\boldsymbol{u}}_{i}\in\tilde{\mathcal{U}}_{i}\left(x_{i}\right)$, where $\boldsymbol{\lambda}\in\mathbb{R}^{Np}$ is the Lagrangian multiplier. Thus, the dual problem can be described as

[TABLE]

which is equivalent to

[TABLE]

The Karush-Kuhn-Tucker (KKT) conditions for depicting the optimal pair $\left(\tilde{\boldsymbol{u}}_{i}^{*},\boldsymbol{\lambda}^{*}\right)$ are

[TABLE]

Remark 3: For each subsystem, we rewrite the dual problem in $\left(\ref{dualProblem2}\right)$ in the following form

[TABLE]

where

[TABLE]

By Dadnskin’s theorem [33], the dual gradient is obtained as $\nabla\Psi_{i}\left(\boldsymbol{\lambda}\right)=-\left(f_{i}\left(x_{i},\tilde{\boldsymbol{u}}_{i}\left(\boldsymbol{\lambda}\right)\right)-\frac{b\left(\epsilon\right)}{l}\right)$.

IV-B Distributed Optimization with Laplacian Consensus

It is worth to point out that the Lagrangian multiplier $\boldsymbol{\lambda}$ in $\left(\ref{dualProblem3}\right)$ is a global variable such that the optimization problem cannot be solved in a distributed way. Resorting to the Laplacian consensus, the optimization problem in $\left(\ref{dualProblem3}\right)$ can be transformed into

[TABLE]

where $\sqrt{\boldsymbol{L}}:=\sqrt{L}\otimes\boldsymbol{I}_{Np}$ with $L$ being the Laplacian matrix corresponding to the topology of the communication graph, $\boldsymbol{\lambda}_{i}$ is the local copy of $\boldsymbol{\lambda}$ for the subsystem $i$, and $\boldsymbol{\Lambda}:=\left[\boldsymbol{\lambda}_{1}^{\mathrm{T}},\boldsymbol{\lambda}_{2}^{\mathrm{T}},\cdots,\boldsymbol{\lambda}_{l}^{\mathrm{T}}\right]^{\mathrm{T}}$ is the vector which stacks the local Lagrangian multipliers $\boldsymbol{\lambda}_{i}$. For the individual subsystem, we can obtain the Lagrangian of optimization problem in $\left(\ref{consensusOpt}\right)$ as follows

[TABLE]

where $\boldsymbol{\mu}_{i}\in\mathbb{R}^{Np}$ is the Lagrangian multiplier of $i$th subsystem. The KKT conditions for the optimal pair $\left(\boldsymbol{\lambda}_{i}^{*},\boldsymbol{\mu}_{i}^{*}\right)$ are described as

[TABLE]

IV-C Primal-Dual Gradient Dynamics

The primal-dual gradient method is adopted to solve the optimization problems in $\left(\ref{dualProblem3}\right)$ and $\left(\ref{consensusOpt}\right)$. The primal-dual gradient dynamics are given as

[TABLE]

where $\alpha>0$ and $\beta>0$ are the step-sizes, $\nabla\Psi_{i}\left(\boldsymbol{\lambda}_{i}^{k}\right)=-\left(f_{i}\left(x_{i},\tilde{\boldsymbol{u}}_{i}\left(\boldsymbol{\lambda}_{i}^{k}\right)\right)-\frac{b\left(\epsilon\right)}{l}\right)$ with

[TABLE]

Remark 4: It is worth to note that the Laplacian matrix $L$ cannot be locally computed such that Eq. $\left(\ref{PDGD1}\right)$ is not distributed. We introduce a new variable to scale the Lagrangian multiplier $\boldsymbol{\mu}_{i}$ as $\boldsymbol{\gamma}_{i}=\sqrt{\boldsymbol{L}}\boldsymbol{\mu}_{i}$. Thus, the PDGD in $\left(\ref{PDGD1}\right)$ is equivalent to

[TABLE]

Notice that from a global perspective, the introduced $\boldsymbol{\gamma}^{0}=\left[\boldsymbol{\gamma}_{1}^{0},\boldsymbol{\gamma}_{2}^{0},\cdots,\boldsymbol{\gamma}_{l}^{0}\right]^{\mathrm{T}}=\sqrt{\boldsymbol{L}}\boldsymbol{\Lambda}^{0}$ requires the initial values to satisfy $\sum_{i=1}^{l}\boldsymbol{\gamma}_{i}^{0}=0$.

IV-D Contraction of the PDGD

Inspired by the Riemannian geometry, some sufficient conditions guaranteeing the exponential convergence of the PDGD in $\left(\ref{PDGD2}\right)$ are given in the following theorem.

Theorem 1: Under Assumptions 2 and 4, the PDGD in $\left(\ref{PDGD2}\right)$ has the exponential convergence rate $\tau=\sqrt{1-\rho\alpha_{i}\beta_{i}d_{i}}$ with Riemannian metric

[TABLE]

if

[TABLE]

where $\underline{\sigma}$ and $\overline{\sigma}$ derived from $\left(\ref{assumption2}\right)$ and $\left(\ref{assumption4}\right)$ are two proper scalars satisfying $0<\underline{\sigma}\leq\overline{\sigma}<\infty$, $0<\rho<\min\left(1,\frac{1}{\alpha_{i}\beta_{i}d_{i}}\right)$ is a user-defined scalar parameter and $\boldsymbol{I}$ is with appropriate dimension.

Proof. Stack the Lagrangian multipliers $\boldsymbol{\lambda}_{i}$ and $\boldsymbol{\gamma}_{i}$ into a vector $y=\left[\boldsymbol{\lambda}_{i}^{\mathrm{T}},\boldsymbol{\gamma}_{i}^{\mathrm{T}}\right]^{\mathrm{T}}$. The optimal vector can be similarly defined as $y^{*}=\left[\left(\boldsymbol{\lambda}_{i}^{*}\right)^{\mathrm{T}},\left(\boldsymbol{\gamma}_{i}^{*}\right)^{\mathrm{T}}\right]^{\mathrm{T}}$. According to the updating dynamics of the primal-dual gradient in $\left(\ref{PDGD2}\right)$, we can obtain the differential dynamics of PDGD as follows

[TABLE]

where $\Theta_{i}=G_{i}H_{i}^{-1}G_{i}^{\mathrm{T}}$ and

[TABLE]

The second equality of $\left(\ref{diffDynPDGD}\right)$ follows from the condition in $\left(\ref{u_i}\right)$ and the KKT condition in $\left(\ref{dualKKT:1}\right)$ under Assumption 2. By constructing the Riemannian metric as $\left(\ref{RieMatrix}\right)$, the difference of the Riemannian energy between the adjacent updates can be written as

[TABLE]

where

[TABLE]

with

[TABLE]

Therefore, to prove $\left(\ref{contractionTheoryDef}\right)$, it suffices to prove that $\Pi\leq\left(\tau^{2}-1\right)M$. Letting $\Upsilon=\left(\tau^{2}-1\right)M-\Pi$, we can obtain

[TABLE]

with

[TABLE]

Resorting to the Schur’s complement, to prove $\Upsilon\geq 0$, it is sufficient to prove $\Upsilon_{3}>0$ and $\Upsilon_{1}-\Upsilon_{2}\Upsilon_{3}^{-1}\Upsilon_{2}^{\mathrm{T}}\geq 0$. By introducing a user-defined parameter $0<\rho<\min\left(1,\frac{1}{\alpha_{i}\beta_{i}d_{i}}\right)$ such that $\tau=\sqrt{1-\rho\alpha_{i}\beta_{i}d_{i}}$, one can get

[TABLE]

By using $\alpha_{i}<\frac{\left(1-\rho\right)\underline{\sigma}}{2\overline{\sigma}^{2}-\left(3+\rho\right)\underline{\sigma}\beta_{i}d_{i}}$ and $\beta_{i}<\frac{\overline{\sigma}}{2d_{i}}$, one have

[TABLE]

In light of Assumptions $2$ and $4$, it suffices to have

[TABLE]

where $\underline{\sigma}$ and $\overline{\sigma}$ are two scalars satisfying $0<\underline{\sigma}\leq\overline{\sigma}<\infty$, which are derived from $\left(\ref{assumption2}\right)$ and $\left(\ref{assumption4}\right)$ appropriately. Thus, we can obtain the last inequality of $\left(\ref{Upsilon1}\right)$ following from

[TABLE]

Therefore, we have $\Upsilon_{1}-\Upsilon_{2}^{\mathrm{T}}\Upsilon_{3}^{-1}\Upsilon_{2}\geq 0$, which completes the proof. $\hfill\blacksquare$

IV-E Convergence Analysis

In light of the contraction of the primal-dual gradient updating dynamics, we analyze the convergence of the distributed primal-dual gradient optimization algorithm. The theoretical results are summarized in the following theorem.

Theorem 2: Suppose Assumptions 1 and 3 hold. For the distributed primal-dual gradient optimization characterized in $\left(\ref{PDGD2}\right)$, it holds for all $k\geq 1$ that

[TABLE]

where $0<\mathcal{Q}_{i}<\infty$ is a scalar depended on the initial and optimal values of $\boldsymbol{\lambda}$.

Proof. By using the triangle inequality $\left\|a+b\right\|\leq\left\|a\right\|+\left\|b\right\|$, we can obtain

[TABLE]

It is seen that the condition in $\left(\ref{u_i}\right)$ is exact which implies $\left\|\boldsymbol{u}_{i}^{k}-\left(\boldsymbol{u}_{i}^{k}\right)^{*}\right\|=0$. Thereafter, according to the $m_{J}$-strong convexity in $\left(\ref{assumption1}\right)$, we have

[TABLE]

According to the KKT conditions $\nabla_{{\boldsymbol{u}}_{i}^{k}}J_{i}\left(x_{i},\left({\boldsymbol{u}}_{i}^{k}\right)^{*}\right)+G_{i}^{\mathrm{T}}\boldsymbol{\lambda}_{i}^{k}=0$ and $\nabla_{{\boldsymbol{u}}_{i}^{*}}J_{i}\left(x_{i},{\boldsymbol{u}}_{i}^{*}\right)+G_{i}^{\mathrm{T}}\boldsymbol{\lambda}_{i}^{*}=0$, we can obtain

[TABLE]

with $0<\mathcal{Q}_{i}<\infty$ being a proper scalar depended on the initial and optimal values of $\boldsymbol{\lambda}$, where the last inequality follows from Assumption 3 and the contraction analysis on the PDGD in Theorem 1. The condition in $\left(\ref{convergenceAnalysisTheorem}\right)$ can be established, and the proof is completed. $\hfill\blacksquare$

IV-F The Stopping Criterion

In this paper, we adopt the constraint tightening approach to provide the capability of early termination for the DMPC optimization with guaranteeing convergence. The following definition is given to characterize the inexact solution resulted from the early termination. Thereafter, to fulfill the requirement on premature termination, the stopping criterion of the DMPC optimization is developed.

Definition 3: [9] For a given $\epsilon>0$, if it holds that

[TABLE]

where $\boldsymbol{u}_{i}\in\tilde{\mathcal{U}}_{i}\left(x_{i}\right)$, $\boldsymbol{u}_{i}$ is called $\epsilon$-feasible solution of the optimization problem in $\left(\ref{erMPC}\right)$.

According to the convergence analysis results on PDGD in $\left(\ref{PDGD2}\right)$, we can establish the stopping criterion in the following theorem.

Theorem 3: Given the initial parameters $\epsilon$ and $\gamma_{i}^{0}$, the control input sequence $\tilde{u}_{i}^{\bar{k}}$ is $\epsilon$-feasible under the distributed primal-dual gradient optimization in $\left(\ref{PDGD2}\right)$, where the superscript $\bar{k}$ represents the stopping iteration defined as

[TABLE]

with $\tau$ given in Theorem 1, and $\boldsymbol{\gamma}_{i}^{0}$ being the initial value of $\boldsymbol{\gamma}_{i}$.

Proof. In light of the PDGD in $\left(\ref{PDGD2}\right)$ and the KKT conditions in $\left(\ref{KKT}\right)$, it holds that

[TABLE]

Therefore, we can rewrite $\left(\ref{stoppingCriterion1}\right)$ in a distributed form as follows

[TABLE]

where the superscript $\bar{k}$ stands for the stopping iteration. According to the contraction analysis results in Theorem 1, one can get

[TABLE]

Following from $\left(\ref{stoppingCriterion3}\right)$ and substituting $\gamma_{i}^{*}=0$ into $\left(\ref{stoppingCriterion4}\right)$, Eq. $\left(\ref{stoppingCriterion}\right)$ can be derived, which completes the proof.$\hfill\blacksquare$

Remark 5: It can be seen that the stopping iteration $\overline{k}$ in $\left(\ref{stoppingCriterion}\right)$ depends on the initial values of $\gamma_{i}$ and the user-defined parameter $\epsilon$. For a given $\epsilon$, to initialize the appropriate values of $\gamma_{i}$ directly impacts the stopping iteration.

V Algorithm Description

V-A Distributed Primal-Dual Gradient Algorithm

In terms of Theorem 3, the distributed primal-dual gradient algorithm terminates at the iteration $k$, which implies the $\epsilon$-feasible solution is obtained. Thus, the distributed primal-dual gradient algorithm based on Laplacian consensus is described in Algorithm 1.

V-B DMPC Algorithm

The overall DMPC algorithm adopting the proposed distributed primal-dual gradient algorithm is formulated in the following Algorithm 2.

VI Recursive Feasibility and Stability

Under the proposed distributed primal-dual gradient algorithm, the recursive feasibility and stability of the DMPC are analyzed. The following theorem summarized the theoretical analysis results.

Theorem 4: Suppose the DMPC optimization problem in $\left(\ref{erMPC}\right)$ is feasible at time instant $t$ for each subsystem $i$ with $i=1,2,\cdots,l$. Under Assumption 3, the following results hold: i) the DMPC optimization problem in $\left(\ref{erMPC}\right)$ has a feasible solution at time instant $t+1$, ii) the state trajectory of the closed-loop system in $\left(\ref{sysForm}\right)$ enters the terminal set $\Omega_{i}$ in finite time and remains in it.

Prrof. The proof of this theorem consists of two parts known as feasibility and stability analysis, respectively.

Feasibility: Recalling the definition on $\epsilon$-feasible solution in $\left(\ref{stoppingCriterion1}\right)$ and the stopping criterion in $\left(\ref{stoppingCriterion3}\right)$, it is shown that a $\epsilon$-feasible solution $\boldsymbol{u}_{i}$ satisfies that

[TABLE]

with $s\in\mathbb{Z}_{\left[0,N-1\right]}$ according to the coupled constraint in $\left(\ref{coupledConstraint}\right)$ where $\left\{u_{i}\left(\cdot|t\right)\right\}=\left\{u_{i}\left(0|t\right),u_{i}\left(1|t\right),\cdots,u_{i}\left(N-1|t\right)\right\}$ and $\left\{x_{i}\left(\cdot|t\right)\right\}=\left\{x_{i}\left(0|t\right),x_{i}\left(1|t\right),\cdots,x_{i}\left(N|t\right)\right\}$. Define a feasible solution at time instant $t+1$ for subsystem $i$ as

[TABLE]

and the corresponding state sequence can be expressed by

[TABLE]

where $A_{K}$ is defined in $\left(\ref{assumption3}\right)$. Thus, we can obtain

[TABLE]

with $s\in\mathbb{Z}_{\left[0,N-2\right]}$. The last inequality follows from the coupled constraint in the terminal set $\left(\ref{localCoupleConstraint}\right)$. Moreover, as the feasibility of the solution at time instant $t$, it holds that $u_{i}\left(s|t\right)\in\mathcal{U}_{i}$ and $x_{i}\left(s|t\right)\in\mathcal{X}_{i}$ for $s\in\mathbb{Z}_{\left[0,N-1\right]}$ according to the constraint in $\left(\ref{erMPC:U}\right)$. It implies that $\tilde{u}_{i}\left(s+1|t+1\right)\in\mathcal{U}_{i}$ and $\tilde{x}_{i}\left(s+1|t+1\right)\in\mathcal{X}_{i}$ for $s\in\mathbb{Z}_{\left[0,N-2\right]}$. On the other hand, the local control law $\tilde{u}_{i}\left(N|t+1\right)=Kx_{i}\left(N|t\right)\in\mathcal{U}_{i}$ in terms of $x_{i}\left(N|t\right)\in\Omega_{i}$. According to the aforementioned analysis results, $\left\{\tilde{u}_{i}\left(\cdot|t+1\right)\right\}$ is a feasible solution to the DMPC optimization problem in $\left(\ref{erMPC}\right)$ at the successor time instant $t+1$.

Stability: The stability is analyzed by taking the objective function in $\left(\ref{erMPC:ObjFcn}\right)$ as a Lyapunov candidate which refers to $V\left(x_{i}\left(t\right)\right):=\sum_{i=1}^{l}J_{i}\left(x_{i}\left(t\right),\boldsymbol{u}_{i}^{*}\left(\cdot|t\right)\right)$. At time instant $t$, denote $\tilde{k}\left(t\right)$ as the iteration $k$ when the stopping criterion in $\left(\ref{stoppingCriterion3}\right)$ is satisfied. In light of Assumption 3, it holds that

[TABLE]

As the inexactness of the optimization solver induced by the premature termination, one can have

[TABLE]

It can be seen that $F\left(x_{i}\left(t\right),u_{i}\left(t\right)\right)\rightarrow 0$ as $t\rightarrow 0$, which implies there exist finite time instants before the state trajectory enters the terminal set. Thereafter, according to the local controller in Assumption 3, the state trajectory of the closed-loop system will remain in the terminal set, which completes the proof. $\hfill\blacksquare$

VII Simulation

A numerical simulation is given to verify the performance of the presented DMPC approach. Consider a coupled double-water tank systems containing four subsystems [31, 34]. The control objective is to regulate the water level towards the given reference by means of the input flow which is subject to a globally coupled rate constraint. As shown in Fig.1, denote the sampled water level as $h_{i,1}$ and $h_{i,2}$, the input flow $q_{i}$ for each subsystem $i$ with $i=1,2,3,4$. Given the reference water level as $r_{i,1}=1$ and $r_{i,1}=0.64$, and the steady-state input flow $\overline{q}_{i}=0.3$, formulate each subsystem $i$ as the following discrete-time linear dynamical system

[TABLE]

where $x_{i,1}=h_{i,1}-r_{i,1}$, $x_{i,2}=h_{i,2}-r_{i,2}$, $u_{i}=q_{i}-\overline{q}_{i}$, and

[TABLE]

Define the objective function for the DMPC optimization problem as

[TABLE]

where $N=8$, $Q_{i}=10\boldsymbol{I}$, $R_{i}=1$,

[TABLE]

and the local control gain matrix $K_{i}=\left[-1.7916\,-0.7337\right]$. Taking the physical constraints into account, each subsystem should consider the local constraints on states as $\mathcal{X}_{i}=\left\{x_{i}\in\mathbb{R}^{2}|-1\leq x_{i,1}\leq 1,-0.64\leq x_{i,2}\leq 0.64\right\}$ and input $\mathcal{U}_{i}=\left\{u_{i}\in\mathbb{R}|-0.3\leq u_{i}\leq 0.3\right\}$, respectively. As the presence of the bound on total input flow rate which is supposed to be 2 in this paper, a global constraint is involved in the DMPC optimization problem as $\sum_{i=1}^{4}q_{i}\leq 2$, which means $\sum_{i=1}^{4}\Phi_{i}^{u}u_{i}\leq 1$ with $\Phi_{i}^{u}=1.25$. In addition, the adjacent matrix and the corresponding Laplacian matrix of the connection network are given as follows

[TABLE]

The initial states of each subsystem are $x_{1}\left(0\right)=[-0.2264\,-0.3981]^{\mathrm{T}}$, $x_{2}\left(0\right)=[0.4520\,-0.3689]^{\mathrm{T}}$, $x_{3}\left(0\right)=[-0.5311\,-0.2828]^{\mathrm{T}}$ and $x_{4}\left(0\right)=[-0.4397\,-0.4897]^{\mathrm{T}}$, respectively. Define the accuracy for coupled constraints as $\epsilon=0.1$, and $\boldsymbol{\gamma}_{i}^{0}$ is initialized accordingly. In terms of Theorem 1, we choose $\alpha=0.2$ and $\beta=0.19$.

The simulation results are shown in Fig. 2-Fig. 5. Fig. 2 gives the control input of each subsystems. Meanwhile, the total control input demonstrates the satisfaction of the global coupled constraints. In Fig. 3, the states of each subsystem are plotted, which shows that the control objective can be achieved. To illustrate the effectiveness of the user-defined parameters introduced in the proposed approach, we give the total control input under different $\rho$ and $\epsilon$ in Fig. 4 and Fig. 5, respectively. It is seen that the changing rate of the control input is increasing as $\rho$ gets bigger. Moreover, a smaller $\epsilon$ leads to more conservative of control input.

VIII Conclusion

In this paper, a novel distributed model predictive control (DMPC) approach is studied for a group of discrete-time linear systems taken the global coupled constraints into account. The DMPC optimization problem is transformed into a dual problem involving all subsystems, which is solved in the framework of Laplacian consensus by using the primal-dual gradient optimization in a fully distributed manner. To reduce the computational burden, a tightening constraint concerning the global coupled constraint is constructed to allow premature termination with guaranteeing the convergence of the optimization. It is worth seeing that the local copies of the Lagrangian multipliers need not reach consensus but within some specified bounds owing to the constraint tightening method. Furthermore, the convergence of the primal-dual gradient optimization is first rigorously analyzed by means of contraction theory in the context of discrete-time nonlinear dynamics. The recursive feasibility and stability of the closed-loop system are established under the inexact solver with rational assumptions. The performance of the proposed approach is demonstrated by a numerical simulation.