Tracking-ADMM for Distributed Constraint-Coupled Optimization

TL;DR

This paper introduces a distributed optimization algorithm combining dynamic average consensus with ADMM to solve constraint-coupled problems, ensuring convergence to optimal solutions in a networked setting.

Contribution

It develops a novel fully distributed scheme that tracks constraint violations dynamically within ADMM, enhancing distributed optimization methods.

Findings

Proves convergence of primal and dual estimates to optimal solutions.

Demonstrates the effectiveness through a numerical case study.

Ensures consensus and asymptotic optimality in a distributed setting.

Abstract

We consider constraint-coupled optimization problems in which agents of a network aim to cooperatively minimize the sum of local objective functions subject to individual constraints and a common linear coupling constraint. We propose a novel optimization algorithm that embeds a dynamic average consensus protocol in the parallel Alternating Direction Method of Multipliers (ADMM) to design a fully distributed scheme for the considered set-up. The dynamic average mechanism allows agents to track the time-varying coupling constraint violation (at the current solution estimates). The tracked version of the constraint violation is then used to update local dual variables in a consensus-based scheme mimicking a parallel ADMM step. Under convexity, we prove that all limit points of the agents' primal solution estimates form an optimal solution of the constraint-coupled (primal) problem. The result is proved by means of a Lyapunov-based analysis simultaneously showing consensus of the dual estimates to a dual optimal solution, convergence of the tracking scheme and asymptotic optimality of primal iterates. A numerical study on optimal charging schedule of plug-in electric vehicles corroborates the theoretical results.

Equations

$$\begin{array}{rl}\underset{x_1,\dots ,x_N}{\min }& \sum _{i=1}^N{f}_i(x_i)\\ \text{subject to:}& \sum _{i=1}^N{A}_i{x}_i=b\\ & x_i\in X_{i}i=1,\dots ,N,\end{array}$$

$$L(\boldsymbol{x},\lambda)=\sum_{i=1}^{N}f_{i}(x_{i})+\lambda^{\top}\bigg{(}\sum_{i=1}^{N}A_{i}x_{i}-b\bigg{)}$$

$$\underset{\lambda \in {\mathbb{R}}^p}{\max }\underset{\mathbf{x}\in X}{\min }L(\mathbf{x},\lambda )=\underset{\lambda \in {\mathbb{R}}^p}{\max }\sum _{i=1}^N{\phi }_i(\lambda ),$$

$${\phi }_i(\lambda )=\underset{x_i\in X_i}{\min }{f}_i(x_i)+{\lambda }^{\top }(A_i{x}_i-b_i),$$

$$\displaystyle\begin{split}&x_{i,k+1}\in\operatorname*{argmin}_{x_{i}\in X_{i}}\Big{\{}f_{i}(x_{i})+\lambda_{k}^{\top}A_{i}x_{i}\\ &\hskip 90.00014pt+\frac{c}{2}\|A_{i}x_{i}-A_{i}x_{i,k}+d_{k}\|^{2}\Big{\}},\end{split}$$

$$d_{k+1}=\frac{1}{N}\sum _{i=1}^N\left(A_i{x}_{i,k+1}-b_i\right)$$

$${\lambda }_{k+1}={\lambda }_k+c{d}_{k+1},$$

$${\lambda }_{i,k+1}=\sum _{j\in {\mathcal{N}}_i}{w}_{ij}{\lambda }_{j,k}+c{d}_{k+1},$$

$$d_{i,k+1}=\sum _{j\in {\mathcal{N}}_i}{w}_{ij}{d}_{j,k}+(A_i{x}_{i,k+1}-b_i)-(A_i{x}_{i,k}-b_i),$$

$${\mathbf{x}}_k$$

$${\mathbf{d}}_k$$

$${\mathbf{\lambda }}_k$$

$$\displaystyle\begin{split}&\boldsymbol{x}_{k+1}\in\operatorname*{argmin}_{\boldsymbol{x}\in X}\Big{\{}f(\boldsymbol{x})+(W\boldsymbol{\lambda}_{k})^{\top}A_{d}\boldsymbol{x}\\ &\hskip 80.00012pt+\frac{c}{2}\|A_{d}\boldsymbol{x}-A_{d}\boldsymbol{x}_{k}+W\boldsymbol{d}_{k}\|^{2}\Big{\}}\end{split}$$

$${\mathbf{d}}_{k+1}=W{\mathbf{d}}_k+A_d{\mathbf{x}}_{k+1}-A_d{\mathbf{x}}_k$$

$${\mathbf{\lambda }}_{k+1}=W{\mathbf{\lambda }}_k+c{\mathbf{d}}_{k+1},$$

$${\delta }_{i,k}^{\prime }={\sum }_{j\in {\mathcal{N}}_i}{w}_{ij}{d}_{j,k}\text{and}{\ell }_{i,k}^{\prime }={\sum }_{j\in {\mathcal{N}}_i}{w}_{ij}{\lambda }_{j,k}$$

$${\delta }_{i,k}={\sum }_{j\in {\mathcal{N}}_i}{w}_{ij}{\delta }_{j,k}^{\prime }\text{and}{\ell }_{i,k}={\sum }_{j\in {\mathcal{N}}_i}{w}_{ij}{\ell }_{j,k}^{\prime }$$

$${\mathbf{\delta }}_k$$

$${\mathbf{\ell }}_k$$

$${\bar{d}}_k=\frac{1}{N}\sum _{i=1}^N{d}_{i,k}\text{and}{\bar{\lambda }}_k=\frac{1}{N}\sum _{i=1}^N{\lambda }_{i,k},$$

$${\bar{d}}_k$$

$${\bar{\lambda }}_{k+1}$$

$${\bar{\mathbf{d}}}_k=\text{mathbbm}{1}_N\otimes {\bar{d}}_k,$$

$${\bar{\mathbf{\lambda }}}_k=\text{mathbbm}{1}_N\otimes {\bar{\lambda }}_k,$$

$${\mathbf{e}}_k^{\mathbf{d}}$$

$${\mathbf{e}}_k^{\mathbf{\lambda }}$$

$${\mathbf{z}}_k=A_d{\mathbf{x}}_k-{\bar{\mathbf{d}}}_k.$$

$${\mathbf{e}}_{k+1}^{\mathbf{d}}$$

$$=W{\mathbf{d}}_k+A_d{\mathbf{x}}_{k+1}-A_d{\mathbf{x}}_k-{\bar{\mathbf{d}}}_{k+1}\pm {\bar{\mathbf{d}}}_k$$

$$\stackrel{(\mathrm{a})}{=}W{\mathbf{d}}_k-{\bar{\mathbf{d}}}_k+[{\mathbf{z}}_{k+1}-{\mathbf{z}}_k]$$

$$\stackrel{(\mathrm{b})}{=}W{\mathbf{e}}_k^{\mathbf{d}}+[{\mathbf{z}}_{k+1}-{\mathbf{z}}_k],$$

$$W{\bar{\mathbf{d}}}_k$$

$$\stackrel{(\mathrm{c})}{=}(\mathcal{W}\text{mathbbm}{1}_N\otimes I_p{\bar{d}}_k)$$

$$\stackrel{(\mathrm{d})}{=}(\text{mathbbm}{1}_N\otimes {\bar{d}}_k)={\bar{\mathbf{d}}}_k,$$

$$(\mathcal{O}\otimes I_p){\mathbf{e}}_k^{\mathbf{d}}$$

$$=\text{mathbbm}{1}_N\otimes \frac{1}{N}\sum _{i=1}^N[d_{i,k}-{\bar{d}}_k]$$

$$=\text{mathbbm}{1}_N\otimes \frac{1}{N}[N{\bar{d}}_k-N{\bar{d}}_k]=0,$$

$${\mathbf{e}}_{k+1}^{\mathbf{d}}=\tilde{W}{\mathbf{e}}_k^{\mathbf{d}}+[{\mathbf{z}}_{k+1}-{\mathbf{z}}_k],$$

Full text

††journal: Automatica

1 Introduction

Given the increasing penetration of network-connected devices and cloud-computing services in the industry sector, systems are becoming more and more large-scale, composed by different interacting agents, and complex. Efficiently operating these systems is thus becoming more challenging: the computational burden associated to the computation of an efficient strategy for the entire system may be prohibitive for a single processing unit and privacy issues may arise if the different agents are forced to disclose sensitive information to the entity in charge of solving the problem. This challenge has been taken up by the control community, which has directed its efforts in devising distributed algorithms for the resolution of such optimal control problems.

In this paper we investigate a set-up that is relevant to several network control applications. We consider agents which can communicate only with their neighbors in a network. Each agent knows its local cost function and constraints, that depend on its individual decision variables. Moreover, the agents’ decision variables are coupled by a linear constraint. Each agent knows only how its local decision variables affect the coupling constraint. The agents collectively aim at minimizing the sum of the local cost functions subject to both local and coupling constraints. The presence of this coupling element makes the problem solution challenging, especially in the considered distributed context in which no central authority (communicating with all the agents) is present.

Even though constraint-coupled problems arise more naturally in practical applications the main effort in distributed optimization, until recently, has been devoted to solve optimization problems in which there is no coupling constraint but the agents are required to agree on a common decision. For this set-up, consensus methods based on (sub)gradient iterations and proximal operators (see [1, 2, 3, 4, 5, 6, 7, 8]), distributed algorithms based on duality (see [9, 10]), and distributed approaches based on the Alternating Direction Method Multipliers (ADMM) (see [11, 12, 13, 14, 15, 16]) have been proposed. We refer the reader to [17, 18] for an in-depth discussion on ADMM. Recently, distributed gradient schemes for problems with common decision variables have also been combined with a technique known as dynamic average consensus (firstly proposed in [19, 20]) to achieve linear convergence rate with constant step-size at the expense of imposing more restrictive assumptions on the cost functions, see [21, 22, 23, 24, 25, 26]. Applying the above mentioned techniques to the constraint-coupled set-up considered in this paper is possible but it would require to define as common decision vector the collection of all local decision vectors. Agents would then have to store and update the entire solution estimate of the problem rather than the portion associated with their decision variables only, and would have all to know the global coupling constraint. These major drawbacks hamper the applicability of those methods to practical applications, thus calling for novel and efficient strategies that takes advantage of the network structure of the problem at hand.

Recently developed approaches tackling constraint-coupled problems directly are mostly based on Lagrangian duality. In [27] a problem with sparse coupling constraints is considered and a distributed dual gradient method achieving linear convergence rate is proposed. A similar structured problem is considered in [28] where a tailored algorithm exploiting such structure is proved to converge to the optimal solution of the problem. In [29, 30], primal-dual approaches are proposed but they need a diminishing step-size to achieve convergence. In [31, 32] distributed dual subgradient algorithms are proposed, in [33] the dual problem is tackled by means of consensus-ADMM and proximal operators, while an alternative approach based on successive duality steps has been investigated in [34]; however, all these algorithms typically exhibit a slow convergence rate for the local decision variables. The tracking mechanism has been also employed in [35] to solve constraint-coupled problems based on an augmented Lagrangian approach for a continuous time setting, but the considered set-up does not allow for nonsmooth costs and local constraints. Finally, in the very recent paper [36] an idea similar to the one proposed in this paper is introduced. However, the approach in [36] requires the agents to perform multiple communication rounds and the number of rounds has to be carefully tuned to bound a non-vanishing steady state error between the asymptotic algorithmic solution and the optimal one.

The contributions of this paper are as follows. We propose a novel, fully distributed optimization algorithm to solve constraint-coupled problems over networks by means of an ADMM-based approach. Differently from distributed ADMM schemes for problems with common decision variables, we design our Tracking-ADMM distributed algorithm by embedding a tracking mechanism into the parallel ADMM designed for constraint-coupled problems.

The resulting algorithm enjoys the following appealing features: (i) no parameter tuning is needed, in fact our algorithm works for all the (constant) choices of a penalty parameter and no other coefficients are necessary; (ii) agents solve optimization problems depending on their local (few) decision variables and asymptotically compute only their portion of an optimal (hence feasible) solution to the given problem; (iii) the local estimate of the coupling constraint violation gives each agent a local assessment on the amount of infeasibility, which can be useful, e.g., in designing distributed (receding horizon) control schemes.

The convergence proof of our Tracking-ADMM for constraint-coupled problems relies on control systems analysis tools. By explicitly relying on Lyapunov theory for linear systems, we are able to find a constrained solution to a suitable discrete Lyapunov equation that leads to an aggregate descent condition allowing us to prove: (i) the (exact) convergence of the dual variables to a dual optimal solution and (ii) that any limit point of the primal sequences is an optimal solution of the original constraint-coupled problem. This novel approach allows us to derive a clean and elegant proof of the asymptotic optimality of our algorithm.

The rest of the paper is organized as follows. In Section 2 we present the problem set-up and in Section 2.2 we revise the ADMM algorithm. In Section 3 we introduce our novel Tracking-ADMM distributed algorithm and analyze its convergence properties, discussing the main steps of the proof. In Section 4 we apply our algorithm on a realistic application related to the optimal charging schedule for a fleet of electric vehicles. In Section 5 we draw some conclusions and finally, in A we report the proofs of all the results stated in the body of the paper.

Notation

The vector in $\mathbb{R}^{n}$ containing all ones is denoted by $\mathbbm{1}_{n}$. The identity matrix of order $n$ is denoted by $I_{n}$. The Kronecker product is denoted by $\otimes$. For a matrix $S$ we write $S^{\top}$ to denote its transpose, $S\succ 0$ to denote that $S$ is positive definite, $\rho(S)$ to denote the spectral radius of $S$, and $\|S\|$ for its spectral norm. We write $\operatorname{blkdiag}(S_{1},\dots,S_{n})$ to refer to the block-diagonal matrix with $S_{1},\dots,S_{n}$ as blocks. For a vector $v$, $\|v\|$ is the Euclidean norm of $v$, and, for any matrix $S\succ 0$, $\|v\|_{S}$ is the weighted norm of $v$, i.e., $\|v\|_{S}^{2}=v^{\top}Sv$.

2 Constraint-Coupled Optimization

In this section we introduce the optimization set-up and recall some preliminaries about the Alternating Direction Method of Multipliers (ADMM).

2.1 Optimization Problem and Assumptions

We consider a system composed of $N$ agents which are willing to cooperatively solve an optimization program formulated over the entire system. Each agent has to set its local decision variables $x_{i}\in\mathbb{R}^{n_{i}}$ so as to minimize the sum of local objective functions $f_{i}:\mathbb{R}^{n_{i}}\rightarrow\mathbb{R}$, while satisfying local constraints $X_{i}\subset\mathbb{R}^{n_{i}}$ as well as a linear constraint that couples the decisions of all the agents. Formally, we address the following mathematical program

[TABLE]

where $A_{i}\in\mathbb{R}^{p\times n_{i}}$ and $b\in\mathbb{R}^{p}$ specify the coupling constraint.

We impose the following regularity conditions on $\mathcal{P}$.

Assumption 1** (Convexity and compactness).**

For all $i=1,\dots,N$, the function $f_{i}$ is convex and the set $X_{i}$ is convex and compact. ∎

Let $\boldsymbol{x}=[x_{1}^{\top}\,\cdots\,x_{N}^{\top}]^{\top}$, consider a vector $\lambda\in\mathbb{R}^{p}$ of Lagrange multipliers and define

[TABLE]

the Lagrangian function obtained by dualizing the coupling constraint $\sum_{i=1}^{N}A_{i}x_{i}=b$. Then, the dual problem of $\mathcal{P}$ is

[TABLE]

where $X=X_{1}\times\cdots\times X_{N}$, and the $i$-th contribution $\varphi_{i}$ is defined as

[TABLE]

the vectors $b_{1},\dots,b_{N}$ being such that $\sum_{i=1}^{N}b_{i}=b$.

The next assumption ensures that $\mathcal{P}$ and $\mathcal{D}$ are well-posed.

Assumption 2** (Existence of optimal solutions).**

Problem $\mathcal{P}$ admits an optimal solution $\boldsymbol{x}^{\star}=[x_{1}^{\star}\!{}^{\top}\,\cdots\,x_{N}^{\star}\!{}^{\top}]^{\top}$ and problem $\mathcal{D}$ admits an optimal solution $\lambda^{\star}$. ∎

Next, we revise the popular ADMM algorithm which provides an effective way to solve $\mathcal{P}$ by splitting the computation over $N$ agents coordinated by a central unit.

2.2 The ADMM Parallel Algorithm

A version of the ADMM algorithm specifically tailored to problem $\mathcal{P}$ is presented in [17, pag. 254, eq. (4.75)] and is reported here with our notation for the reader’s convenience. Given initial values $x_{i,0}\in X_{i}$, $d_{0}=(1/N)\sum_{i=1}^{N}(A_{i}x_{i,0}-b_{i})$, and $\lambda_{0}\in\mathbb{R}^{p}$, at each iteration $k\geq 0$, a set of $N$ agents and a central unit perform the following two steps. First, all agents, $i=1,\dots,N$, compute (in parallel) a minimizer of the following optimization problem

[TABLE]

The evolution of (3) is analyzed in [17, pagg. 254-256] and its convergence property is reported below.

Proposition** ([17, Proposition 4.2]).**

Let Assumptions 1 and 2 hold. Then, any limit point of the primal sequence ${\{[x_{1,k}^{\top}\,\cdots\,x_{N,k}^{\top}]^{\top}\}}_{k\geq 0}$ generated by (3a), is an optimal solution of $\mathcal{P}$, and the dual sequence ${\{\lambda_{k}\}}_{k\geq 0}$, generated by (3c), converges to an optimal solution of $\mathcal{D}$. ∎

Notice that, since all limit points of the primal sequence are optimal, they are necessarily feasible for the coupling constraint. It thus follows that the sequence ${\{d_{k}\}}_{k\geq 0}$ generated by (3b) converges to zero. Moreover, we shall stress that no requirement on the penalty parameter $c$ is necessary for the convergence result to hold.

Finally, note that the algorithm described by (3) requires a central unit to compute (3b) and (3c). This hampers the applicability of ADMM to a distributed computation framework, where agents communicate only with neighbors according to a (typically sparse) communication graph. Schemes requiring a central unit are usually referred to as parallel schemes. For this reason in the rest of the paper we will refer to (3) as the parallel ADMM.

3 Tracking-ADMM Distributed Algorithm

In this section we propose our novel distributed optimization algorithm for constraint-coupled problems. We first introduce the distributed framework and then describe the algorithm along with its convergence properties.

3.1 Distributed Computation Framework

Assume that, at each iteration $k$, the $N$ agents communicate according to a graph $\mathcal{G}=(\mathcal{V},\mathcal{E})$, where $\mathcal{V}=\{1,\dots,N\}$ is the set of nodes, representing the agents, and $\mathcal{E}\subseteq\mathcal{V}\times\mathcal{V}$ is the set of edges, representing the communication links. The presence of edge $(i,j)$ in $\mathcal{E}$ models the fact that agent $i$ receives information from agent $j$. We assume that the communication graph does not change across iterations and, consequently, $\mathcal{E}$ does not depend on the iteration index $k$. We denote by $\mathcal{N}_{i}=\{j\in\mathcal{V}\mid(i,j)\in\mathcal{E}\}$ the set of neighbors of agent $i$ in $\mathcal{G}$, assuming that $(i,i)\in\mathcal{E}$ for all $i=1,\dots,N$. We impose the following connectivity property on $\mathcal{G}$.

Assumption 3** (Connectivity).**

The graph $\mathcal{G}$ is undirected and connected, i.e., $(i,j)\in\mathcal{E}$ if and only if $(j,i)\in\mathcal{E}$ and for every pair of vertices in $\mathcal{V}$ there exists a path of edges in $\mathcal{E}$ that connects them. ∎

Each edge $(i,j)\in\mathcal{E}$ has an associated weight $w_{ij}$, which measures how much agent $i$ values the information received by agent $j$. For those $(i,j)\notin\mathcal{E}$ we set $w_{ij}=0$, which models the fact that agent $i$ does not receive any information from agent $j$. We impose the following assumption on the network weights.

Assumption 4** (Balanced information exchange).**

For all $i,j=1,\dots,N$, $w_{ij}\in[0,1)$ and $w_{ij}=w_{ji}$. Furthermore

• •

$\sum_{i=1}^{N}w_{ij}=1$* for all $j=1,\dots,N$,*

• •

$\sum_{j=1}^{N}w_{ij}=1$* for all $i=1,\dots,N$,*

and $w_{ij}>0$ if and only if $(i,j)\in\mathcal{E}$. ∎

Let $\mathcal{W}\in\mathbb{R}^{N\times N}$ be the matrix whose $(i,j)$-th entry is $w_{ij}$, often referred to as the consensus matrix. Assumption 4 translates into requiring $\mathcal{W}$ to be symmetric and doubly stochastic, i.e., $\mathcal{W}=\mathcal{W}^{\top}$ and $\mathcal{W}\mathbbm{1}_{N}=\mathcal{W}^{\top}\mathbbm{1}_{N}=\mathbbm{1}_{N}$. We should point out that Assumptions 3 and 4 are common in the consensus-based distributed optimization literature, see, e.g., [2, 3].

Finally, we impose the following additional assumption on the consensus matrix.

Assumption 5**.**

$\mathcal{W}$* is positive semidefinite. ∎*

Note that this assumption is not too restrictive. Indeed, in Section 3.3.2 we will show that if the agents perform two consecutive communication in the same iteration with any consensus matrix satisfying Assumption 4, then this is equivalent to a single communication with a consensus matrix satisfying both Assumptions 4 and 5.

3.2 Algorithm Description

In this section, we start from the parallel ADMM and gradually introduce the reader to our proposed algorithm to jointly gain insights about the underlying mechanism and motivate the role of the consensus and tracking schemes.

The update (3c) for $\lambda_{k}$ in the parallel ADMM resembles a gradient step aiming at maximizing the cost function of $\mathcal{D}$. Distributed implementations of gradient-based approaches to solve optimization problems with common decision variables in the form of $\mathcal{D}$ are well known, see e.g. [2]. Typically, each agent $i$ maintain a vector $\lambda_{i,k}\in\mathbb{R}^{p}$, representing a local version (or copy) of $\lambda_{k}$, which is iteratively updated according to a consensus-based scheme to force agreement of the local copies. If, at each iteration, the quantity $d_{k+1}$ were available to all agents, then we could propose the following local update step of $\lambda_{i,k}$ for agent $i$

[TABLE]

for all $i=1,\dots,N$.

However, update (4) cannot be implemented in a fully distributed scheme since $d_{k+1}$ is not locally available and should be computed by a central unit (cf. (3b)). In the same spirit of $\lambda_{i,k}$, a distributed counterpart for (3b) can be obtained by equipping each agent $i$ with a local auxiliary quantity $d_{i,k}\in\mathbb{R}^{p}$, which serves as a local estimate of $d_{k}$. Since $d_{k}$ is the average of $A_{i}x_{i,k}-b_{i}$, $i=1,\dots,N$ (with $\sum_{i=1}^{N}b_{i}=b$), we propose to update $d_{i,k}$ according to a (distributed) dynamic average consensus mechanism, see [19, 20]:

[TABLE]

initialized with $d_{i,0}=A_{i}x_{i,0}-b_{i}$, for all $i=1,\dots,N$.

In this way, $d_{i,k}$ acts as a distributed tracker of the (time-varying) signal $(1/N)\sum_{i=1}^{N}\left(A_{i}x_{i,k}-b_{i}\right)$. Using $d_{i,k+1}$ in place of $d_{k+1}$ in (4) makes the update of $\lambda_{i,k}$ fully distributed.

Clearly, since we propose to replace the centralized quantities $\lambda_{k}$ and $d_{k}$ with their local counterparts, the local minimization in (3a) has to be adjusted accordingly. The Tracking-ADMM is formally summarized in Algorithm 1 from the perspective of agent $i$. Specifically, the adapted local minimization to compute $x_{i,k+1}$ is shown in Step 8, where $\lambda_{k}$ and $d_{k}$ in the original centralized update (3a), are replaced by the local averages $\ell_{i,k}$ and $\delta_{i,k}$, respectively (cf. Steps 6 and 7).

Some remarks are in order. First, we shall stress that all update steps (cf. Steps 8-10) are fully distributed, as they only use quantities that are either locally known to agent $i$ or collected by agent $i$ via neighboring communications (cf. Steps 6 and 7). Moreover, the minimization in Step 8 is always well defined in view of Assumption 1.

While the initialization of $x_{i,k}$ and $\lambda_{i,k}$ can be arbitrary, the correct initialization of $d_{i,k}$ as per Step 3 is crucial, consistently with other tracking-based approaches as, e.g., the ones mentioned in the introduction (cf. [21, 22, 23, 24, 25, 26]). Sensible values for initializing Algorithm 1 are given by $x_{i,0}\in\operatorname*{argmin}_{x_{i}\in X_{i}}f_{i}(x_{i})$, $d_{i,0}=A_{i}x_{i,0}-b/N$, and $\lambda_{i,0}=0$.

Finally, the parameter $c>0$ in Step 8 and 10 is constant and is similar to the step-size of a gradient-based method, but, differently from gradient-like approaches, its value can be arbitrary.

3.3 Algorithm Analysis

In this subsection we analyze the proposed Tracking-ADMM algorithm and state its convergence properties. To ease the exposition, we state and discuss the main steps of the convergence analysis, deferring their proofs to A.

3.3.1 Aggregate Reformulation of Tracking-ADMM

First of all we shall rewrite the execution of Algorithm 1 by all agents in a compact form.

Recalling that $\boldsymbol{x}=[x_{1}^{\top}\,\cdots\,x_{N}^{\top}]^{\top}$, let us consistently denote with bold symbols the vectors collecting the corresponding quantity of all agents, i.e.,

[TABLE]

Then, the evolution of Algorithm 1 over the whole multi-agent network can be compactly written as

[TABLE]

where $f(\boldsymbol{x})=\sum_{i=1}^{N}f_{i}(x_{i})$, $A_{d}=\operatorname{blkdiag}(A_{1},\dots,A_{N})$, $X=X_{1}\times\cdots\times X_{N}$, and $W=\mathcal{W}\otimes I_{p}$. Note that, in (6), we used $\boldsymbol{\delta}_{k}=W\boldsymbol{d}_{k}$ and $\boldsymbol{\ell}_{k}=W\boldsymbol{\lambda}_{k}$, which represent the network-wide formulations of Steps 6 and 7 of Algorithm 1, respectively.

As can be seen from (6), Tracking-ADMM describes a dynamical system composed of two parts: a nonlinear dynamics given by the optimization step (cf. (LABEL:eq:all_primal_update)) and a linear dynamics given by a twofold consensus-based step (cf. (6b) and (6c)). In the forthcoming analysis, we will investigate the structural properties of these two parts separately. Then we show how to combine them to prove convergence of the proposed algorithm to an optimal solution of the original problem $\mathcal{P}$.

3.3.2 Discussion on Assumption 5

Before starting with the convergence analysis, we now show that Assumption 5 is not restrictive.

Suppose to substitute Steps 6 and 7 with

[TABLE]

and then perform the updates

[TABLE]

before Step 8. These two steps can then be compactly written as

[TABLE]

which are equivalent to a single update round with $\mathcal{W}^{2}$ in place of $\mathcal{W}$. Since, under Assumption 4, $\mathcal{W}$ is symmetric and doubly stochastic, then $\mathcal{W}^{2}$ is symmetric, doubly stochastic, and positive semidefinite, thus satisfying Assumption 5. Letting $\mathcal{E}^{2}$ be the set of edges $(i,j)$ such that there exists a $\kappa\in\mathcal{V}$ for which $(i,\kappa),(\kappa,j)\in\mathcal{E}$, it is easy to show that the graph $\mathcal{G}^{2}=(\mathcal{V},\mathcal{E}^{2})$ satisfies Assumption 3 and $\mathcal{W}^{2}$ satisfies Assumption 4 with $\mathcal{E}^{2}$ in place of $\mathcal{E}$.

Therefore, Assumption 5 can be easily satisfied by simply performing two update rounds rather than one. We choose not to alter Algorithm 1 to keep the notation light and avoid writing $\mathcal{W}^{2}$ in place of $\mathcal{W}$ everywhere.

3.3.3 Network Average Quantities

We start the analysis of Tracking-ADMM by stating some important properties regarding the network averages of the local variables $d_{i,k}$ and $\lambda_{i,k}$, which are defined as

[TABLE]

respectively. The following lemmas highlight key features of $\bar{d}_{k}$ and $\bar{\lambda}_{k}$.

Lemma 1** (Tracking property).**

Under Assumption 4 it holds that

[TABLE]

for all $k\geq 0$. ∎

Lemma 2** (Average dual update).**

Under Assumption 4, it holds that

[TABLE]

for all $k\geq 0$. ∎

As can be seen by a formal comparison, (8) and (9) mimic the update step (3b) and (3c) of the parallel ADMM. Therefore, by virtue of Lemmas 1 and 2, running Algorithm 1 results in the agents network behaving, on average, like the parallel ADMM.

3.3.4 Consensus Errors as Linear Dynamical Systems

The introduction of $\bar{d}_{k}$ and $\bar{\lambda}_{k}$ is also useful to analyze the behavior of $d_{i,k}$ and $\lambda_{i,k}$. Specifically, as usually done in consensus-based optimization algorithms, we shall study the consensus error, i.e., distance between the agents’ local quantities $d_{i,k}$ and $\lambda_{i,k}$ and their respective network average $\bar{d}_{k}$ and $\bar{\lambda}_{k}$.

To ease the notation, let us introduce the following vectors containing $N$ copies of the respective average quantities

[TABLE]

which are useful when comparing the aggregate vectors $\boldsymbol{d}_{k}$ and $\boldsymbol{\lambda}_{k}$ with their averages $\bar{d}_{k}$ and $\bar{\lambda}_{k}$. Then, define the consensus errors as

[TABLE]

and the auxiliary sequence

[TABLE]

From the expression of $\boldsymbol{d}_{k+1}$ in (6b), we can write

[TABLE]

where in (a) we used the definition of $\boldsymbol{z}_{k+1}$ and $\boldsymbol{z}_{k}$ (cf. (10)) and in (b) we exploited the following

[TABLE]

which holds due to the mixed-product property of the Kronecker product in (c) and the doubly stochasticity of $\mathcal{W}$ under Assumption 4 in (d).

Finally, letting $\mathcal{O}$ be the $N\times N$ matrix with all entries equal to $1/N$ and noticing that

[TABLE]

we can subtract the quantity $(\mathcal{O}\otimes I_{p})\boldsymbol{e}^{\boldsymbol{d}}_{k}=0$ from the right hand side of (11) to get

[TABLE]

where we set $\tilde{W}=W-(\mathcal{O}\otimes I_{p})$.

With an analogous reasoning, we can also show that

[TABLE]

where, in (a), we used (6c) and (9) pre-multiplied by $\mathbbm{1}_{N}\otimes$.

Under Assumptions 3 and 4, we have that $\|\mathcal{W}-\mathcal{O}\|<1$, then

[TABLE]

Thus, all eigenvalues of $\tilde{W}$ lie within the unit circle. This implies that the dynamical systems in (12) and (13), describing the evolution of the consensus errors, are asymptotically stable. The importance of this fact is twofold: i) (12) and (13) enjoys the input-to-state stability property: if the sequence ${\{\boldsymbol{z}_{k+1}-\boldsymbol{z}_{k}\}}_{k\geq 0}$ is bounded then also ${\{\boldsymbol{e}^{\boldsymbol{d}}_{k}\}}_{k\geq 0}$ and ${\{\boldsymbol{e}^{\boldsymbol{\lambda}}_{k}\}}_{k\geq 0}$ are bounded; ii) if $\boldsymbol{z}_{k+1}-\boldsymbol{z}_{k}$ vanishes, then also $\boldsymbol{e}^{\boldsymbol{d}}_{k}$ and $\boldsymbol{e}^{\boldsymbol{\lambda}}_{k}$ do.

Point i) is formalized in the following lemma.

Lemma 3** (Bounded sequences).**

Under Assumptions 1, 3 and 4 we have that

1. (i)

the sequences ${\{\boldsymbol{x}_{k}\}}_{k\geq 0}$ and ${\{\boldsymbol{z}_{k}\}}_{k\geq 0}$ are bounded; 2. (ii)

the sequences ${\{\boldsymbol{e}^{\boldsymbol{d}}_{k}\}}_{k\geq 0}$ and ${\{\boldsymbol{e}^{\boldsymbol{\lambda}}_{k}\}}_{k\geq 0}$ are bounded. ∎

The results in Lemmas 1, 2, and 3 can be interpreted as follows. The sequences ${\{\boldsymbol{d}_{k}\}}_{k\geq 0}$ and ${\{\boldsymbol{\lambda}_{k}\}}_{k\geq 0}$ generated by (6b) and (6c) evolve in a (possibly large) neighborhood of their averages $\bar{d}_{k}$ and $\bar{\lambda}_{k}$, which, in turn, evolve according to a parallel ADMM iteration.

3.3.5 Optimality Characterization

We now focus our attention to the analysis of the (nonlinear) update step in (LABEL:eq:all_primal_update). The following result represents the distributed counterpart of the key inequality to prove convergence of the parallel ADMM.

Proposition 1** (Local optimality).**

Under Assumptions 1 and 2 we have that

[TABLE]

for any optimal solution pair $(\boldsymbol{x}^{\star},\lambda^{\star})$ for $\mathcal{P}$ and $\mathcal{D}$, where we set $\boldsymbol{\lambda}^{\star}=\mathbbm{1}_{N}\otimes\lambda^{\star}$. ∎

If we ignore the cross term $2c[\boldsymbol{z}_{k+1}-A_{d}\boldsymbol{x}^{\star}]^{\top}\boldsymbol{e}^{\boldsymbol{\lambda}}_{k+1}$, then (14) tells us that as long as $\boldsymbol{\bar{\lambda}}_{k}$ changes from one iteration to the next one, the network of agents is, on average, getting closer to the optimal solution of $\mathcal{D}$. The additional cross term depends on the consensus error $\boldsymbol{e}^{\boldsymbol{\lambda}}_{k+1}=\boldsymbol{\lambda}_{k+1}-\boldsymbol{\bar{\lambda}}_{k+1}$ and is clearly a distinctive feature of Tracking-ADMM (as opposed to the parallel ADMM) arising from its distributed nature. In the parallel ADMM, a condition similar to (14) can be derived and leveraged to prove the algorithm convergence. However, as expected, this is not the case here: (14) is not sufficient to prove convergence of the proposed scheme and we are required to study the interplay between the minimization step in (LABEL:eq:all_primal_update) and the two consensus steps in (6b) and (6c).

3.3.6 Tracking-ADMM Convergence

The convergence of Algorithm 1 is divided into two results. The first theorem proves that the two consensus mechanism are successful, i.e., agents eventually agree on a common value for the tracker and for the dual variables (cf. points i and ii), and the algorithm converges to a solution that is feasible for $\mathcal{P}$ (cf. point iii); the second theorem proves that the primal and dual solution obtained by Algorithm 1 are optimal for $\mathcal{P}$ and $\mathcal{D}$, respectively.

Theorem 1** (Convergence).**

Under Assumptions 1-5, the sequences generated by Tracking-ADMM satisfy:

1. (i)

$\lim_{k\to\infty}\|\boldsymbol{e}^{\boldsymbol{d}}_{k}\|=0$, 2. (ii)

$\lim_{k\to\infty}\|\boldsymbol{e}^{\boldsymbol{\lambda}}_{k}\|=0$, 3. (iii)

$\lim_{k\to\infty}\|\bar{d}_{k}\|=0$, 4. (iv)

${\{\|\boldsymbol{\bar{\lambda}}_{k}-\boldsymbol{\lambda}^{\star}\|^{2}+c^{2}\|\boldsymbol{z}_{k}-A_{d}\boldsymbol{x}^{\star}\|^{2}\}}_{k\geq 0}$* is convergent,*

for any optimal solution pair $(\boldsymbol{x}^{\star},\lambda^{\star})$ for $\mathcal{P}$ and $\mathcal{D}$, with $\boldsymbol{\lambda}^{\star}=\mathbbm{1}_{N}\otimes\lambda^{\star}$. ∎

The main idea behind the proof of Theorem 1 is the following. We start by considering (12) and (13) together as a single dynamical system for the consensus updates (the consensus system). Since the consensus system is asymptotically stable, then we can build a positive-definite quadratic Lyapunov function which monotonically decreases along its free evolution. However, the consensus system is not autonomous, thus the inequality describing the variation of such Lyapunov function across iterations contains terms that are not defined in sign, so that it is no longer necessarily decreasing. The idea is to properly combine such inequality with (14) of Proposition 1 so as to “balance out” the terms that are not defined in sign. In this way, we obtain an aggregate descent condition that allows us to prove the asymptotic stability of the overall (nonlinear) dynamical system modeling our Tracking-ADMM algorithm.

The reader should note that point i of Theorem 1 implies that, for all $i=1,\dots,N$, $d_{i,k}$ tends to $\bar{d}_{k}$ as $k\to\infty$. In light of (8), this ensures that each agent is able to locally assess the amount of infeasibility of the current primal iterates for the coupling constraint. This feature is of great importance in those applications where feasibility up to a given tolerance is sufficient.

Finally, the following result shows that by running Algorithm 1 agents are able to compute optimal solutions to $\mathcal{P}$ and $\mathcal{D}$ in a fully distributed way.

Theorem 2** (Optimality).**

Under Assumptions 1-5, the sequences generated by Tracking-ADMM are such that

1. (i)

any limit point of the primal sequence ${\{\boldsymbol{x}_{k}\}}_{k\geq 0}$ is an optimal solution $\boldsymbol{x}^{\star}$ of $\mathcal{P}$; 2. (ii)

each dual sequence ${\{\lambda_{i,k}\}}_{k\geq 0}$, $i=1,\dots,N$, converges to the same optimal solution $\lambda^{\star}$ of $\mathcal{D}$. ∎

We point out that Theorem 2 guarantees that the sequences ${\{x_{i,k}\}}_{k\geq 0}$ are asymptotically optimal, hence feasible for both local and the coupling constraints. Notice also that Algorithm 1 works with (not necessarily strictly) convex primal cost functions and does not require any primal recovery procedure.

4 Numerical Study

To demonstrate the effectiveness of our approach, we test Tracking-ADMM on a modified version of the Plug-in Electric Vehicles Optimal Charging Schedule problem described in [37]. A fleet of $N$ electric vehicles has to select an overnight charging schedule so as to minimize the sum of the electricity costs for charging the internal battery of each vehicle while complying with vehicle-level constraints (user requirements on the final state of charge and battery physical limitations) and a network-wide constraint (maximum power that the grid can deliver). For simplicity we focus on the “only charging” case of the charging problem in [37]. Differently from [37], where the vehicles can decide only whether to charge or not their internal battery at a fixed charging rate, here we allow each vehicle to optimize also the charging rate. The resulting optimization program can then be formalized as follows

[TABLE]

where $X_{i}$ are bounded polyhedral sets. We refer the reader to [37] for the precise formulations of the quantities in (15).

The coupling constraints in (15) is an inequality, which makes (15) not fitting the formulation in $\mathcal{P}$. However, we can easily turn the coupling constraint into an equality by introducing $p$ additional slack variables for each vehicle. Problem (15) thus becomes

[TABLE]

which now fits $\mathcal{P}$ with the identifications $x_{i}=[\xi_{i}^{\top}\;s_{i}^{\top}]^{\top}$, $f_{i}(x_{i})=\gamma_{i}^{\top}\xi_{i}$, $A_{i}=[\tilde{A}_{i}\;I_{p}]$, and $X_{i}=\Xi_{i}\times S_{i}$ where $S_{i}=\{s_{i}:0\leq s_{i}\leq\bar{s}\}$, for some $\bar{s}\geq b$, for all $i=1,\dots,N$.

In our simulation we considered $N=100$ vehicles. Each vehicle $i$ has 24 decision variables $\xi_{i}$ plus 24 additional slack variables $s_{i}$ ($n_{i}=48$) and the local constraint set $X_{i}$ is defined by 245 inequalities. The number of coupling constraint, and consequently the number of Lagrange multipliers, is $p=24$. The communication network and the weight matrix $\mathcal{W}$ have been generated at random so as to satisfy Assumptions 3, 4, and 5.

We run our Tracking-ADMM for 200 iterations, initialized as suggested in Section 3.2, setting $c=10^{-4}$.

In Figure 1 we report the behavior across iterations of the difference $\sum_{i=1}^{N}f_{i}(x_{i,k})-f^{\star}$ between the value of the cost function achieved by the primal tentative solution $x_{i,k}$ and the optimal cost $f^{\star}$ computed by a centralized solver (upper plot), and the infinite norm $\|\sum_{i=1}^{N}A_{i}x_{i,k}-b\|_{\infty}$ of the joint constraints violation (lower plot). As it can be seen from the picture, the proposed algorithm reaches feasibility and optimality. Note that, differently from the approaches based on dual decomposition, the primal iterates generated by Algorithm 1 achieve feasibility of the coupling constraint without requiring any recovery procedure (which is known to slow down the algorithm convergence).

In Figure 2 we also plot the evolution of the agents local estimates ${\{\lambda_{i,k}\}}_{k\geq 0}$ of the optimal Lagrange multipliers $\lambda^{\star}$ (red triangles). From the picture the reader can easily see that local estimates reach consensus and optimality.

5 Conclusions

In this paper we have proposed a novel distributed method to solve constraint-coupled convex optimization problems in which the sum of local cost functions needs to be minimized while satisfying both individual constraints (involving one component of the decision variables) and a common linear coupling constraint (involving all the components). The distributed algorithm combines the (parallel) ADMM algorithm tailored for this class of optimization problems with a dynamic tracking mechanism. We proved that each agent asymptotically computes an optimal dual solution and its portion of an optimal solution to the target (primal) problem. Moreover, the tracking scheme allows agents to obtain a local measure of the coupling-constraint violation. Numerical computations corroborated the theoretical results.

Appendix A Proofs

Proof of Lemma 1

We prove (8) by induction. For $k=0$, given the initialization in Step 3, we have that

[TABLE]

which proves that (8) holds for $k=0$. Assume now that (8) holds up to $k$. If we can prove that (8) holds also for $k+1$, then the proof is completed. Indeed,

[TABLE]

where in (a) we used Steps 9 and 6 in (b) we exchanged the summations and we used column stochasticity of the weights (cf. Assumption 4), in (c) we used the definition of $\bar{d}_{k}$ and the induction step in (d). ∎

Proof of Lemma 2

By definition of $\bar{\lambda}_{k+1}$, we can write

[TABLE]

which can be further elaborated as follows

[TABLE]

where in (a) we used the definition of $\lambda_{i,k+1}$ in Step 10, in (b) we switched the summations, and in (c) we used column stochasticity of the weights (cf. Assumption 4) together with the definition of $\bar{\lambda}_{k}$ and $\bar{d}_{k+1}$. ∎

Proof of Lemma 3

By (LABEL:eq:all_primal_update) we have that $\boldsymbol{x}_{k}\in X$ for all $k$ and, under Assumption 1, $X$ is bounded. Therefore the sequence ${\{\boldsymbol{x}_{k}\}}_{k\geq 0}$ is bounded. Under Assumption 4, from (8) in Lemma 1 we have that $\bar{d}_{k}=(1/N)(A\boldsymbol{x}_{k}-b)$ and, therefore, also ${\{\boldsymbol{\bar{d}}_{k}\}}_{k\geq 0}$ is bounded since $\boldsymbol{\bar{d}}_{k}=\mathbbm{1}_{N}\otimes\bar{d}_{k}$. Moreover, since $\boldsymbol{z}_{k}=A_{d}\boldsymbol{x}_{k}-\boldsymbol{\bar{d}}_{k}$, then also ${\{\boldsymbol{z}_{k}\}}_{k\geq 0}$ is bounded. This proves point i.

Being the sequence ${\{\boldsymbol{z}_{k}\}}_{k\geq 0}$ bounded, we have that ${\{\boldsymbol{z}_{k+1}-\boldsymbol{z}_{k}\}}_{k\geq 0}$ is also bounded. Therefore, exploiting the input-to-state stability property of the dynamical system in (12) under Assumptions 3 and 4, we can conclude that ${\{\boldsymbol{e}^{\boldsymbol{d}}_{k}\}}_{k\geq 0}$ is a bounded sequence. Finally, owing to the input-to-state stability of the dynamical system in (13) under Assumptions 3 and 4 together with boundedness of ${\{\boldsymbol{e}^{\boldsymbol{d}}_{k}\}}_{k\geq 0}$, we also have that ${\{\boldsymbol{e}^{\boldsymbol{\lambda}}_{k}\}}_{k\geq 0}$ is bounded, thus proving point ii and concluding the proof. ∎

Proof of Proposition 1

By [17, Lemma 4.1] applied to (LABEL:eq:all_primal_update), we have that $\boldsymbol{x}_{k+1}$ satisfies

[TABLE]

for any $\boldsymbol{x}\in X$. Recalling that $\boldsymbol{d}_{k+1}=W\boldsymbol{d}_{k}+A_{d}\boldsymbol{x}_{k+1}-A_{d}\boldsymbol{x}_{k}$ and $\boldsymbol{\lambda}_{k+1}=W\boldsymbol{\lambda}_{k}+c\,\boldsymbol{d}_{k+1}$ (cf. (6b) and (6c), respectively), we have

[TABLE]

for any $\boldsymbol{x}\in X$. Under Assumption 2, we can set $\boldsymbol{x}$ equal to some optimal solution $\boldsymbol{x}^{\star}$ of $\mathcal{P}$, to obtain

[TABLE]

Under Assumption 2, by the Saddle Point Theorem in [17, pag. 665] we have that

[TABLE]

for all $\boldsymbol{x}\in X$ and for all $\lambda$. Given (A.17) and the fact that $L(\boldsymbol{x}^{\star},\lambda^{\star})=f(\boldsymbol{x}^{\star})$, we have that

[TABLE]

for all $\boldsymbol{x}\in X$ and any $\lambda^{\star}$. Writing $b=A\boldsymbol{x}^{\star}$ and noticing that ${\lambda^{\star}}^{\top}[A\boldsymbol{x}-A\boldsymbol{x}^{\star}]={\boldsymbol{\lambda}^{\star}}^{\top}[A_{d}\boldsymbol{x}-A_{d}\boldsymbol{x}^{\star}]$, it holds

[TABLE]

Setting $\boldsymbol{x}=\boldsymbol{x}_{k+1}$, simplifying $f(\boldsymbol{x}_{k+1})$, and bringing everything on the left-hand side, we have

[TABLE]

Adding and subtracting $\boldsymbol{\bar{\lambda}}_{k+1}^{\top}[A_{d}\boldsymbol{x}_{k+1}-A_{d}\boldsymbol{x}^{\star}]$ yields

[TABLE]

Focusing on the first term in (A.18), we have

[TABLE]

where (a) is due to the block structure of the vectors involved in the product, (b) is given by the fact that $[\bar{\lambda}_{k+1}-\lambda^{\star}]$ does not depend on $i$, (c) is given by the definition of $A$, (d) is given by the fact that $A\boldsymbol{x}^{\star}=b$, (e) is due to (8) in Lemma 1, (f) is given by (9) in Lemma 2, and (g) uses the definitions of $\boldsymbol{\bar{\lambda}}_{k}$ and $\boldsymbol{\lambda}^{\star}$. Owing to the fact that

[TABLE]

and considering the second term in (A.18), we have

[TABLE]

where in (a) we used the definition of $\boldsymbol{z}_{k+1}$ (cf. (10)). Using (A.19) and (A.20) in (A.18) and multiplying by $2c$ we get

[TABLE]

Finally, using

[TABLE]

with the identifications $a=\boldsymbol{\bar{\lambda}}_{k+1}$, $b=\boldsymbol{\lambda}^{\star}$ and $c=\boldsymbol{\bar{\lambda}}_{k}$, (14) follows, thus concluding the proof. ∎

Proof of Theorem 1

The dynamical evolution of $\boldsymbol{e}^{\boldsymbol{d}}_{k}$ and $\boldsymbol{e}^{\boldsymbol{\lambda}}_{k}$, given in equations (12) and (13), respectively, can be arranged as

[TABLE]

with

[TABLE]

and $\boldsymbol{u}_{k}=c[\boldsymbol{z}_{k}-A_{d}\boldsymbol{x}^{\star}]$. Let $P=P^{\top}\succ 0$ and consider

[TABLE]

where in (a) we used the dynamics (A.22) and added $\pm\boldsymbol{e}_{k}$ inside the norm, in (b) we expanded the square, while in (c) we performed some algebraic manipulations. Summing (A.23) to the optimality-based inequality (14) (cf. Proposition 1) yields

[TABLE]

with $H=[I\quad 0]$. Notice that, choosing $P$ such that $P-F^{\top}PF\succ 0$, if the terms at $k+1$ were analogous to the terms at $k$, then (A.24) would describe a descent condition. To enforce this, we need to find a $P$ satisfying the following conditions

[TABLE]

Satisfying (A.25) translates into (A.24) being a non-expansive condition for the sequence ${\{\|\boldsymbol{\bar{\lambda}}_{k+1}-\boldsymbol{\lambda}^{\star}\|^{2}+2\boldsymbol{u}_{k+1}^{\top}H\boldsymbol{e}_{k+1}+\|G\boldsymbol{u}_{k+1}-\boldsymbol{e}_{k+1}\|_{P}^{2}\}}_{k\geq 0}$, i.e., a stability condition for the nonlinear dynamical system representing Algorithm 1.

Take $P$ partitioned in blocks as follows

[TABLE]

Owing to the fact that all eigenvalues of $\tilde{W}$ lies in the open unit circle, thus $(I-\tilde{W})$ is invertible and, from the equality $MM^{-1}=M^{-1}M=I$ with $M=(I-\tilde{W})$, we have the following identities

[TABLE]

Noticing that

[TABLE]

condition (A.25a) can be rewritten as $G^{\top}P=H(I-F)^{-1}$ and translates into the following constraints on the blocks of $P$:

[TABLE]

where in (a) we used (A.26a), while (A.27c) follows by simply summing (A.27a) and (A.27b) and will be useful in the forthcoming discussion. Note that, from (A.27a), we can further deduce that $P_{2}=P_{2}^{\top}$ since both $P_{1}$ and $(I-\tilde{W})^{-1}$ are symmetric matrices owing to $P$ being symmetric and Assumption 4, respectively.

We can thus express the matrix $P-F^{\top}PF$ as

[TABLE]

where in (a) we used (A.27) and in (b) we leveraged identities in (A.26) together with (A.27a) and (A.27b) to express $P_{2}$ and $P_{3}$ as a function of $P_{1}$ only.

Next, we find a value for $P_{1}\succ 0$ ensuring that $P-F^{\top}PF$ is positive definite. This requirement can be posed in terms of its blocks by means of the Schur complement, i.e., $P-F^{\top}PF\succ 0$ if and only if

[TABLE]

Take $P_{1}=2I$ and let $T\Lambda T^{\top}=\tilde{W}$ be the eigendecomposition of the symmetric matrix $\tilde{W}$, with $TT^{\top}=I$ and $\Lambda$ being a diagonal matrix containing the eigenvalues of $\tilde{W}$. Condition (A.29a) is trivially satisfied by this choice of $P_{1}$, while (A.29b) becomes

[TABLE]

which is equivalent to $I+2\Lambda-3\Lambda^{2}\succ 0$ and always holds true when all the eigenvalues of $\tilde{W}$ lies within $(-1/3,1)$. In view of Assumption 5, this condition is guaranteed.

Finally, using again Schur’s complement, it is also easy to prove that $P_{1}=2I$ satisfies condition (A.25c). In fact, we have that

[TABLE]

if and only if

[TABLE]

Condition (A.30b) is satisfied since $\tilde{W}$ has all eigenvalues in the open unit circle and hence $I-\tilde{W}\succ 0$.

From (A.25a) we can rewrite (A.24) as

[TABLE]

where we set $Q=P-F^{\top}PF\succ 0$. Rearranging some terms, and summing (A.31) from $k=0$ to $M-1$, with $M\in\mathbb{N}$, we have

[TABLE]

where in the second inequality we neglected the terms $-\|\boldsymbol{\bar{\boldsymbol{\lambda}}}_{M}-\boldsymbol{\lambda}^{\star}\|^{2}$ and $-\|G\boldsymbol{u}_{M}-\boldsymbol{e}_{M}\|_{P}^{2}$ in the right hand side since they are non-positive (recall $P\succ 0$) and we used the fact that $|2\boldsymbol{u}_{k}^{\top}H\boldsymbol{e}_{k}|\leq\mathcal{C}$ for all $k\geq 0$ owing to the results of Lemma 3. Taking the limit as $M\to\infty$ finally yields

[TABLE]

which, recalling that $Q\succ 0$, implies that $\lim_{k\to\infty}\|\boldsymbol{\bar{\lambda}}_{k+1}-\boldsymbol{\bar{\lambda}}_{k}\|=0$ and $\lim_{k\to\infty}\|\boldsymbol{e}_{k}\|=0$. Hence, $\lim_{k\to\infty}\|\boldsymbol{e}^{\boldsymbol{d}}_{k}\|=0$ and $\lim_{k\to\infty}\|\boldsymbol{e}^{\boldsymbol{\lambda}}_{k}\|=0$, thus proving points i and ii, respectively.

Since $\boldsymbol{\bar{d}}_{k+1}=\frac{1}{c}(\boldsymbol{\bar{\lambda}}_{k+1}-\boldsymbol{\bar{\lambda}}_{k})$, we also have that $\lim_{k\to\infty}\|\boldsymbol{\bar{d}}_{k}\|=0$ and therefore $\lim_{k\to\infty}\|\bar{d}_{k}\|=0$, thus proving point iii.

From (A.31) we also have that the sequence

[TABLE]

is non-increasing, bounded below since $|2\boldsymbol{u}_{k}^{\top}H\boldsymbol{e}_{k}|\leq\mathcal{C}$ and $\|G\boldsymbol{u}_{k}-\boldsymbol{e}_{k}\|_{P}^{2}\geq 0$ (recall $P\succ 0$), and, therefore, convergent. Since $\|\boldsymbol{e}_{k}\|$ is vanishing and ${\{\boldsymbol{u}_{k}\}}_{k\geq 0}$ is bounded, we have that also the sequence ${\{\|\boldsymbol{\bar{\lambda}}_{k}-\boldsymbol{\lambda}^{\star}\|^{2}+\|G\boldsymbol{u}_{k}\|_{P}^{2}\}}_{k\geq 0}$ is convergent. A straightforward computation using (A.27c) shows that $G^{\top}PG=(I-\tilde{W})^{-2}\succ 0$ which, recalling the definition of $\boldsymbol{u}_{k}=c\,(\boldsymbol{z}_{k}-A_{d}\boldsymbol{x}^{\star})$, implies that the sequence

[TABLE]

is convergent, which finally implies that

[TABLE]

is convergent due to norm equivalence, thus proving point iv and concluding the proof. ∎

Proof of Theorem 2

Following the first steps in the proof of Proposition 1, for any optimal solution $\boldsymbol{x}^{\star}$ for $\mathcal{P}$, we can write

[TABLE]

By adding and subtracting $\boldsymbol{\bar{\lambda}}_{k+1}^{\top}[A_{d}\boldsymbol{x}_{k+1}-A_{d}\boldsymbol{x}^{\star}]$ to the left-hand side, we have

[TABLE]

Noticing that

[TABLE]

where the last inequality is due to Lemma 1, we have that

[TABLE]

From Theorem 1iv we have that the sequence

[TABLE]

is bounded. Since ${\{\boldsymbol{z}_{k}\}}_{k\geq 0}$ is also bounded by Lemma 3, we have that ${\{\boldsymbol{\bar{\lambda}}_{k}\}}_{k\geq 0}$ and, therefore, also ${\{\bar{\lambda}_{k}\}}_{k\geq 0}$ are both bounded. Recalling that ${\{\boldsymbol{x}_{k}\}}_{k\geq 0}$ is bounded, by Theorem 1 points iii and ii respectively, it holds

[TABLE]

and therefore

[TABLE]

Since $\boldsymbol{x}_{k}\in X$ for all $k\geq 0$ and, by Lemma 1 and Theorem 1iii,

[TABLE]

we know that any limit point $\boldsymbol{\tilde{x}}$ of the sequence ${\{\boldsymbol{x}_{k}\}}_{k\geq 0}$ is feasible for $\mathcal{P}$ and, thus, satisfies $f(\boldsymbol{\tilde{x}})\geq f(\boldsymbol{x}^{\star})$. Combining this fact with (A.35) we can conclude that $f(\boldsymbol{\tilde{x}})=f(\boldsymbol{x}^{\star})$ for any limit point $\boldsymbol{\tilde{x}}$. Since all limit points of ${\{\boldsymbol{x}_{k}\}}_{k\geq 0}$ are feasible and achieve the optimal value, they are optimal solutions for $\mathcal{P}$, proving Theorem 2i.

Next, we prove the convergence of the local dual variables estimates to an optimal dual solution. Similarly to the proof of Proposition 1, by applying [17, Lemma 4.1] to Step 8 we have that $x_{i,k+1}$ satisfies

[TABLE]

for all $x_{i}\in X_{i}$. Noticing that $d_{i,k+1}=\delta_{i,k}+A_{i}x_{i,k+1}-A_{i}x_{i,k}$ and $\lambda_{i,k+1}=\ell_{i,k}+c\,d_{i,k+1}$ from Steps 9 and 10 respectively and subtracting $\lambda_{i,k+1}^{\top}b_{i}$, we have

[TABLE]

for all $x_{i}\in X_{i}$, which is equivalent to

[TABLE]

and, thus, by definition of $\varphi_{i}(\lambda)$ in (2), we have

[TABLE]

if $b_{1},\dots,b_{N}$ are such that $\sum_{i=1}^{N}b_{i}=b$. Setting $b_{i}=A_{i}x_{i}^{\star}$, summing over $i=1,\dots,N$, adding and subtracting $\boldsymbol{\bar{\lambda}}_{k+1}^{\top}[A_{d}\boldsymbol{x}_{k+1}-A_{d}\boldsymbol{x}^{\star}]$ to the right hand side, and exploiting again (A.32), we obtain

[TABLE]

Since the sequence ${\{\bar{\lambda}_{k+1}\}}_{k\geq 0}$ is bounded, it admits a convergent subsequence $\{(\bar{\lambda}_{k+1},\boldsymbol{x}_{k+1})\}_{k\in\mathcal{K}}$ with $\mathcal{K}\subseteq\mathbb{N}$. Let $(\tilde{\lambda},\boldsymbol{\tilde{x}})$ be its limit point. Now, taking the limit of (A.36) across $\mathcal{K}$, we have

[TABLE]

where in (a) we used (A.33) and (A.34), the fact that limit points of ${\{\boldsymbol{x}_{k}\}}_{k\geq 0}$ achieve the optimal value in (b), and weak duality in (c). On the other hand, from (2), we have

[TABLE]

for all $x_{i}\in X_{i}$, where the first equality is due to $\lim_{k\to\infty}\lambda_{i,k}=\bar{\lambda}_{k}$, $i=1,\dots,N$ (recall $\lim_{k\to\infty}\|\boldsymbol{e}^{\boldsymbol{\lambda}}_{k}\|=0$), and in the last equality we plugged the limit point $\tilde{\lambda}$ of $\{\bar{\lambda}_{k+1}\}_{k\in\mathcal{K}}$. Setting $x_{i}$ in (A.38) equal to the minimizer of $\varphi_{i}(\tilde{\lambda})$ and summing over $i=1,\dots,N$ yields

[TABLE]

Combining the latter with (A.37) gives

[TABLE]

which means that $\tilde{\lambda}$ is optimal for $\mathcal{D}$.

Recall that, from Theorem 1iv, we have that the sequence

[TABLE]

is convergent for any pair $\boldsymbol{x}^{\star}$ and $\lambda^{\star}$ which are optimal for $\mathcal{P}$ and $\mathcal{D}$, respectively. We can thus take $\boldsymbol{\lambda}^{\star}=\mathbbm{1}_{N}\otimes\tilde{\lambda}$ and $\boldsymbol{x}^{\star}=\boldsymbol{\tilde{x}}$ and conclude that

[TABLE]

But since ${\{\|\boldsymbol{\bar{\lambda}}_{k}-\boldsymbol{\lambda}^{\star}\|^{2}+c^{2}\|\boldsymbol{z}_{k}-A_{d}\boldsymbol{x}^{\star}\|^{2}\}}_{k\geq 0}$ is convergent, all its limit points are the same, and therefore

[TABLE]

for some primal optimal solution $\boldsymbol{x}^{\star}$ and some optimal dual solution $\lambda^{\star}$, which implies that, for all $i=1,\dots,N$,

[TABLE]

since $\boldsymbol{e}^{\boldsymbol{\lambda}}_{k}$ is vanishing. So that Theorem 2ii follows. ∎