text:
\section{Introduction}
This paper considers  time-varying optimization problems for network systems, where objective and constraints evolve over time~\cite{popkov2005gradient,simonetto2014double,Shahrampour2018,hall2015online,bernstein2018online,simonetto2016class}. The applicability of time-varying optimization problems is evident in a number of domains including power grids~\cite{dall2018optimal,tang2017real,hauswirth2016projected}, communication systems~\cite{Low1999Flow,Chen12}, and online methods in signal processing \cite{asif2014sparse}, just to name a few; see also the representative works~\cite{rahili2017distributed,fazlyab2016self,neely2017online} and~\cite{dall2019optimization} for additional time-varying models and application examples.

In particular, assume that the temporal domain is discretized as $\{k \tau, k \in \mathbb{N}\}$, with $\tau>0$ being a given interval time~\cite{simonetto2016class,dall2019optimization}. This paper focuses on time-varying optimization problems in the following form~\cite{popkov2005gradient,simonetto2014double}:
\begin{subequations}
\label{P1}
\begin{align}
\textrm{(P1)~~} \min_{{\bf x}\in\mathbb{R}^m,{\bf y}\in\mathbb{R}^n}~&f^{(k)}({\bf x}) + g^{(k)}({\bf y})\label{P1C}\\
\text{s.t.}~&{\bf A}^{(k)}{\bf x}+{\bf B}^{(k)}{\bf y}={\bf b}^{(k)}\label{constraint_intro}\\
&{\bf x}\in\mathcal{X}^{(k)}, {\bf y}\in\mathcal{Y}^{(k)},
\end{align}
\end{subequations}
where $k$ is the index for the time steps $\{k \tau, k \in \mathbb{N}\}$;  $f^{(k)}(\cdot)$ and $g^{(k)}(\cdot)$ are strongly convex functions (for all times $k \tau$); ${\bf A}^{(k)}\in\mathbb{R}^{\ell\times m}, {\bf B}^{(k)}\in\mathbb{R}^{\ell\times n}$ are time-varying matrices; and, $\mathcal{X}^{(k)},\mathcal{Y}^{(k)}$ are convex constraint sets. For optimization and learning problems with streams of data, the interval $\tau$ coincides with the inter-arrival date of data points~\cite{dall2019optimization}; when problem \eqref{P1} is associated with a network, the functions $f^{(k)}(\cdot)$ and $g^{(k)}(\cdot)$ can capture performance objectives that evolve over time, whereas~\eqref{constraint_intro} can capture time-varying physical or logical interactions in the network. Denoting as ${\bf x}^{*,(k)}$ an optimal solution of~\eqref{P1} at time $k \tau$, the optimization model~\eqref{P1} leads to an optimal \emph{trajectory}. The problem addressed in this paper pertains to the development of algorithms that enable \emph{tracking of the optimal trajectory} $\{{\bf x}^{*, (k)}\}_{k \in \mathbb{N}}$. It is worth noticing that, for  problems with inequality constraints, one can always add slack variables to re-write them as equality constraints -- thus fitting  the formulation \eqref{P1}.

Previous efforts\cite{kar2011gossip,olfati2007consensus} have addressed dynamic (consensus-type) problems as a series of static problems, and have assumed a time-scale separation between algorithms and variability of the problem,  so that convergence is reached for each time. However, this might not be  when the system parameters and the problem inputs change fast, at a time scale that is comparable with the execution of one (or a few) algorithmic steps; \cite{rahili2017distributed,fazlyab2016self, cherukuri2017saddle} successfully approach dynamic problem in continuous time, but only for isolated systems where time-varying exogenous inputs are available at a central processor.

In recent years, an intensive research has focused on real-time implementation: \cite{rahili2017distributed} presents a control algorithm for real-time multi-agent systems with the ability to track optimal trajectory, however, only the cost function is time-varying; \cite{tang2017real} proposes an online algorithm for optimal power flow problem based on quasi-Newton method. It can be shown that proposed algorithm is able to provide suboptimal solution at a fast timescale. The tracking ability hinges on the accurate estimation of second order information. For the same application, \cite{dall2018optimal} leverages dual subgradient method and system feedback (or measurements) to design a tracking algorithm based on a double smoothing strategy. Regularization terms are added in both primal and dual subproblems to prove Q-linear convergence to a neighborhood of optimal solution for each time instance; \cite{bernstein2018online} further extend double smoothing algorithm to more general settings and provide a regret analysis. The authors in~\cite{onlineSaddle} presents a saddle-point method for networked online convex optimization, and provide a regret analysis for problems with a Lipshitz-continuous function. Lastly,~\cite{colombino2019online,zheng2019implicit} considered online optimization methods for output regulation problems in dynamical systems. It is also worth mentioning that prediction-correction methods have been utilized to solve time-varying convex optimization problems; see e.g.,~\cite{simonetto2016class,fazlyab2017prediction} and the recent survey paper~\cite{simonetto2020time}.

This paper proposes the development of an online algorithm for time-varying convex problems based on the  alternating direction method of multipliers (ADMM) method\cite{boyd2011distributed}. Using a quadratic regularization term, ADMM can allow one to deal with nonsmooth terms,  and it exhibits improved convergence properties relative to dual (sub)gradient methods, especially for problems with ill-conditioned dual functions \cite{deng2016global}\cite{hong2017linear}. The paper present a new algorithm that has following characteristics: i) at each step the primal subproblems are solved via proximal gradient descent -- providing favorable scalability to large-scale problems and accommodating non-smooth objectives; ii) a dual perturbation method is utilized, where the dual variables are suitably perturbed at every iteration to gain in convergence rate. Related works along this line include the following: \cite{ling2014decentralized} leverages ADMM to solve a real-time multi-agent problem. But it differs from the present work because  it  considers only consensus constraints (a special case of our general formulation); \cite{cao2019dynamic} considers a dynamic sharing problem, and convergence to a neighborhood is provided under standard assumptions; however, the constraint is also a special case of our formulation.

We note that relative to existing online primal-dual methods~\cite{onlineSaddle,simonetto2014double,bernstein2018online}, the proposed method can handle non-differentiable costs. We also note that the choice of ADMM as opposed to, e.g. an Arrow-Hurwicz method, is due to the following two reasons: i) problem \eqref{P1} has a two-block structure (i.e., ${\bf x},{\bf y}$) which can be effectively handled by ADMM; and, ii) the objective function has nonsmooth terms, which would nevertheless require modifications of the Arrow-Hurwicz method.

To summarize, this paper has the following main contributions:
\begin{itemize}
    \item We develop an online proximal-ADMM algorithm for solving time-varying optimization problems of  form~\eqref{P1};
    \item We provide  convergence analysis, which shows that the proposed algorithm can track the optimal solution trajectory of~\eqref{P1} under mild assumptions; in particular, our methodology does not require smoothness of the cost and does not rely on the full-rankness of the constraint matrix  (see Table \ref{tab:conditions} for detailed comparison between the convergence conditions of a few algorithms).
\end{itemize}
\begin{table}
\begin{center}
\begin{tabular}{ |p{1.0cm}||p{1.3cm}|p{1.5cm}|p{1.7cm}| p{1.1cm}| }
 \hline
 & {\bf Strong Convexity} & {\bf Lipschitz Continuity} & {\bf Full Row Rank} & {\bf Optimality} \\
 \hline
 \hline
 \multirow{4}{4em}{\bf Classic ADMM} & $f^{(k)}$ & $\nabla f^{(k)}$&${\bf A}^{(k)},({\bf B}^{(k)})^T$ & \multirow{4}{4em}{Optimal solution} \\
 \cline{2-4}&$f^{(k)}, g^{(k)}$ & $ \nabla f^{(k)} $ & ${\bf A}^{(k)}$ &\\
 \cline{2-4}&$f^{(k)}$ & $\nabla f^{(k)}, \nabla g^{(k)}$&$({\bf B}^{(k)})^T$ &  \\
 \cline{2-4}&$f^{(k)}, g^{(k)}$ & $\nabla f^{(k)}, \nabla g^{(k)}$&\xmark &\\
 \hline
 {\bf Proposed Algorithm} &   $f^{(k)}, g^{(k)}$  & \xmark   & \xmark & Perturbed solution [cf.~\eqref{AKKT}] \\
 \hline
\end{tabular}
\end{center}
\caption{Trade off between optimality and conditions for linear convergence \cite{deng2016global}. }\label{tab:conditions}
\end{table}

Our previous work \cite{zhang2017dynamic} focuses on ADMM-based online algorithms  to track a solution of a domain-specific linearized AC optimal power flow (OPF) problem in power grids. In this application domain, this work significantly extends \cite{zhang2017dynamic} in the following ways:

\begin{itemize}
    \item We consider more general linearized OPF formulations, which can be used to
    deal with, for example, OPF problems in a distributed setting, where the power system is divided into areas~\cite{autonomous_grids}; and,
    \item A different algorithm which works under  milder conditions and has wider applicability is proposed.
\end{itemize}

The remainder of paper is organized as follows. Section \RomanNumeralCaps 2 will give the general time-varying problem formulation. Section \RomanNumeralCaps 3 will introduce our online algorithm. Section \RomanNumeralCaps 4 will apply proposed algorithm to two applications, one is in power systems, the other one is route selection. Tracking ability is shown in \RomanNumeralCaps 5 and \RomanNumeralCaps 6 in the form of convergence analysis and simulation, respectively.

\section{Problem Formulation}

Consider the time-varying  problem \eqref{P1}\footnote{Throughout this paper, boldface characters denote vectors or matrices; characters with superscript ${(k)}$ denote time varying iterates and parameters; for a given vector ${\bf x}$ and matrix ${\bf G}$, $\|{\bf x}\|^2_{{\bf G}} := {\bf x}^T{\bf G}{\bf x}$; $<{\bf x}, {\bf y}>$ denotes the inner product between the vectors ${\bf x}$ and ${\bf y}$. Given a non-differentiable function $h$, the proximal operator is defined as
$\text{prox}_{h}({\bf x}) = \arg\min_{\bf z}\|{\bf z}-{\bf x}\|^2 + h({\bf z}).$}.  At time $k$, if problem \eqref{P1} is solved to global optimality, then we say that the {\it perfect tracking} is achieved. However, in many applications~\cite{dall2019optimization} such perfect tracking may not be possible because before the problem at time $k$ is solved, it may have already evolved to a new problem. Specifically, iterative algorithms often involve multiple iterations of computing and communication, and by the time algorithms converge for time $k$, problem parameters such as  ${\bf A}^{(k)},{\bf B}^{(k)}, {\bf b}^{(k)}$ might have already changed. Therefore it is desirable to design algorithms with certain ``tracking ability", which means that the iterates can be continuously steered to stay close to the time-varying optimal solutions.

Let us reformulate problem \eqref{P1} as follows. First, we rewrite the time-varying constraint sets $\mathcal{X}^{(k)},\mathcal{Y}^{(k)}$ into indicator functions in the objective; and then we separate objective into non-differential functions $f_0^{(k)}({\bf x}), g_0^{(k)}({\bf y})$ and differential functions $f_1^{(k)}({\bf x}), g_1^{(k)}({\bf y})$. At time $k$ we consider the following  time-varying problem:
\begin{subequations}
\label{P2}
\begin{align}
\textrm{(P2)~~}\min_{{\bf x}\in\mathbb{R}^m,{\bf y}\in\mathbb{R}^n}~&f^{(k)}({\bf x}) + g^{(k)}({\bf y}) \\
\text{s.t.}~&{\bf A}^{(k)}{\bf x}+{\bf B}^{(k)}{\bf y}={\bf b}^{(k)}\label{P2C}
\end{align}
\end{subequations}
where
\begin{align}
    & f^{(k)}({\bf x}) \hspace{-1mm}:=\hspace{-1mm} f_0^{(k)}({\bf x})+f_1^{(k)}({\bf x}), g^{(k)}({\bf y})\hspace{-1mm} := \hspace{-1mm} g_0^{(k)}({\bf y})+g_1^{(k)}({\bf y}),\nonumber\\
    & {\bf A}^{(k)} := {\bf A}(t_k),{\bf B}^{(k)} := {\bf B}(t_k),{\bf b}^{(k)} := {\bf b}(t_k),\nonumber\\
    & \mathcal{X}^{(k)}:=\mathcal{X}(t_k),\mathcal{Y}^{(k)}:=\mathcal{Y}(t_k).\nonumber
\end{align}
Throughout the paper we will assume that the following assumption holds.
\begin{assumption}
\label{Assum_uniform}
For each time $k$, $f^{(k)}, g^{(k)}$ satisfy
\vspace{-1mm}
\begin{align}
\hspace{-4mm}\langle\partial f^{(k)}({\bf x}_1)\hspace{-1mm}-\hspace{-1mm}\partial f^{(k)}({\bf x}_2), {\bf x}_1\hspace{-1mm}-\hspace{-1mm}{\bf x}_2\rangle&\hspace{-1mm}\geq\tilde{v}_f\|{\bf x}_1\hspace{-1mm}-\hspace{-1mm}{\bf x}_2\|, \forall {\bf x}_1,{\bf x}_2\label{eqn:convexity_f}\\
\hspace{-3mm}\langle\partial g^{(k)}({\bf y}_1)\hspace{-1mm}-\hspace{-1mm}\partial g^{(k)}({\bf y}_2), {\bf y}_1\hspace{-1mm}-\hspace{-1mm}{\bf y}_2\rangle&\geq\tilde{v}_g\|{\bf y}_1\hspace{-1mm}-\hspace{-1mm}{\bf y}_2\|,\forall {\bf y}_1,{\bf y}_2\label{eqn:convexity_g}
\end{align}
where $\tilde{v}_f,\tilde{v}_g$ are uniform lower bounds of strongly convex constants for $f^{(k)}, g^{(k)}$. Functions $f_1^{(k)}, g_1^{(k)}$ have Lipschitz-continuous gradients, i.e.,
\vspace{-1mm}
\begin{align}
    \|\nabla f_1^{(k)}({\bf x}_1)\hspace{-1mm} -\hspace{-1mm} \nabla f_1^{(k)}({\bf x}_2)\|&\leq \tilde{L}_f\|{\bf x}_1\hspace{-1mm}-\hspace{-1mm}{\bf x}_2\|,\forall {\bf x}_1,{\bf x}_2 \label{eqn:Lp_f}\\
    \|\nabla g_1^{(k)}({\bf y}_1) \hspace{-1mm}- \hspace{-1mm}\nabla g_1^{(k)}({\bf y}_2)\|&\leq\tilde{L}_g\|{\bf y}_1\hspace{-1mm}-\hspace{-1mm}{\bf y}_2\|,\forall{\bf y}_1,{\bf y}_2\label{eqn:Lp_g}
\end{align}
where $\tilde{L}_f,\tilde{L}_g$ are uniform upper bounds of Lipschitz constants for $\nabla f_1^{(k)}, \nabla g_1^{(k)}$.
\end{assumption}
\begin{assumption}\label{Assum_slater}
For each time $k$, the functions $f^{(k)}({\bf x}), g^{(k)}({\bf y})$ are coercive;  i.e.,
\begin{align*}
    f^{(k)}({\bf x})\rightarrow\infty~\text{as}~\|{\bf x}\|\rightarrow\infty,~g^{(k)}({\bf y})\rightarrow\infty~\text{as}~\|{\bf y}\|\rightarrow\infty \, .
\end{align*}
\end{assumption}
Assumption \ref{Assum_slater} will be instrumental to ensure that the iterates are bounded. Since continuous coercive functions' level sets $\{{\bf x}| f({\bf x}) \leq \mu_1, \forall \mu_1\}, \{{\bf y}| g({\bf y}) \leq \mu_2, \forall \mu_2\} $ are always compact, the optimal solutions to problem \eqref{P2}, defined as ${\bf x}^{\text{opt},(k)}, {\bf y}^{\text{opt},(k)}$, are bounded, i.e.,
\begin{align*}
    &\|{\bf x}^{\text{opt},(k)}\|\leq \sigma_1, ~\|{\bf y}^{\text{opt},(k)}\|\leq \sigma_2
\end{align*}
for some positive constants  $\sigma_1, \sigma_2$.

\section{Online Proximal-ADMM using Perturbations}

This section presents an ADMM-based algorithm to track an optimal solution trajectory of the time-varying problem \eqref{P1}. As summarized in Table I, the proposed algorithm exhibits linear convergence guarantees under less stringent  conditions relative to existing ADMM-based methods (even for static problems). In fact, although classic ADMM is conceptually simple and easy to implement, the conditions under which it is convergent is shown to be quite restrictive \cite{deng2016global}. We propose a new algorithm by leveraging the idea of dual perturbation \cite{hajinezhad2019perturbed,koshal2011multiuser} and gradient steps; this will provide a way to demonstrate convergence for a larger family of problems. However, a linear convergence rate at milder conditions comes at the cost of ensuring tracking of an approximate Karush-Kuhn-Tucker (KKT) point~\cite{hajinezhad2019perturbed,koshal2011multiuser,bernstein2018online}.

Accordingly, we propose to add a small perturbation to the dual variable ${\mbox{\boldmath $\lambda$}}$ in the form of $1-\beta\gamma$, where $\gamma>0$ is the perturbation parameter and $\beta \gamma \in(0,1)$. The perturbed augmented Lagrangian function is then defined as
\begin{align}
&\mathcal{L}^{(k)}({\bf x},{\bf y};{\mbox{\boldmath $\lambda$}}) =\mathcal{L}_1^{(k)}({\bf x},{\bf y};{\mbox{\boldmath $\lambda$}})+  f_0^{(k)}({\bf x})+g_0^{(k)}({\bf y}) \label{AL2}
\end{align}
where
\begin{align}
&\mathcal{L}_1^{(k)}({\bf x},{\bf y};{\mbox{\boldmath $\lambda$}})\hspace{-1mm} = \hspace{-1mm}f^{(k)}_1({\bf x})\hspace{-1mm}+\hspace{-1mm}g^{(k)}_1({\bf y}) +\frac{\beta}{2}\|{\bf A}^{(k)}{\bf x}\hspace{-1mm}+\hspace{-1mm}{\bf B}^{(k)}{\bf y}\hspace{-1mm}-\hspace{-1mm}{\bf b}^{(k)}\|^2\nonumber\\
&- (1-\beta\gamma){\mbox{\boldmath $\lambda$}}^T\hspace{-1mm}({\bf A}^{(k)}{\bf x}+{\bf B}^{(k)}{\bf y}\hspace{-1mm}-\hspace{-1mm}{\bf b}^{(k)}) \nonumber.
\end{align}
Mirroring~\cite{chambolle2011first},  to update ${\bf x}$ and ${\bf y}$, one can performs the following steps in an online fashion (where $k$ is here the time index):
\begin{align*}
{\bf y}^{(k+1)} \hspace{-1mm}&=\hspace{-1mm}\arg\min_{{\bf y}}\;\left\langle\frac{\partial\mathcal{L}_1^{(k+1)}({\bf x}^{(k)},{\bf y}^{(k)};{\mbox{\boldmath $\lambda$}}^{(k)})}{\partial{\bf y}^{(k)}},{\bf y}-{\bf y}^{(k)}\right\rangle\\
&+ g_0^{(k+1)}({\bf y})  + \frac{1}{2\alpha_2}\|{\bf y}-{\bf y}^{(k)}\|^2,\\
{\bf x}^{(k+1)} \hspace{-1mm}&=\hspace{-1mm}\arg\min_{{\bf x}}\; \left\langle\frac{\partial\mathcal{L}_1^{(k+1)}({\bf x}^{(k)},{\bf y}^{(k+1)};{\mbox{\boldmath $\lambda$}}^{(k)})}{\partial{\bf x}^{(k)}},{\bf x}-{\bf x}^{(k)}\right\rangle\\
&+ f_0^{(k+1)} ({\bf x}) + \frac{1}{2\alpha_1}\|{\bf x}-{\bf x}^{(k)}\|^2,
\end{align*}
where $\alpha_1,\alpha_2$ are step sizes. We are now ready to outline the \emph{online  proximal-ADMM algorithm} with perturbations, whose steps are the following (we stress again that $k$ is  the time index):
\begin{subequations}
\label{perturb}
\begin{align}
\hspace{-1cm}&\hspace{-2mm}{\bf y}^{(k+1)} \hspace{-1mm}= \hspace{-1mm}\underset{g_0^{(k+1)}}{\textrm{prox}}\hspace{-1mm}\left({\bf y}^{(k)} \hspace{-1mm}-\hspace{-1mm} \alpha_2\frac{\partial\mathcal{L}_1^{(k+1)}({\bf x}^{(k)},{\bf y}^{(k)};{\mbox{\boldmath $\lambda$}}^{(k)})}{\partial {\bf y}^{(k)}}\right),\label{itr:y}\\
&\hspace{-2mm}{\bf x}^{(k+1)} \hspace{-2mm}= \hspace{-1mm}\underset{f_0^{(k+1)}}{\textrm{prox}}\hspace{-1mm}\left({\bf x}^{(k)} \hspace{-2mm}- \hspace{-1mm}\alpha_1\frac{\partial\mathcal{L}_1^{(k+1)}({\bf x}^{(k)},{\bf y}^{(k+1)};{\mbox{\boldmath $\lambda$}}^{(k)})}{\partial {\bf x}^{(k)}}\hspace{-1mm}\right),\label{itr:x}\\
&\hspace{-2mm}{\mbox{\boldmath $\lambda$}}\hspace{-0.5mm}^{(k+1)} \hspace{-2mm}= \hspace{-1mm}(1\hspace{-1mm}-\hspace{-1mm}\beta\gamma){\mbox{\boldmath $\lambda$}}\hspace{-0.5mm}^{(k)}\hspace{-1mm} - \hspace{-1mm}\beta\hspace{-1mm}\left({\bf A}\hspace{-0.5mm}^{(k+1)}{\bf x}\hspace{-0.5mm}^{(k+1)}\hspace{-1.5mm}+\hspace{-1mm}{\bf B}\hspace{-0.5mm}^{(k+1)}{\bf y}\hspace{-0.5mm}^{(k+1)}\hspace{-1.5mm}-\hspace{-1mm}{\bf b}\hspace{-0.5mm}^{(k+1)}\hspace{-1mm}\right)\label{perturb_dual}.
\end{align}
\end{subequations}
The underlying assumption here is that for each time $k$, the interval $\tau$ is sufficient  to run at least one iteration of \eqref{perturb}.

Compared to classical ADMM-based algorithms (for both static and time-varying optimization),  key differences here are in the proximal gradient steps in the primal update and the perturbation added to ${\utwi{\lambda}}$. The proximal gradient steps may provide favorable computational gains when applied to a large-scale problem; it also facilitate ones to develop measurement-based algorithms as discussed in, e.g.,~\cite{dall2019optimization}. The perturbation added to ${\utwi{\lambda}}$ emerges when considering a regularized Lagrangian function of the form $\mathcal{L}^{(k)}({\bf x},{\bf y};{\utwi{\lambda}}) - \frac{\gamma}{2}\|{\utwi{\lambda}}\|^2$ as in e.g.,~\cite{bernstein2018online,koshal2011multiuser,hajinezhad2019perturbed}. This additional term renders the regularized Lagrangian strongly concave in ${\utwi{\lambda}}$. Adding a (small) perturbation in dual variable (or, equivalently, considering a regularized Lagrangian) is a very useful technique to ensure convergence of the ADMM, and even obtain a linear convergence behavior as explained shortly. To gain   intuition,  let us consider a toy example as follows:
\begin{align}
\min_{{\bf x}} \quad 0, \quad \mbox{s.t.} \quad {\bf A} {\bf x} =0
\end{align}
where ${\bf A}$ is some fixed matrix, not necessarily positive semidefinite.
The optimality condition for the above problem can be written down as the following saddle point problem
\begin{align}
\min_{{\bf x}}\max_{{\mbox{\boldmath $\lambda$}}}~ {\bf x}^T{\bf A}{\mbox{\boldmath $\lambda$}}\label{perturb_prob}.
\end{align}
One can apply the alternating gradient descent/ascent method for solving problem \eqref{perturb_prob}, whose steps are similar as \eqref{perturb} and are given below
\begin{align*}
{\bf x}^{(k+1)} &= {\bf x}^{(k)} - \alpha({\bf A}{\mbox{\boldmath $\lambda$}}^{(k)}),\\
{\mbox{\boldmath $\lambda$}}^{(k+1)} &= {\mbox{\boldmath $\lambda$}}^{(k)} + \beta({\bf A}^T{\bf x}^{(k)}).
\end{align*}
In Figure \ref{fig:perturb}, we plot ${\bf x}^T{\bf A}{\utwi{\lambda}}$ using a random matrix ${\bf A}$. An interesting observation  is that the algorithm will diverge if no perturbation is added to ${\bf y}$ as shown in Figure \ref{c1}; also see \cite{daskalakis2018limit} for a formal proof. However, once a small perturbation is added to ${\bf y}$ in both primal and dual updates, i.e.
\begin{align*}
{\bf x}^{(k+1)} &= {\bf x}^{(k)} - \alpha({\bf A}{\utwi{\lambda}}^{(k)}(1-\gamma\beta)),\\
{\mbox{\boldmath $\lambda$}}^{(k+1)} &= {\mbox{\boldmath $\lambda$}}^k(1-\gamma\beta) + \beta({\bf A}^T{\bf x}^{(k)}),
\end{align*}
where $\gamma>0$ is a small number, the algorithm will converge as shown in Figure \ref{c2}. This example serves as a motivation to use the perturbation technique.

\begin{figure}
        \centering
        \begin{subfigure}[b]{0.40\textwidth}
            \centering
            \includegraphics[width=\textwidth]{no_perturb-eps-converted-to.pdf}
            \caption[]%
            {{\small Performance without perturbation}}
            \label{c1}
        \end{subfigure}
        \hfill
        \begin{subfigure}[b]{0.40\textwidth}
            \centering
            \includegraphics[width=\textwidth]{with_perturb-eps-converted-to.pdf}
            \caption[]%
            {{\small Performance with perturbation.}}
            \label{c2}
        \end{subfigure}
        {\small }
        \caption{Example of trends of the objective value of \eqref{perturb_prob} for methods with and without perturbation.}
        \label{fig:perturb}
\end{figure}

\section{Convergence Analysis}

In this section we provide analytical results for convergence and tracking ability of the proposed algorithm.

From \cite{deng2016global}, it is known that existing ADMM has relatively strict conditions for linear convergence and these conditions may not hold true in some applications; for example, the coefficient matrices in constraints
\eqref{P2C} might not have full row rank (this is the case for the application presented later in the paper). Further, in some applications, the objective function of \eqref{P2} may also contain non-smooth terms, which can jeopardize the Lipschitz continuity property.
In contrast, the proposed algorithm could be utilized in a wider range of time-varying optimization problems. It is also worth pointing out that \cite{deng2016global} deals with static optimization problems; here, the focus is on time-varying settings. We begin by first making the following assumption.
\begin{assumption}
\label{Assum1}
For each time $k$, there exists a saddle point ${\bf w}^{\text{opt},(k)} = ({\bf x}^{\text{opt,(k)}},{\bf y}^{\text{opt,(k)}},{\utwi{\lambda}}^{\text{opt,(k)}})$ to problem \eqref{P2} that satisfies the KKT condition:
\begin{align*}
    ({\bf B}^{(k)})^T{\utwi{\lambda}}^{\text{opt},(k)} &\in \partial g^{(k)}({\bf y}^{\text{opt},(k)}),\\
    ({\bf A}^{(k)})^T{\utwi{\lambda}}^{\text{opt},(k)} &\in \partial f^{(k)}({\bf x}^{\text{opt},(k)}),\\
    {\bf A}^{(k)}{\bf x}^{\text{opt},(k)}+&{\bf B}^{(k)}{\bf y}^{\text{opt},(k)} = {\bf b}^{(k)},
\end{align*}
where ${\utwi{\lambda}}^{\text{opt,(k)}}$ is dual variable associated with \eqref{P2C}. The optimal dual variable has a uniform bound, i.e. $\|{\utwi{\lambda}}^{\text{opt},(k)}\|\leq \mathcal{M}$, where $\mathcal{M}$ is a constant.
\end{assumption}
Assumption \ref{Assum1} is a standard assumption to for convergence analysis \cite{deng2016global}. If Assumption \ref{Assum1} does not hold at a time  $k$, the problem formulation would not be well posed  since there is no solution trajectory to track.  \\
Next,  we analyze the convergence of the algorithm. To proceed, we concatenate primal and dual optimizer as $\{{\bf w}^*\} = \{{\bf x}^*;{\bf y}^*;{\mbox{\boldmath $\lambda$}}^*\}$ (for static case) as the optimizer of $\max\limits_{{\utwi{\lambda}}}\min\limits_{{\bf x},{\bf y}}\mathcal{L}^{(k)}$ at time $k$. For notation simplicity we neglect superscript $k$ for static case and we have:
\begin{subequations}\label{AKKT}
\begin{align}
{\bf A}^T{\mbox{\boldmath $\lambda$}}^*-\nabla f_1({\bf x}^*)\in\partial f_0({\bf x}^*)\label{x_opt}\\
{\bf B}^T{\mbox{\boldmath $\lambda$}}^*- \nabla g_1({\bf y}^*) \in \partial g_0({\bf y}^*) \label{y_opt}\\
{\bf A}{\bf x}^*+{\bf B}{\bf y}^*-{\bf b}+\gamma{\mbox{\boldmath $\lambda$}}^* = 0\label{dual_opt}.
\end{align}
\end{subequations}
Condition \eqref{AKKT} is a perturbed version of KKT conditions, related to approximate KKT (AKKT)\cite{dutta2013approximate,andreani2011sequential}. Basically, optimizer ${\bf w}^*$ is not necessarily the KKT point of original problem \eqref{P1}, but rather an approximate solution.
Let $\{{\bf x}^{\text{opt},(k)}, {\bf y}^{\text{opt},(k)}, {\utwi{\lambda}}^{\text{opt},(k)}\}$ be a KKT point of problem \eqref{P1}, which is also the solution to \eqref{AKKT} when $\gamma = 0$. From Assumption \ref{Assum1}, we know that any optimal dual solution ${\utwi{\lambda}}^{\text{opt},(k)}$ is bounded. Let ${\bf v}^{\text{opt},(k)} = \{{\bf x}^{\text{opt},(k)}, {\bf y}^{\text{opt},(k)}\}, {\bf v}^{*,(k)} = \{{\bf x}^{*,(k)}, {\bf y}^{*,(k)}\}$; we can then directly use the result in \cite[Proposition 3.1]{koshal2011multiuser} to show that the distance between ${\bf v}^{*,(k)}$ and ${\bf v}^{\text{opt},(k)}$ is bounded. that is,\footnote{The original result is for smooth strongly convex function, but it is easy to check that it still holds true to our problem.},
\begin{align}\label{nedic}
    \kappa\|{\bf v}^{*,(k)} - {\bf v}^{\text{opt},(k)}\|^2 + \frac{\gamma}{2}\|{\utwi{\lambda}}^{*,(k)}\|^2 \leq \frac{\gamma}{2}\|{\utwi{\lambda}}^{\text{opt},(k)}\|^2,
\end{align}
where $\kappa>0$ is a constant. We can further derive that
\begin{align}\label{nedic2}
    \|{\bf v}^{*,(k)} - {\bf v}^{\text{opt},(k)} \|\leq \sqrt{\gamma}\max\|{\utwi{\lambda}}^{\text{opt},(k)}\|\cdot\textrm{c}, \; \forall~k
\end{align}
where $c>0$ is some constant; see \cite[Proposition 3.1]{koshal2011multiuser} for a detailed discussion and the  proof of the result. Since every optimal dual variable ${\utwi{\lambda}}^{\text{opt},(k)}$ to the original problem \eqref{P1} is bounded, i.e. $\|{\utwi{\lambda}}^{\text{opt},(k)}\|\leq \mathcal{M}$, it follows that the distance between the AKKT point and a KKT point is bounded too, and such a bound depends on the choice of $\gamma$ (the smaller the $\gamma$, the smaller the distance since $\|{\utwi{\lambda}}^{\text{opt},(k)}\|$ is independent of $\gamma$).
In the following analysis, we focus on bounding the distance between iterates generated by \eqref{perturb} and AKKT points.
To proceed, we first state the assumptions on ${\bf x}^{*,(k)}, {\bf y}^{*,(k)}$ and on problem parameters as follows.
\begin{assumption}
\label{Assum_Bounded}
The successive difference between AKKT points is bounded:
\begin{align}
\|{\bf x}^{*,(k+1)} - &{\bf x}^{*,(k)}\|\leq \sigma_{{\bf x}}, \,\,\,\, \|{\bf y}^{*,(k+1)} - {\bf y}^{*,(k)}\|\leq \sigma_{{\bf y}},\label{diff}
\end{align}
where ${\bf x}^{*,(k)},{\bf y}^{*,(k)}$ are primal optimizer of \eqref{AKKT} at time $k$; $\sigma_{\bf x}>0, \sigma_{\bf y} >0$ are some constants. Also, the variation of the problem parameters is bounded as:
\begin{align}
    &\|{\bf A}^{(k+1)} - {\bf A}^{(k)}\| \leq \sigma_{\bf A}, ~\|{\bf B}^{(k+1)} - {\bf B}^{(k)}\| \leq \sigma_{\bf B} \label{diff_ab}\\
    &\|{\bf A}^{(k)}\|\leq\tilde{\sigma}_{{\bf A}},~\|{\bf B}^{(k)}\|\leq\tilde{\sigma}_{{\bf B}},~\|{\bf b}^{(k+1)}-{\bf b}^{(k)}\|\leq\sigma_{\bf b}\label{parameter}
\end{align}
where $\sigma_{\bf A},\sigma_{\bf B},\tilde{\sigma}_{{\bf A}},\tilde{\sigma}_{{\bf B}},\sigma_{\bf b}$ are some given positive constants.
\end{assumption}
Assumption \ref{Assum_Bounded} is  common in time-varying optimization~\cite{dall2019optimization, simonetto2014double, fazlyab2016self,bernstein2018online,Shahrampour2018,hall2015online,madden2020bounds}; worst-case bounds for \eqref{diff} can be obtained assuming that  the  sets $\mathcal{X}^{(k)},\mathcal{Y}^{(k)}$ are compact uniformly in time. Another approach is to measure the  distance based on the \textit{optimal drift}, without assuming a specific bound; see, e.g.,~\cite{ling2014decentralized,cao2019dynamic}. The parameters  $\sigma_{{\bf x}}$ and $\sigma_{{\bf y}}$ quantify the maximum variation of the optimal solutions over two consecutive time steps. Since the paper deals with a tracking problem, conventional wisdom would suggest that better tracking performance can be achieved when  \eqref{P1} is not changing rapidly; this will be confirmed in the convergence results presented later (see also ~\cite{dall2019optimization}).

The main result of the paper is stated next. For notation simplicity, let $\{{\bf w}^{(k)}\}=\{{\bf x}^{(k)};{\bf y}^{(k)};{\utwi{\lambda}}^{(k)}\}$ be the iterates generated by~\eqref{perturb}, and  let  $\{{\bf w}^{*,(k)}\}=\{{\bf x}^{*,(k)},{\bf y}^{*,(k)},{\utwi{\lambda}}^{*,(k)}\}$ be an optimizer of $\max\limits_{{\utwi{\lambda}}}\min\limits_{{\bf x},{\bf y}}\mathcal{L}^{(k)}$.

\begin{theorem}\label{main}
Suppose that Assumptions~\ref{Assum1}--\ref{Assum_Bounded} hold for each time $k$.  Let ${\bf G} =
\textrm{diag}(\frac{1}{\alpha_1}{\bf I}, \frac{1}{\alpha_2}{\bf I}, \frac{1}{\beta}{\bf I})$ be a  positive definite matrix. Assume that the  step size $\beta$ and the perturbation constant $\gamma$ satisfy: $\beta\gamma + \beta \leq 1, \beta \leq 1$. Finally, assume that the step sizes satisfy the  following:
\begin{align*}
    0 & < \alpha_1 \leq \left((1+\beta\gamma)\tilde{\sigma}^2_{\bf A} + \frac{\tilde{L}_f^2}{\tilde{v}_f} \right)^{-1}  \nonumber \\
    0 & < \alpha_2 \leq \left(\frac{2\beta^2\max\tilde{\sigma}^4_{\bf B}}{\tilde{v}_g}+ \frac{\tilde{L}_g^2}{\tilde{v}_g} \right)^{-1} \, .
\end{align*}

Then, at every time $k$, the tracking error of the algorithm~\eqref{perturb} evolves as:
\begin{align}
	\|{\bf w}^{(k)} \hspace{-1mm}-\hspace{-1mm} {\bf w}^{*,(k)}\|_{\bf G} & \leq
	 \frac{1}{1+\delta} \|{\bf w}^{(k-1)}-{\bf w}^{*,(k-1)}\|_{\bf G} \nonumber \\
	 & \,\, + \frac{1}{1+\delta} \left(\frac{\sigma_{\bf x}^2}{\alpha_1} +\frac{\sigma_{\bf y}^2}{\alpha_2} + 2\sigma_{{\utwi{\lambda}}}^2\right)^{\frac{1}{2}}
	 \label{eq:linear_conv}
	\end{align}
where $ \sigma_{{\utwi{\lambda}}} :=\tilde{\sigma}_{{\bf A}}\sigma_{\bf x}+\tilde{\sigma}_{{\bf B}}\sigma_{\bf y}+\sigma_{\bf b}+\sigma_{\bf A}\mathcal{J}(\sigma_1) + \sigma_{\bf B}\mathcal{J}(\sigma_2)$, $\mathcal{J}(\sigma_1) = \sigma_1 + \sqrt{\gamma}\mathcal{M}c, \mathcal{J}(\sigma_2) = \sigma_2 + \sqrt{\gamma}\mathcal{M}c$ and  $\delta$ satisfies the condition:
\begin{align}
0 < \delta \leq \min\left(\frac{\tilde{v}_f}{(1+\beta\gamma)\tilde{\sigma}^2_{\bf A} + \frac{\tilde{L}_f^2}{\tilde{v}_f}},\frac{\tilde{v}^2_g}{4\beta^2\tilde{\sigma}^4_{\bf B} + 2\tilde{L}_g^2}, \beta\gamma\right).
\end{align}
\end{theorem}

\begin{corollary}
\label{linear}
Under the assumptions of Theorem~\ref{main}, one has the following asymptotic behavior for the tracking error:
\begin{align}
	\label{eq:result_theorem}
		\limsup_{k\rightarrow\infty} \|{\bf w}^{(k)} - {\bf w}^{*,(k)}\|_{\bf G} \leq
\frac{1}{\delta} \left(\frac{\sigma_{\bf x}^2}{\alpha_1} +\frac{\sigma_{\bf y}^2}{\alpha_2} + 2\sigma_{{\utwi{\lambda}}}^2\right)^{\frac{1}{2}} .
\end{align}
\end{corollary}

The proofs of Theorem~\ref{main} and Corollary~\ref{linear} are provided in the Appendix.\\
\textit{Remark.} We also provide a comment on the boundedness of the dual iterates; specifically, ${\utwi{\lambda}}^{(k)}$ can be bounded as follows:
\begin{align*}
    \|{\utwi{\lambda}}^{(k)}\| &= \|{\utwi{\lambda}}^{(k)} - {\utwi{\lambda}}^{*,(k)} + {\utwi{\lambda}}^{*,(k)}\| \\&\leq \|{\utwi{\lambda}}^{(k)} - {\utwi{\lambda}}^{*,(k)}\| + \|{\utwi{\lambda}}^{*,(k)}\| \, .
\end{align*}
From Corollary 1, one has that $\|{\utwi{\lambda}}^{(k)} - {\utwi{\lambda}}^{*,(k)}\|$ is bounded; further, from \cite[Proposition 3.1]{koshal2011multiuser} and \eqref{nedic} it follows that $\|{\utwi{\lambda}}^{*,(k)}\|\leq \|{\utwi{\lambda}}^{\text{opt},(k)}\|$, where ${\utwi{\lambda}}^{\text{opt},(k)}$ is the optimal dual solution (without perturbation), which is assumed to be bounded as $\|{\utwi{\lambda}}^{\text{opt},(k)}\|\leq\mathcal{M}$. Therefore, $\|{\utwi{\lambda}}^{(k)}\|$ is bounded too. The distance between ${\bf w}^{(k)}$ and the optimal solution ${\bf w}^{\text{opt},(k)} = \{{\bf x}^{\text{opt},(k)},{\bf y}^{\text{opt},(k)},{\utwi{\lambda}}^{\text{opt},(k)}\}$ can be bounded as follows:
\begin{align*}
    &\|{\bf w}^{(k)} - {\bf w}^{\text{opt},(k)}\|_{\bf G} = \|{\bf w}^{(k)} - {\bf w}^{*,(k)} \hspace{-1mm}+\hspace{-1mm} {\bf w}^{*,(k)}\hspace{-1mm} -\hspace{-1mm} {\bf w}^{\text{opt},(k)}\|_{\bf G} \\
   &\leq \|{\bf w}^{(k)} - {\bf w}^{*,(k)}\|_{\bf G} + \|{\bf w}^{*,(k)} - {\bf w}^{\text{opt},(k)}\|_{\bf G}\\
   &\stackrel{\eqref{eq:result_theorem}}{\leq}\frac{1}{\delta} \left(\frac{\sigma_{\bf x}^2}{\alpha_1} +\frac{\sigma_{\bf y}^2}{\alpha_2} + 2\sigma_{{\utwi{\lambda}}}^2\right)^{\frac{1}{2}} + \|{\bf w}^{*,(k)} - {\bf w}^{\text{opt},(k)}\|_{\bf G}.
\end{align*}
From \eqref{nedic2} and $\|{\utwi{\lambda}}^{\text{opt},(k)}\|\leq\mathcal{M}$, we know that
\begin{align*}
    &\|{\bf w}^{*,(k)} - {\bf w}^{\text{opt},(k)}\|^2_{\bf G}\\
    =&\|{\bf v}^{*,(k)} - {\bf v}^{\text{opt},(k)}\|^2_{(\frac{1}{\alpha_1};\frac{1}{\alpha_2})}+\|{\utwi{\lambda}}^{*,(k)} - {\utwi{\lambda}}^{\text{opt},(k)}\|^2_{\frac{1}{\beta}}\\
    \leq&\max(\frac{1}{\alpha_1}, \frac{1}{\alpha_2})\gamma\mathcal{M}^2+\frac{2}{\beta}\mathcal{M}^2.
\end{align*}
This eventually gives us
\begin{align*}
    &\|{\bf w}^{(k)} - {\bf w}^{\text{opt},(k)}\|^2_{\bf G} \\
    \leq& \frac{2}{\delta^2} \left(\frac{\sigma_{\bf x}^2}{\alpha_1} +\frac{\sigma_{\bf y}^2}{\alpha_2} + 2\sigma_{{\utwi{\lambda}}}^2\right) + 2 {\max(\frac{1}{\alpha_1}, \frac{1}{\alpha_2})\gamma\mathcal{M}^2+\frac{4}{\beta}\mathcal{M}^2} \, .
\end{align*}
The result~\eqref{eq:linear_conv} asserts that the proposed algorithm exhibits linear convergence with a contraction coefficient of $1/(1-\delta)$~\cite{dall2019optimization}; the evolution of the tracking error depends on the temporal variability of the optimal trajectory, which is bounded by the second term on the right-hand-side of~\eqref{eq:linear_conv}. It is worth pointing out that, if the problem~\eqref{P1} is static, then~\eqref{eq:linear_conv} boils down to
$$
\|{\bf w}^{(k)} \hspace{-1mm}-\hspace{-1mm} {\bf w}^{*}\|_{\bf G} \leq
	 \frac{1}{1+\delta}\|{\bf w}^{(k-1)}-{\bf w}^{*}\|_{\bf G}
$$
showing linear convergence of the proximal-ADMM method in batch optimization (${\bf w}^{*}$ is in this case the solution of the static problem).

The asymptotic result~\eqref{eq:result_theorem} matches existing results in online methods for time-varying optimization~\cite{bernstein2018online,dall2019optimization,dixit2019online}. In particular, the bound depends on the $\delta$ (which affects the contraction coefficient) and the maximum variation of the optimal trajectory over two consecutive time steps, and it shows  how the variation of the problem parameters and optimal solutions can affect the  tracking performance.

Although~\eqref{eq:result_theorem} asserts that the maximum tracking error is bounded, its tightness is to be investigated on a case-by-case basis (i.e., based on the particular evolution of the solution)\footnote{For example, ~\cite{madden2020bounds} showed that a bound of the form~\eqref{eq:result_theorem} is actually tightly met for online gradient and proximal-gradient descent  for a particular sequence of adversarial cost functions.}.

The following corollary presented  shows how to maximize $\delta$ (and, hence, how to minimize the worst-case tracking bound).

\begin{corollary}
\label{stepsize}
If the step sizes are selected as
\begin{align*}
    \alpha_1 =\hspace{-1mm} \left(\hspace{-1mm}(1+\beta\gamma)\tilde{\sigma}^2_{\bf A} + \frac{\tilde{L}_f^2}{\tilde{v}_f} \right)^{-1}\hspace{-2mm},
    \alpha_2 =\hspace{-1mm} \left(\frac{2\beta^2\max\tilde{\sigma}^4_{\bf B}}{\tilde{v}_g}+ \frac{\tilde{L}_g^2}{\tilde{v}_g} \right)^{-1}
\end{align*}
then one has that the contraction coefficient $\frac{1}{1+ \delta}$ can be computed using the following expression:
\begin{align}
\delta = \min\left(\frac{\tilde{v}_f}{(1+\beta\gamma)\tilde{\sigma}^2_{\bf A} + \frac{\tilde{L}_f^2}{\tilde{v}_f}},\frac{\tilde{v}^2_g}{4\beta^2\tilde{\sigma}^4_{\bf B} + 2\tilde{L}_g^2}, \beta\gamma\right).
\end{align}
\end{corollary}

Notice that, for example, if we specify $\beta=0.5$ and $ \gamma=1$, then $\alpha_1,\alpha_2, \delta$ depend only  on the problem itself; i.e,
\vspace{-2mm}
\begin{align}
    \delta &= \min\left(\frac{\tilde{v}_f}{\frac{3}{2}\tilde{\sigma}^2_{\bf A} + \frac{\tilde{L}_f^2}{\tilde{v}_f}},\frac{\tilde{v}^2_g}{\tilde{\sigma}^4_{\bf B} + 2\tilde{L}_g^2}, \frac{1}{2}\right),\label{eqn:delta1}\\
    \alpha_1 &= \frac{1}{\frac{3}{2}\tilde{\sigma}^2_{\bf A} + \frac{\tilde{L}_f^2}{\tilde{v}_f}}, ~\alpha_2 = \frac{1}{\frac{\max\tilde{\sigma}^4_{\bf B}}{2\tilde{v}_g}+ \frac{\tilde{L}_g^2}{\tilde{v}_g}}.\label{eqn:a12}
\end{align}
As long as one picks $\delta$ as in \eqref{eqn:delta1}, there exist suitable $\alpha_1,\alpha_2$ to ensure convergence (see \eqref{eqn:alpha1}--\eqref{eqn:alpha2} in the proof).

\section{Example of Motivating Applications}\label{simulation}

\subsection{Multi-area power grid optimization}\label{power}

In this section, we briefly outline an example in power grids. We consider a distribution network featuring distributed energy resources (DERs), and we apply the proposed methodology to drive the DER output powers to the solution of an optimization problem encapsulating voltage constraints and given performance objectives. We demonstrate that the proposed methodology is amenable to settings where the distribution system is partitioned in areas; each area is autonomously controlled, and it ``trades'' power with adjacent areas based on given economic objectives~\cite{autonomous_grids}. In contrast,  previous works in the context of real-time optimal power flow involve centralized algorithms~\cite{tang2017real,hauswirth2016projected} or algorithms with a gather-and-broadcast architecture~\cite{dall2018optimal}.

Similar to~\cite{autonomous_grids}, consider partitioning a power distribution network into $C$ clusters, and denote as ${\mbox{$\mathcal{C}$}}_{i}$ the set of electrical nodes within cluster $i = 1, \ldots, C$. Two clusters $i$ and $j$ are adjacent if there is at least an electrical node $i$ such that $i \in {\mbox{$\mathcal{C}$}}_{i}$ and $i \in {\mbox{$\mathcal{C}$}}_{j}$. Let ${\cal B}_{i,j} := {\mbox{$\mathcal{C}$}}_{i} \cap {\mbox{$\mathcal{C}$}}_{j}$ be the set of \emph{boundary nodes} connecting cluster $i$ to cluster $j$, and define ${\cal B}_{i} := \cup_{j\neq i} {\cal B}_{i,j}$. Further, let ${\cal I}_{i} := {\mbox{$\mathcal{C}$}}_{i} \backslash {\cal B}_{i}$ be the set of \emph{internal nodes} for cluster $i$. For future developments, let $N_{i} := |{\cal I}_{i}|$ be the number of internal nodes if cluster $i$, and let ${\cal N}_{i} \subset \{1, \ldots, C\}$ be the set of neighboring clusters of the $i$th one (i.e., cluster connected to the $i$th one).

Let ${\bf x}_{j}^{i} := [P_{j}^{i}, Q_{j}^{i}]^\textsf{T} \in \mathbb{R}^2$ collect the net active and reactive powers injected by DERs at the node $j \in {\cal I}_{i}$ of cluster $i$. Particularly, ${\bf x}_{j}^{i} $ can represent the powers injected by \emph{one} DER located at node $j$, or the aggregate net power injections of \emph{a group of} DERs located at node $j$ (e.g., a household with multiple controllable devices) and we stack the setpoints $\{{\bf x}_{j}^{i}\}_{ j \in {\cal I}_{i}}$ in the vector ${\bf x}^{i} \in \mathbb{R}^{2 N_{i}}$. If no controllable DERs are present at a given location, the corresponding vector  ${\bf x}_{j}^{i}$ is  set to $\mathbf{0}$\footnote{For notation simplicity, the model is outlined for balanced systems and for the case where one household/building with DERs is located at a node. However, the model can be trivially extended to multiphase networks~\cite{linModels} and for the case where multiple households/buildings with DERs are located at a node (at the cost of increasing the complexity of the notation).}. On the other hand, ${\utwi{\ell}}_{j}^{i} \in \mathbb{R}^2$ denotes the net non-controllable loads at node $j \in {\cal I}_{i}$, and ${\utwi{\ell}}^{i} \in \mathbb{R}^{2 N_{i}}$ stacks the loads $\{{\utwi{\ell}}_{j}^{i}\}_{ j \in {\cal I}_{i}}$. It is assumed that no DERs and no non-controllable loads are located at the boundary nodes ${\cal B}_{i,j}$.

Let $V_{j}^{i} \in \mathbb{C}$ denote the complex line-to-ground voltage phasor at node $j$ of cluster $i$, and let ${\bf v}^{i} := [\{|V_{j}^{i}|, j \in {\cal I}_{i}\}]^\textsf{T}$ be the vector of voltage  magnitudes of the internal nodes ${\cal I}_{i}$. For each pair of neighboring clusters $(i,j)$, let ${\bf x}_{n}^{j \rightarrow i} := [P_{n}^{j \rightarrow i}, Q_{n}^{j \rightarrow i}]^\textsf{T} \in \mathbb{R}^2$ represent the active and reactive powers flowing into area $i$ from area $j$ through node $n \in {\cal B}_{i,j} $; on the other hand, ${\bf x}_{n}^{i \rightarrow j} \in \mathbb{R}^2$ contains the active and reactive powers flowing into area $j$ from area $i$ through node $n \in {\cal B}_{i,j} $. From Kirchhoff's Law, it holds that ${\bf x}_{n}^{j \rightarrow i} + {\bf x}_{n}^{i \rightarrow j} = \mathbf{0}$. To facilitate the syntheses of computationally-affordable algorithms, we leverage the following approximate linear relationship between net injected power and voltage magnitude (see e.g.,~\cite{dall2018optimal} and references therein):
\begin{subequations} \label{relation1}
\begin{align}
\tilde{{\bf v}}^{i}  &:=  \hspace{-2mm} \sum_{j \in {\cal I}_{i}} {\bf A}_{j}^{i} ({\bf x}^{i}_j - {\utwi{\ell}}^{i}_j) +\hspace{-2mm}  \sum_{j \in {\cal N}_{i}} \sum_{n \in {\cal B}_{i,j}}\hspace{-2mm} {\bf A}_{n}^{j \rightarrow i} {\bf x}_{n}^{j \rightarrow i} \hspace{-1mm}+\hspace{-1mm} {\bf a}, \label{eqn:lin_v} \\
 &=  {\bf A}^{i} ({\bf x}^{i} - {\utwi{\ell}}^{i}) +  \sum_{j \in {\cal N}^{i}}  {\bf A}^{j \rightarrow i} {\bf x}^{j \rightarrow i} + {\bf a}, \label{eqn:lin_v2}
\end{align}
\end{subequations}
where ${\bf A}^i = [{\bf A}_j^i]_{j\in{\cal I}_i}, {\bf A}^{j\rightarrow i} = [{\bf A}^{j\rightarrow i}_n]_{n\in{\cal B}_{i,j}}, {\bf a}$ are time-varying problem parameters derived from linearized power flow equation.
Another linear relationship between net injected power and power between clusters is captured in the following equation:
\begin{subequations}\label{relation2}
\begin{align}
{\bf x}^{j \rightarrow i} &:=  \sum_{k \in {\cal I}_{i}} {\bf M}^{j \rightarrow i}_{k}  ({\bf x}^{i}_k - {\utwi{\ell}}^{i}_k)  + \mathbf{m}^{j \rightarrow i}, \label{eqn:lin_s0} \nonumber\\  + & \sum_{k \in {\cal N}_{i} \backslash \{j\}} \sum_{n \in {\cal B}_{i,k}} {\bf M}^{k, j \rightarrow i}_{n}{\bf x}_{n}^{k \rightarrow i}\\
&=   {\bf M}^{j \rightarrow i}  ({\bf x}^{i}\hspace{-1mm} - \hspace{-1mm}{\utwi{\ell}}^{i}) + \mathbf{m}^{j \rightarrow i}+ \hspace{-4mm} \sum_{k \in {\cal N}_{i} \backslash \{j\}} \hspace{-3mm}{\bf M}^{k, j \rightarrow i} {\bf x}^{j \rightarrow i}, \label{eqn:lin_s01}
\end{align}
\end{subequations}
where ${\bf M}^{j\rightarrow i} \hspace{-1mm}=\hspace{-1mm} [{\bf M}^{j\rightarrow i}_k]_{k\in{\cal I}_i}, {\bf M}^{k,j\rightarrow i} \hspace{-1mm}=\hspace{-1mm} [{\bf M}^{k,j\rightarrow i}_{n}]_{n\in{\cal B}_{i,k}}$,${\bf m}^{j\rightarrow i}$ are also time-varying problem parameters depending on the actual network physics.
All model parameters in \eqref{relation1}--\eqref{relation2} can be obtained as shown in~\cite{linModels}.
Now we are ready to state our real-time OPF problem as follows:
\begin{subequations}
\begin{align}
&\hspace{-5mm}\min_{\substack{\{{\bf x}^{i}\}, \{{\bf x}_{n}^{j \rightarrow i}, {\bf x}_{n}^{i \rightarrow j}\}}}  \hspace{.2cm} \sum_{i=1}^C  [f^{i}({\bf x}^{i}) + g^i(\{{\bf x}^{j \rightarrow i}\})] \tag{P3} \label{obj} \\
 \text{s.t.} ~&{\bf x}_{j}^{i} \in {\cal Y}_{j}^{i}, \forall \,\, j \in {\cal I}^{i}, \,\, i = 1, \ldots C \label{eqn:constr_Y_exact}  \\
& v^{\text{min}} \mathbf{1} \leq \tilde{{\bf v}}^{i} \leq v^{\text{max}} \mathbf{1}, \forall \,\, i = 1, \ldots C \label{cons}\\
& {\bf x}^{j \rightarrow i} =   {\bf M}^{j \rightarrow i}  ({\bf x}^{i} - {\utwi{\ell}}^{i}) + \mathbf{m}^{j \rightarrow i}+\hspace{-5mm}  \sum_{k \in {\cal N}_{i} \backslash \{j\}} \hspace{-5mm}{\bf M}^{k, j \rightarrow i} {\bf x}^{j \rightarrow i} \nonumber\\
&\hspace{3cm}, \, \forall \, j \in {\cal N}_{i}, \, \, i = 1, \ldots, C \label{between}\\
& {\bf x}^{j \rightarrow i} + {\bf x}^{i \rightarrow j} = \mathbf{0} , \,\, \forall \,\, \textrm{neighboring~areas~} (i,j)
\end{align}
\end{subequations}
where the time-varying objective function models the amount of real power curtailed and the amount of reactive power injected or absorbed (which leads to non-smooth term in the objective, e.g. $\ell_1$ term). For notation simplicity, we write objective function in \eqref{obj} as $\Psi({\bf x})$. Consider ${\bf M}^{j\rightarrow i}$ consists of $1,0$, with $1$ for real power, $0$ for reactive power. Putting \eqref{eqn:lin_s01} back to \eqref{eqn:lin_v2}, adding slack variables ${\utwi{\gamma}}^i,{\utwi{\beta}}^i$ to \eqref{cons} formulate equality constraints, and adding strongly convex term w.r.t ${\utwi{\gamma}} = \{{\utwi{\gamma}}^i\},{\utwi{\beta}} = \{{\utwi{\beta}}^i\}$ we have the following formulation:
\begin{subequations}
\begin{align}
 &\hspace{-0.5cm}\min_{\substack{\{\{{\bf x}^{i}\}, \{{\bf x}_{n}^{j \rightarrow i}, {\bf x}_{n}^{i \rightarrow j}\}\},\{{\utwi{\gamma}}^{i},{\utwi{\beta}}^{i}\geq 0\}\} } } \Psi({\bf x})+ a\|{\utwi{\gamma}}\|^2 + b\|{\utwi{\beta}}\|^2\tag{P4} \label{obj2} \\
 \text{s.t.~} & {\bf x}_{j}^{i} \in {\cal Y}_{j}^{i},  \forall \,\, j \in {\cal I}_{i}, \,\, i = 1, \ldots C \label{eqn:constr_Y_exact}  \\
& v^{\text{min}}\mathbf{1} -\tilde{{\bf v}}^{i} + {\utwi{\gamma}}^{i}= \mathbf{0},  \forall \,\, i = 1, \ldots C \label{cons2}\\
&\tilde{{\bf v}}^{i} + {\utwi{\beta}}^{i} - v^{\text{max}}\mathbf{1}=\mathbf{0} ,  \forall \,\, i = 1, \ldots C \label{cons3}\\
& {\bf x}^{j \rightarrow i} =   {\bf M}^{j \rightarrow i}  ({\bf x}^{i} - {\utwi{\ell}}^{i}) + \mathbf{m}^{j \rightarrow i}+\hspace{-2mm}  \sum_{k \in {\cal N}_{i} \backslash \{j\}}\hspace{-2mm} {\bf M}^{k, j \rightarrow i} {\bf x}^{j \rightarrow i} \nonumber\\
&\hspace{3cm}, \, \forall \, j \in {\cal N}_{i}, \, \, i = 1, \ldots, C \label{eqn:between}\\
& {\bf x}^{j \rightarrow i} + {\bf x}^{i \rightarrow j} = \mathbf{0} , \,\, \forall \,\, \textrm{neighboring~areas~} (i,j).\label{eqn:consensus}
 \end{align}
\end{subequations}
We can now clearly see a mapping from \eqref{obj2} to \eqref{P1}: objective functions
are $\Psi({\bf x})$ and $a\|{\utwi{\gamma}}\|^2 + b\|{\utwi{\beta}}\|^2$(where $a,b > 0$ are small); two blocks of variables are $\{{\bf x}^{i},{\bf x}^{j\rightarrow i}\}$ and $\{{\utwi{\gamma}}^{i},{\utwi{\beta}}^{i}\}$; constraints are all linear and separable w.r.t each network node. \textcolor{black}{Problem \eqref{obj2} is time varying in both objective function and constraint parameters.} In order to better illustrate how the proposed algorithm can be applied, we use a 4-cluster network (see Figure \ref{fig:net}) as an example. First, we substitute ${\bf x}^{j\rightarrow i}$ in \eqref{eqn:lin_v2} with \eqref{eqn:between}; then, we substitute $\tilde{{\bf v}}^i$ in \eqref{cons2}--\eqref{cons3} with \eqref{eqn:lin_v2};
last, we define the corresponding augmented Lagrangian function as follows:
\vspace{-2mm}
\begin{align*}
&\mathcal{L}({\bf x},{\utwi{\gamma}},{\utwi{\beta}},{\utwi{\lambda}}) = \Psi({\bf x}) + a\|{\utwi{\gamma}}\|^2 + b\|{\utwi{\beta}}\|^2 \nonumber\\
+ &\sum_{k}\hspace{-1mm}\sum_{i\in C_k}\frac{\rho}{2}\|\textcolor{black}{\sum_{j\in C_k}({\bf A}^{i}_j+{\bf A}_j^{3}){\bf x}^{i}_j  \hspace{-1mm}+\hspace{-1mm} {\bf a}_k }\hspace{-1mm}+\hspace{-1mm}{\utwi{\beta}}^{i} \hspace{-1mm}-\hspace{-1mm} v^{\text{max}} \hspace{-1mm}+ \hspace{-1mm}\frac{{\utwi{\lambda}}_1(1-\gamma)}{\rho}\|^2\\
+&\sum_{k}\hspace{-1mm}\sum_{i\in C_k} \frac{\rho}{2}\|v^{\text{min}} \hspace{-1mm}\textcolor{black}{-\hspace{-1mm}\sum_{j\in C_k} ({\bf A}^{i}_j\hspace{-1mm}+\hspace{-1mm}{\bf A}^{3}_j){\bf x}^{i}_j \hspace{-1mm}- \hspace{-1mm}{\bf a}_k} \hspace{-1mm}+\hspace{-1mm} {\utwi{\gamma}}^{i}\hspace{-1mm}+\hspace{-1mm}\frac{{\utwi{\lambda}}_2(1-\gamma)}{\rho}\|^2\\
+& \sum_{k}\frac{\rho}{2}\|\sum_{i}{\bf x}^{i\rightarrow j} - \textcolor{black}{\sum_{j\in C_k} {\bf x}^{j}} + \frac{{\utwi{\lambda}}_4(1-\gamma)}{\rho}\|^2\\
+&\frac{\rho}{2}\|\sum_{i\neq j}{\bf x}^{i\rightarrow j}+\frac{{\utwi{\lambda}}_3(1-\gamma)}{\rho}\|^2
\end{align*}

The detailed updates follow the same way as \eqref{perturb} and from Table I we know that linear convergence to AKKT is guaranteed. To further improve our algorithm for this particular application. We incorporate system measurements in both primal and dual updates in the following way:
\begin{align}
\sum_{j\in C_k}({\bf A}^{i}_j+{\bf A}_j^{3}){\bf x}^{i}_j +{\bf a}_k &\rightarrow \phi({\bf x}),~ \sum_{j\in C_k}{\bf x}^{j}\rightarrow \psi({\bf x}),
\end{align}
where $\phi({\bf x}), \psi({\bf x})$ are measurements. This is beneficial in that: i) A natural distributed computing scheme is achieved while without feedback it is not clear whether the algorithm can be implemented in a distributed way; ii) feedback terms are much less than uncontrollable terms, which essentially shrinks the measuring time; iii) it is easier to satisfy power flow equations with the help of system measurements.

\begin{figure}
\centering
\includegraphics[width=8cm]{net.png}
\caption{power distribution network with 4 clusters, node 3 is a boundary node that belongs to all 4 clusters.}
\label{fig:net}
\end{figure}

\subsection{Simulations}
In this section, we test our algorithm using the same power systems settings. We consider a similar system as in \cite{dall2018optimal}, where a modified IEEE 37-node test feeder is utilized. The network is obtained by considering a single phase equivalent, and by replacing the loads on phase ``c" specified in the original dataset with real load data measured from feeders in a neighborhood called Anatolia in California during a week in August 2012. It is assumed that the aggregations of photovoltaic systems are located at nodes 4, 7, 10, 13, 17, 20, 22, 23, 26, 28, 29, 30, 31, 32, 33, 34, 35, and 36. The rating of  these inverters are 300kVA for $i=3$, 350kVA for $i=15, 16$ and 200kVA for the remaining ones. The objective is set to be $f^i({\bf x}^i) = c_p(P_{\text{av},i} - P_i)^2 + c_q(Q_i)^2 + \bar{c}_q|Q_i|, g^i({\bf x}^{j\rightarrow i}) = 0$ where $P_{\text{av},i}$ is the maximum real power available from the PV system $i$, and $c_p = 3, c_q = 1, \bar{c}_q = 0.1$. The voltage limits are set to be $V^{\text{min}} = 0.95$pu, $V^{\text{max}} = 1.05$pu. The generation profiles are simulated based on real solar irradiance data and have a granularity of 1 second.
\begin{figure}
	\centering
	\includegraphics[width=8cm]{a2_p-eps-converted-to.pdf}
	\caption{Real power at feeder head during 12:00-12:30.}
	\vspace{-.2cm}
	\label{f2}
\end{figure}
First we specify a given trajectory for the power at the common coupling, which is color-coded in red in Fig. \ref{f2} (negative power indicates reverse power flows). It can be seen that our algorithm is able to regulate $P_0^k$ close to $P_{0,\text{set}}^k$ in real time. Figure \ref{f3} illustrates the voltage profiles for selected nodes. From 10:00 to 12:00 we observe a few flickers, which is caused by rapid variations of the solar irradiance. Other than that, it can be seen that voltage regulation is enforced and a flat voltage profile is obtained. Note that even there are some relatively large jumps from around 12:00 to 14:00, our algorithm is still able to track the optimal trajectory. A comparison with double smoothing algorithm \cite{dall2018optimal} is presented in Figure \ref{f4}. The proposed strategy has potentially better voltage regulation ability, especially for extreme cases e.g. the two spikes from 10:00 to 12:00.
\begin{figure}
	\centering
	\includegraphics[width=8cm]{a2_v-eps-converted-to.pdf}
	\caption{Voltage  profile achieved (only some nodes are considered for illustration purposes). }
	\vspace{-.2cm}
	\label{f3}
\end{figure}
\begin{figure}
	\centering
	\includegraphics[width=8cm]{voltage_sum_measurements-eps-converted-to.pdf}
	\caption{Index for the overall voltage violation across the system $\sum\limits_{n\in\mathcal{N}}\left(\max(|V_n^k| - v^{\text{max}} ,0)+\max(v^{\text{min}} - |V_n^k|,0)\right)$}
	\vspace{-.2cm}
	\label{f4}
\end{figure}
\begin{figure}[t]
        \centering
        \includegraphics[width=8cm]{cons_0916_c1-eps-converted-to.pdf}
        \caption{Power violation of Cluster 1: power violation for each cluster is defined as $\|\sum\limits_{i\in\mathcal{N}^{(j)}}{\bf x}^{i\rightarrow j} - \sum\limits_{j\in C}{\bf x}^{j}\|^2$}
        \label{power_vio}
    \end{figure}%
\begin{figure}[t]
        \centering
        \includegraphics[width=8cm]{cons_0916_cv-eps-converted-to.pdf}
        \caption{Consensus violation: $\|{\bf x}^{(i\rightarrow j)}+{\bf x}^{(j\rightarrow i)}\|^2$}
        \label{cons_power}
    \end{figure}
We proceed to test in the same setting except we are adding consensus constraints. In Figure \ref{power_vio} we can see that for all 4 clusters, power violation decreases dramatically in first a few minutes and remains at a  low level of $10^{-10}$. The power consensus violation is shown in Figure \ref{cons_power}, where a steep drop at the begging and flat low line after that are observed.

\section{Conclusion}
This paper gives a general online optimization problem formulation and proposes a online algorithm based on alternating direction method of multipliers that can continuously track optimal solution in real time. The steps of ADMM are proximal gradient steps with modification of adding perturbation to dual variable and incorporating system feedback for certain applications. The resulting algorithm is proved to converge to a neighborhood of optimal solution for each time instance. Numerical results for power systems applications also demonstrate the practicality of the proposed algorithm. Our future research will focus on general online nonconvex optimization problems.

\appendices

\subsection{Proof of Theorem \ref{main}}

To show the result of Theorem~\ref{main}, we start from the following lemma.

\begin{lemma}
\label{lem:linear}
Under the assumptions of Theorem~\ref{main}, it holds that:
\begin{align}
\label{eqn:lin}
\|{\bf w}^{(k)}-{\bf w}^{*,(k)}\|^2_{\bf G} \geq (1+\delta)\|{\bf w}^{(k+1)}-{\bf w}^{*,(k)}\|^2_{\bf G} .
\end{align}

\end{lemma}

\begin{proof}
Since the lemma focuses on a particular time instant $k$, we replace ${\bf w}^{*,(k)}$ with ${\bf w}^{*}$ for notation simplicity. Also, to better fit the equations in the columns, we replace the subscript $^{(k)}$ and $^{(k+1)}$ with $^{k}$ and $^{k+1}$, respectively, in the proof.

From the optimality condition of subproblem \eqref{itr:y}, one has that:
\begin{align}
&{\bf A}^T{\mbox{\boldmath $\lambda$}}^k(1-\beta\gamma) - \beta {\bf A}^T({\bf A}{\bf x}^{k} + {\bf B}{\bf y}^{k+1}-{\bf b})\nonumber\\-&\nabla f_1({\bf x}^k) + \frac{1}{\alpha_1}{\bf I}({\bf x}^k\hspace{-1mm}-{\bf x}^{k+1})\in\partial f_0({\bf x}^{k+1}).\nonumber
\end{align}
Rearrange terms we get
\begin{align}
& {\bf A}^{T}{\mbox{\boldmath $\lambda$}}^{k+1}\hspace{-1mm}+\hspace{-1mm}\beta{\bf A}^T{\bf A}({\bf x}^{k+1}\hspace{-1mm} - \hspace{-1mm}{\bf x}^{k})\hspace{-1mm} +\hspace{-1mm} \nabla f_1({\bf x}^{k+1})\hspace{-1mm} -\hspace{-1mm} \nabla f_1({\bf x}^k) \nonumber\\
+& \frac{1}{\alpha_1}{\bf I}({\bf x}^k-{\bf x}^{k+1})\in\partial f_0({\bf x}^{k+1}) + \nabla f_1({\bf x}^{k+1}).\label{x_k}
\end{align}
Also from optimality condition of subproblem and \eqref{itr:x}, one has that
\begin{align*}
&{\bf B}^T{\mbox{\boldmath $\lambda$}}^k(1\hspace{-1mm}-\hspace{-1mm}\beta\gamma)\hspace{-1mm} - \hspace{-1mm} \beta {\bf B}^{T}({\bf A}{\bf x}^{k+1}\hspace{-1mm}+\hspace{-1mm} {\bf B}{\bf y}^{k+1}\hspace{-1mm} -\hspace{-1mm} {\bf b})\hspace{-1mm}\nonumber\\
+&\beta{\bf B}^T{\bf A}({\bf x}^{k+1}-{\bf x}^k) + \beta{\bf B}^T{\bf B}({\bf y}^{k+1}-{\bf y}^k)\nonumber\\
+& \frac{1}{\alpha_2}{\bf I}({\bf y}^k\hspace{-1mm} -\hspace{-1mm} {\bf y}^{k+1}) -\nabla g_1({\bf y}^k)\in\partial g_0({\bf y}^{k+1}) .\nonumber
\end{align*}
Rearrange terms we get
\begin{align}
& {\bf B}^T({\mbox{\boldmath $\lambda$}}^{k+1} + \beta {\bf A}({\bf x}^{k+1}-{\bf x}^k) +\beta{\bf B}({\bf y}^{k+1}-{\bf y}^k))\nonumber\\
&+ \nabla g_1({\bf y}^{k+1}) - \nabla g_1({\bf y}^k)+ \frac{1}{\alpha_2}{\bf I}({\bf y}^k-{\bf y}^{k+1}) \nonumber\\
\in &\partial g_0({\bf y}^{k+1}) + \nabla g_1({\bf y}^{k+1})\label{y_k}.
\end{align}
Furthermore, from the dual update \eqref{perturb_dual}, one can obtain:
\begin{align}
\frac{1}{\beta}({\mbox{\boldmath $\lambda$}}^k - {\mbox{\boldmath $\lambda$}}^{k+1}) = \gamma{\mbox{\boldmath $\lambda$}}^k + ({\bf A}{\bf x}^{k+1} + {\bf B}{\bf y}^{k+1} - {\bf b}),\label{dual}
\end{align}
and, together with optimality condition, one obtains:
\begin{align}\label{dual_opt2}
\frac{1}{\beta}({\mbox{\boldmath $\lambda$}}^k\hspace{-1mm} - \hspace{-1mm}{\mbox{\boldmath $\lambda$}}^{k+1})\hspace{-1mm} =\hspace{-1mm} \gamma({\mbox{\boldmath $\lambda$}}^k\hspace{-1mm}-\hspace{-1mm}{\mbox{\boldmath $\lambda$}}^*) \hspace{-1mm}+\hspace{-1mm} {\bf A}({\bf x}^{k+1}\hspace{-2mm}-\hspace{-1mm}{\bf x}^{*})\hspace{-1mm} + \hspace{-1mm}{\bf B}({\bf y}^{k+1}\hspace{-2mm}-\hspace{-1mm}{\bf y}^{*})
\end{align}
Since the functions $f = f_0 + f_1$ and $g = g_0+g_1$ are strongly convex, we leverage \eqref{eqn:convexity_f} and \eqref{eqn:convexity_g}
and, by plugging the optimality condition \eqref{x_k} and \eqref{y_k}, we have
\begin{align}
&\langle{\bf A}^{T}({\mbox{\boldmath $\lambda$}}^{k+1}\hspace{-1mm}-\hspace{-1mm}{\utwi{\lambda}}^{*})+\beta{\bf A}^T{\bf A}({\bf x}^{k+1}\hspace{-1mm} -\hspace{-1mm} {\bf x}^{k})\hspace{-1mm} +\hspace{-1mm} \nabla f_1({\bf x}^{k+1}) \hspace{-1mm}- \hspace{-1mm}\nabla f_1({\bf x}^k)\nonumber\\
&+ \frac{1}{\alpha_1}{\bf I}({\bf x}^k-{\bf x}^{k+1}),{\bf x}^{k+1}-{\bf x}^*\rangle\geq v_f\|{\bf x}^{k+1}-{\bf x}^*\|^2\label{x}\\
&\langle{\bf B}^T({\mbox{\boldmath $\lambda$}}^{k+1}-{\utwi{\lambda}}^* + \beta {\bf A}({\bf x}^{k+1}-{\bf x}^k) +\beta{\bf B}({\bf y}^{k+1}-{\bf y}^k)) \nonumber\\
&+ \nabla g_1({\bf y}^{k+1}) - \nabla g_1({\bf y}^k)+\frac{1}{\alpha_2}{\bf I}({\bf y}^k-{\bf y}^{k+1}),{\bf y}^{k+1}-{\bf y}^*\rangle \nonumber\\
&\geq v_g\|{\bf y}^{k+1}-{\bf y}^*\|^2\label{y}
\end{align}
Next, add \eqref{x} and \eqref{y} together, define $\Phi = v_f\|{\bf x}^{k+1}-{\bf x}^*\|^2+v_g\|{\bf y}^{k+1}-{\bf y}^*\|^2$ and plug in \eqref{dual_opt2} to obtain:
\begin{align}
&\langle\beta{\bf A}({\bf x}^{k+1}\hspace{-1mm}-\hspace{-1mm}{\bf x}^k)\hspace{-1mm}+\hspace{-1mm}{\utwi{\lambda}}^{k+1}\hspace{-1mm}-\hspace{-1mm}{\utwi{\lambda}}^*,\frac{1}{\beta}({\utwi{\lambda}}^k\hspace{-1mm}-\hspace{-1mm}{\utwi{\lambda}}^{k+1}))\hspace{-1mm}-\hspace{-1mm}\gamma({\utwi{\lambda}}^k\hspace{-1mm}-\hspace{-1mm}{\utwi{\lambda}}^*) \rangle \nonumber\\
+ &\langle\frac{1}{\alpha_1}({\bf x}^k\hspace{-1mm}-\hspace{-1mm}{\bf x}^{k+1})+\nabla f_1({\bf x}^{k+1})-\nabla f_1({\bf x}^k), {\bf x}^{k+1}-{\bf x}^*\rangle \nonumber\\
+& \langle\nabla g_1({\bf y}^{k+1})\hspace{-1mm}-\hspace{-1mm}\nabla g_1({\bf y}^k)+\beta{\bf B}^T{\bf B}({\bf y}^{k+1}\hspace{-1mm}-\hspace{-1mm}{\bf y}^k), {\bf y}^{k+1}\hspace{-1mm}-\hspace{-1mm}{\bf y}^* \rangle \nonumber\\
+&\langle \frac{1}{\alpha_2}({\bf y}^k-{\bf y}^{k+1}), {\bf y}^{k+1}-{\bf y}^*\rangle
\geq \Phi. \label{cross_all}
\end{align}
We then proceed to split the cross terms and one can notice that there are similar terms for ${\bf x},{\bf y},{\utwi{\lambda}}$ in the following form:
\begin{align*}
    &\langle{\utwi{\lambda}}^{k+1}-{\utwi{\lambda}}^*,\frac{1}{\beta}({\utwi{\lambda}}^k-{\utwi{\lambda}}^{k+1})\rangle + \langle \frac{1}{\alpha_1}({\bf x}^k\hspace{-1mm}-\hspace{-1mm}{\bf x}^{k+1}), {\bf x}^{k+1}\hspace{-1mm}-\hspace{-1mm}{\bf x}^*\rangle\\
    &+\frac{1}{\alpha_2}({\bf y}^k-{\bf y}^{k+1}), {\bf y}^{k+1}-{\bf y}^*\rangle
\end{align*}
We group them together and define the following quantities:
$${\bf G} = \left(\begin{matrix}
\frac{1}{\alpha_1}{\bf I}&0&0\\0&\frac{1}{\alpha_2}{\bf I}&0\\0&0&\frac{1}{\beta}{\bf I}
\end{matrix}\right),{\bf w}^k = \left(\begin{matrix}
{\bf x}^k\\{\bf y}^k\\ {\mbox{\boldmath $\lambda$}}^k
\end{matrix}\right) , \,\, {\bf w}^* = \left(\begin{matrix}
{\bf x}^*\\{\bf y}^*\\ {\mbox{\boldmath $\lambda$}}^*
\end{matrix}\right),$$
so that one can rewrite the inequality as follows:
\begin{align}
&({\bf w}^{k+1} - {\bf w}^*)^T{\bf G}({\bf w}^{k} - {\bf w}^{k+1}) + \gamma\langle {\mbox{\boldmath $\lambda$}}^*-{\mbox{\boldmath $\lambda$}}^k, {\mbox{\boldmath $\lambda$}}^{k+1} - {\mbox{\boldmath $\lambda$}}^*\rangle\nonumber\\
+&\langle {\mbox{\boldmath $\lambda$}}^{k}\hspace{-1mm} - \hspace{-1mm}{\mbox{\boldmath $\lambda$}}^{k+1}, {\bf A}({\bf x}^{k+1}-{\bf x}^{k})\rangle + \beta\gamma\langle {\mbox{\boldmath $\lambda$}}^*\hspace{-1mm}-\hspace{-1mm}{\mbox{\boldmath $\lambda$}}^k, {\bf A}({\bf x}^{k+1}-{\bf x}^{k})\rangle \nonumber\\
+&\langle\beta{\bf B}^T{\bf B}({\bf y}^{k+1} - {\bf y}^k),{\bf y}^{k+1}-{\bf y}^* \rangle\nonumber\\
+&\langle{\bf x}^{k+1}\hspace{-1mm}-\hspace{-1mm}{\bf x}^*,\nabla f_1(x^{k+1})\hspace{-1mm} -\hspace{-1mm} \nabla f_1({\bf x}^k)\rangle\nonumber\\
+&\langle{\bf y}^{k+1}\hspace{-1mm}-\hspace{-1mm}{\bf y}^*,\nabla g_1(y^{k+1})\hspace{-1mm} -\hspace{-1mm} \nabla g_1({\bf y}^k)\rangle\hspace{-1mm}\geq \Phi \, . \label{eqn:crossx}
\end{align}
To tackle the cross terms related only to ${\utwi{\lambda}}$, we consider the following equality
\begin{align}
\|{\bf a}\hspace{-1mm}-\hspace{-1mm}{\bf c}\|^2_{\bf G}\hspace{-1mm} - \hspace{-1mm}\|{\bf b}\hspace{-1mm}-\hspace{-1mm}{\bf c}\|^2_{\bf G}
\hspace{-1mm}=\hspace{-1mm} 2({\bf a}\hspace{-1mm}-\hspace{-1mm}{\bf c})^T{\bf G}({\bf a}\hspace{-1mm}-\hspace{-1mm}{\bf b}) - \|{\bf a}\hspace{-1mm}-\hspace{-1mm}{\bf b}\|^2_{\bf G},\label{eqn:fact}
\end{align}
and use it in \eqref{eqn:crossx} to arrive at the following inequality:
\begin{align}
&({\bf w}^{k+1} - {\bf w}^*)^T{\bf G}({\bf w}^{k} - {\bf w}^{k+1}) \geq \frac{\gamma}{2}\|{\utwi{\lambda}}^{k+1}-{\utwi{\lambda}}^*\|^2\nonumber\\
-& \frac{\gamma}{2}\|{\utwi{\lambda}}^{k}\hspace{-1mm}-\hspace{-1mm}{\utwi{\lambda}}^{k+1}\|^2\hspace{-1mm} + \hspace{-1mm}\frac{\gamma}{2}\|{\utwi{\lambda}}^{k}\hspace{-1mm}-\hspace{-1mm}{\utwi{\lambda}}^*\|^2
\hspace{-1mm}+\hspace{-1mm}\langle {\mbox{\boldmath $\lambda$}}^{k+1}\hspace{-1mm} - \hspace{-1mm}{\mbox{\boldmath $\lambda$}}^{k}, {\bf A}({\bf x}^{k+1}\hspace{-1mm}-\hspace{-1mm}{\bf x}^{k})\rangle\nonumber\\
+& \beta\gamma\langle {\mbox{\boldmath $\lambda$}}^k\hspace{-1mm}-\hspace{-1mm}{\mbox{\boldmath $\lambda$}}^*\hspace{-1mm}, {\bf A}({\bf x}^{k+1}\hspace{-1mm}-\hspace{-1mm}{\bf x}^{k})\rangle\hspace{-1mm}+\hspace{-1mm}\langle{\bf x}^{*}\hspace{-1mm}-\hspace{-1mm}{\bf x}^{k+1}\hspace{-1mm},\nabla\hspace{-1mm} f_1(x^{k+1})\hspace{-1mm} -\hspace{-1mm} \nabla\hspace{-1mm} f_1({\bf x}^k)\rangle\hspace{-1mm} \nonumber\\
+&\langle{\bf y}^{*}\hspace{-1mm}-\hspace{-1mm}{\bf y}^{k+1},\nabla g_1(y^{k+1})\hspace{-1mm} -\hspace{-1mm} \nabla g_1({\bf y}^k)\rangle\hspace{-1mm}\nonumber\\
+&\langle\beta{\bf B}^T{\bf B}({\bf y}^{k+1} - {\bf y}^k),{\bf y}^{k+1}-{\bf y}^* \rangle+\Phi.\label{eqn:crossy}
\end{align}
Then, we utilize the Cauchy-Schwarz inequality to bound the following term:
\begin{align}
&\langle {\bf A}({\bf x}^{k+1} -{\bf x}^k),{\mbox{\boldmath $\lambda$}}^{k+1}-{\mbox{\boldmath $\lambda$}}^k\rangle \nonumber\\
\geq& -\frac{1}{2\rho_1}\|{\bf A}({\bf x}^{k+1}\hspace{-1mm}-\hspace{-1mm}{\bf x}^k)\|^2 \hspace{-1mm}-\hspace{-1mm}\frac{\rho_1}{2}\|{\mbox{\boldmath $\lambda$}}^{k+1} \hspace{-1mm}-\hspace{-1mm}{\mbox{\boldmath $\lambda$}}^k\|^2, \forall \rho_1 >0\label{cross3}
\end{align}
We implement same process for the rest of cross terms:
\begin{align*}
&\beta\gamma\langle {\bf A}({\bf x}^{k+1} -{\bf x}^k),{\mbox{\boldmath $\lambda$}}^k-{\mbox{\boldmath $\lambda$}}^*\rangle \nonumber\\
\geq& -\frac{\beta\gamma}{2\rho_2}\|{\bf A}({\bf x}^{k+1}\hspace{-1mm}-\hspace{-1mm}{\bf x}^k)\|^2\hspace{-1mm} -\hspace{-1mm}\frac{\beta\gamma\rho_2}{2}\|{\mbox{\boldmath $\lambda$}}^{k} \hspace{-1mm}-\hspace{-1mm} {\mbox{\boldmath $\lambda$}}^*\|^2, \forall \rho_2 >0\label{cross4}\\
&\langle\beta{\bf B}^T{\bf B}({\bf y}^k- {\bf y}^{k+1}),{\bf y}^{k+1}-{\bf y}^* \rangle\nonumber\\
\geq&-\frac{\beta}{\rho_3}\|{\bf B}({\bf y}^{k+1} - {\bf y}^{k})\|^2 - \beta\rho_3\|{\bf B}({\bf y}^{k+1}-{\bf y}^*)\|^2\nonumber\\
\geq& -\frac{\beta\sigma_{\max}^2({\bf B})}{\rho_3}\|{\bf y}^{k+1}\hspace{-1mm} -\hspace{-1mm} {\bf y}^{k}\|^2-\beta\rho_3\sigma_{\max}^2({\bf B})\|{\bf y}^{k+1}\hspace{-1mm}-\hspace{-1mm}{\bf y}^*\|^2 \, ,
\end{align*}
where $\sigma_{\max}({\bf B})$ is the maximal singular value of ${\bf B}$. The remaining terms in~\eqref{eqn:crossy} are related only to gradients of $f_1, g_1$, which can be bounded as follows:
\begin{align}
&\langle{\bf x}^{*}-{\bf x}^{k+1},\nabla f_1(x^{k+1}) - \nabla f_1({\bf x}^k)\rangle
\nonumber\\
&+\langle{\bf y}^{*}-{\bf y}^{k+1},\nabla g_1(y^{k+1}) - \nabla g_1({\bf y}^k)\rangle\nonumber \\
\geq& -\frac{L_f^2}{2\rho_4}\|{\bf x}^{k+1}-{\bf x}^k\|^2 - \frac{\rho_4}{2}\|{\bf x}^{k+1}-{\bf x}^*\|^2 \nonumber\\
&-\frac{L_g^2}{2\rho_5}\|{\bf y}^{k+1}-{\bf y}^k\|^2 - \frac{\rho_5}{2}\|{\bf y}^{k+1}-{\bf y}^*\|^2,
\end{align}
where we have used the Cauchy-Schwarz inequality and we leveraged the Lipschitz continuity of $f_1,g_1$. Also, from \eqref{eqn:fact}, it can be noticed that:
\begin{align*}
&\|{\bf w}^k-{\bf w}^*\|^2_{\bf G} - \|{\bf w}^{k+1}-{\bf w}^*\|^2_{\bf G} \\
= ~&2({\bf w}^k-{\bf w}^*)^TG({\bf w}^k-{\bf w}^{k+1}) - \|{\bf w}^{k} - {\bf w}^{k+1}\|^2_{\bf G} \, .
\end{align*}
It therefore follows that:
\begin{align*}
&\|{\bf w}^k\hspace{-1mm}-\hspace{-1mm}{\bf w}^*\|^2_{\bf G} \hspace{-1mm}-\hspace{-1mm} \|{\bf w}^{k+1}\hspace{-1mm}-\hspace{-1mm}{\bf w}^*\|^2_{\bf G}  \nonumber\\
\geq& \|{\bf w}^k\hspace{-1mm}-\hspace{-1mm}{\bf w}^{k+1}\|^2_{\bf G} + \gamma\|{\utwi{\lambda}}^{k+1}-{\utwi{\lambda}}^*\|^2 \nonumber\\
-& \gamma\|{\utwi{\lambda}}^{k}-{\utwi{\lambda}}^{k+1}\|^2 + \gamma\|{\utwi{\lambda}}^{k}-{\utwi{\lambda}}^*\|^2\nonumber\\
-&\frac{\sigma_{\max}^2({\bf A})}{\rho_1}\|{\bf x}^{k+1}\hspace{-1mm}-\hspace{-1mm}{\bf x}^k\|^2 \hspace{-1mm}-\hspace{-1mm}\rho_1\|{\mbox{\boldmath $\lambda$}}^{k+1} \hspace{-1mm}-\hspace{-1mm}{\mbox{\boldmath $\lambda$}}^k\|^2\nonumber\\
-&\frac{\beta\gamma\sigma_{\max}^2({\bf A})}{\rho_2}\|{\bf x}^{k+1}\hspace{-1mm}-\hspace{-1mm}{\bf x}^k\|^2\hspace{-1mm} -\hspace{-1mm}\beta\gamma\rho_2\|{\mbox{\boldmath $\lambda$}}^{k} \hspace{-1mm}-\hspace{-1mm} {\mbox{\boldmath $\lambda$}}^*\|^2\nonumber\\
-&\frac{L_f^2}{\rho_4}\|{\bf x}^{k+1}-{\bf x}^k\|^2 - \rho_4\|{\bf x}^{k+1}-{\bf x}^*\|^2 \nonumber\\
-&\frac{L_g^2}{\rho_5}\|{\bf y}^{k+1}-{\bf y}^k\|^2 - \rho_5\|{\bf y}^{k+1}-{\bf y}^*\|^2+ \Phi\nonumber \\
-&\frac{\beta\sigma_{\max}^2({\bf B})}{\rho_3}\|{\bf y}^{k+1}\hspace{-1mm} -\hspace{-1mm} {\bf y}^{k}\|^2\hspace{-1mm}-\hspace{-1mm}\beta\rho_3\sigma_{\max}^2({\bf B})\|{\bf y}^{k+1}\hspace{-1mm}-\hspace{-1mm}{\bf y}^*\|^2
\end{align*}
Rearranging the terms in a suitable way, we arrive at the following inequality:
\begin{align}
&\|{\bf w}^k\hspace{-1mm}-\hspace{-1mm}{\bf w}^*\|^2_{\bf G} \hspace{-1mm}-\hspace{-1mm} \|{\bf w}^{k+1}\hspace{-1mm}-\hspace{-1mm}{\bf w}^*\|^2_{\bf G}
\nonumber\\
\geq& (\frac{1}{\alpha_1}\hspace{-1mm}-\hspace{-1mm} \frac{\sigma_{\max}^2({\bf A})}{\rho_1} \hspace{-1mm}-\hspace{-1mm} \frac{\beta\gamma\sigma_{\max}^2({\bf A})}{\rho_2} - \frac{L_f^2}{\rho_4})\|{\bf x}^k-{\bf x}^{k+1}\|^2\nonumber\\
+&(\frac{1}{\alpha_2} - \beta\frac{\sigma_{\max}^2({\bf B})}{\rho_3} - \frac{L_g^2}{\rho_5})\|{\bf y}^k-{\bf y}^{k+1}\|^2\nonumber\\
+&(\frac{1}{\beta}-\gamma-\rho_1)\|{\utwi{\lambda}}^k-{\utwi{\lambda}}^{k+1}\|^2+(2v_f-\rho_4)\|{\bf x}^{k+1}-{\bf x}^*\|^2\nonumber\\
+&(2v_g-\rho_5-\beta\rho_3\sigma_{\max}^2({\bf B}))\|{\bf y}^{k+1}-{\bf y}^*\|^2 \nonumber\\
+&\gamma\|{\utwi{\lambda}}^{k+1}-{\utwi{\lambda}}^*\|^2 + (\gamma-\beta\gamma\rho_2)\|{\utwi{\lambda}}^k-{\utwi{\lambda}}^*\|^2\label{cross2}.
\end{align}
Recall that the goal is to prove the following inequality:
\begin{align}
\|{\bf w}^k-{\bf w}^*\|^2_{\bf G} \geq (1+\delta)\|{\bf w}^{k+1}-{\bf w}^*\|^2_{\bf G},\label{linconv1}
\end{align}
where $\delta>0$ is a constant. For brevity, denote the right-hand-side of \eqref{cross2} as $C$; then it is sufficient to prove that:
\begin{align*}
&C\geq \delta\|{\bf w}^{k+1}\hspace{-1mm}-\hspace{-1mm}{\bf w}^*\|^2_{\bf G} \nonumber\\
&\hspace{-3mm}=\hspace{-1mm}\frac{\delta}{\alpha_1}\|{\bf x}^{k+1}\hspace{-1mm}-\hspace{-1mm}{\bf x}^*\|^2 \hspace{-1mm}+\hspace{-1mm} \frac{\delta}{\alpha_2}\|{\bf y}^{k+1}\hspace{-1mm}-\hspace{-1mm}{\bf y}^*\|^2\hspace{-1mm}+\hspace{-1mm}\frac{\delta}{\beta}\|{\mbox{\boldmath $\lambda$}}^{k+1}\hspace{-1mm}-\hspace{-1mm}{\mbox{\boldmath $\lambda$}}^*\|^2 \,
\end{align*}
which  requires  the following to hold true:
\begin{align*}
&\frac{1}{\alpha_1} - \frac{\max\sigma^2({\bf A})}{\rho_1} - \frac{\beta\gamma\max\sigma^2({\bf A})}{\rho_2} - \frac{L_f^2}{\rho_4} \geq 0\\
&\frac{1}{\alpha_2} - \frac{\beta\max\sigma^2({\bf B})}{\rho_3} - \frac{L_g^2}{\rho_5}\geq 0,~ \frac{1}{\beta}-\gamma-\rho_1 \geq 0 \\
&2v_f-\rho_4-\frac{\delta}{\alpha_1} \geq 0,~ 2v_g - \rho_5-\beta\rho_3\max\sigma^2({\bf B}) -\frac{\delta}{\alpha_2} \geq 0\\
&\gamma-\frac{\delta}{\beta} \geq 0,~ \gamma-\beta\gamma\rho_2\geq 0.
\end{align*}
From the inequalities above, one can notice  that the constant $\alpha_1,\alpha_2$ is closely related to various constants as well as the  singular values of ${\bf A}$ and ${\bf B}$, denoted as $\sigma({\bf A})$ and $\sigma({\bf B})$, respectively. Specifically, this leads to the following conditions for the step sizes:
\begin{align}
&\frac{\delta}{2v_f\hspace{-1mm}-\hspace{-1mm}\rho_4} \leq \alpha_1 \leq \frac{1}{(\frac{1}{\rho_1}\hspace{-1mm}+\hspace{-1mm}\frac{\beta\gamma}{\rho_2})\max\sigma^2({\bf A})\hspace{-1mm} + \hspace{-1mm}\frac{L_f^2}{\rho_4}},  \label{eqn:alpha1}\\
&\frac{\delta}{2v_g\hspace{-1mm}-\hspace{-1mm}\rho_5\hspace{-1mm}-\hspace{-1mm}\beta\rho_3\max\sigma^2({\bf B})}\leq \alpha_2\leq \frac{1}{\frac{\beta\max\sigma^2({\bf B})}{\rho_3} \hspace{-1mm}+\hspace{-1mm} \frac{L_g^2}{\rho_5}} \, . \label{eqn:alpha2}
\end{align}
To ensure that there exists step sizes $\alpha_1,\alpha_2$ that satisfy the condition above, one can choose $\rho_1 = 1, \rho_2 = 1, \rho_3 = \frac{v_g}{2\beta\max\sigma^2({\bf B})}, \rho_4 = v_f, \rho_5 = v_g$, and from Assumption \ref{Assum_Bounded},\ref{Assum_uniform} we know we have all problem dependent parameters uniform bounded.
As for other parameters, one has that:
\begin{align*}
\frac{1}{\beta}-\gamma-1 \geq 0,~ \gamma-\beta\gamma\geq 0\Rightarrow~ \beta\gamma+\beta \leq 1,~ \beta\leq 1.
\end{align*}
This completes the proof.
\end{proof}
Using  Lemma~\ref{lem:linear}, we can now proven the result of  Theorem~\ref{main}.  For notation simplicity, we define $r = \frac{1}{1+\delta}$. Using~\eqref{eqn:lin} and the triangle inequality, we have that
	\begin{align*}
	&\|{\bf w}^{(k)} \hspace{-1mm}-\hspace{-1mm} {\bf w}^{*,(k)}\|_{\bf G} \leq r\|{\bf w}^{(k-1)}\hspace{-1mm} - \hspace{-1mm} {\bf w}^{*,(k-1)} \hspace{-1mm}+\hspace{-1mm} {\bf w}^{*,(k-1)} \hspace{-1mm}-\hspace{-1mm} {\bf w}^{*,(k)}\|_{\bf G} \\
	& \leq r\|{\bf w}^{(k-1)}\hspace{-1mm} - \hspace{-1mm} {\bf w}^{*,(k-1)}\|_{\bf G} +  r\|{\bf w}^{*,(k-1)} - {\bf w}^{*,(k)}\|_{\bf G} \, .
	\end{align*}
We now find a bound on $\|{\bf w}^{*,(k-1)} \hspace{-1mm}-\hspace{-1mm} {\bf w}^{*,(k)}\|_{\bf G}$.
From \eqref{diff}, \eqref{dual_opt} and the fact that we have chosen $\gamma=1$, we know that
	   \begin{align*}
		&{\bf A}^{(k+1)}{\bf x}^{*,(k+1)}\hspace{-1mm}-\hspace{-1mm}{\bf A}^{(k)}{\bf x}^{*,(k)}\hspace{-1mm} +\hspace{-1mm} {\bf B}^{(k+1)}{\bf y}^{*,(k+1)}\hspace{-1mm}-\hspace{-1mm}{\bf B}^{(k)}{\bf y}^{*,(k)}\nonumber\\
		+&{\bf b}^{(k)}-{\bf b}^{(k+1)}+ {\utwi{\lambda}}^{*,(k+1)} -{\utwi{\lambda}}^{*,(k)} = 0.
		\end{align*}
Move ${\utwi{\lambda}}$ terms to the other side of the equation and take norm for both hand sides, it follows:
		\begin{align*}
		&\|{\utwi{\lambda}}^{*,(k+1)} - {\utwi{\lambda}}^{*,(k)}\| \\
        \leq& \|{\bf A}^{(k+1)}({\bf x}^{*,(k+1)} \hspace{-1mm}-\hspace{-1mm} {\bf x}^{*,(k)}) + ({\bf A}^{(k+1)}-{\bf A}^{(k)}){\bf x}^{*,(k)}\| \\
        &+\|{\bf B}^{(k+1)}({\bf y}^{*,(k+1)} \hspace{-1mm}-\hspace{-1mm} {\bf y}^{*,(k)}) + ({\bf B}^{(k+1)}-{\bf B}^{(k)}){\bf y}^{*,(k)}\|\\
        &+\|{\bf b}^{(k)}-{\bf b}^{(k+1)}\|.
        \end{align*}
        From triangle inequality of norm we know
        \begin{align*}
        &\|{\utwi{\lambda}}^{*,(k+1)} - {\utwi{\lambda}}^{*,(k)}\| \\
        \leq& \|{\bf A}^{(k+1)}\|\|{\bf x}^{*,(k+1)} \hspace{-1mm}-\hspace{-1mm} {\bf x}^{*,(k)}\|\hspace{-1mm} +\hspace{-1mm} \|{\bf B}^{(k+1)}\|\|{\bf y}^{*,(k+1)}\hspace{-1mm} -\hspace{-1mm} {\bf y}^{*,(k)}\| \\ &+ \|{\bf b}^{(k)}-{\bf b}^{(k+1)}\|\hspace{-1mm}+\hspace{-1mm} \|{\bf A}^{(k+1)}-{\bf A}^{(k)}\|\|{\bf x}^{*,(k)}\|\\
		&+\|{\bf B}^{(k+1)}-{\bf B}^{(k)}\|\|{\bf y}^{*,(k)}\|
		\end{align*}
		From \eqref{nedic2} we know that
		\begin{align*}
		   \|{\bf x}^{*,(k)}\| &\leq \sigma_1 + \sqrt{\gamma}\max\|{\utwi{\lambda}}^{\text{opt}, (k)}\|\cdot c \\
		    &\leq \sigma_1 + \sqrt{\gamma}\mathcal{M}c:= \mathcal{J}(\sigma_1), \\
		    \|{\bf y}^{*,(k)}\|&\leq \sigma_2 + \sqrt{\gamma}\max\|{\utwi{\lambda}}^{\text{opt}, (k)}\|\cdot c\\
		    &\leq \sigma_2 + \sqrt{\gamma}\mathcal{M}c:= \mathcal{J}(\sigma_2).
		\end{align*}
		Combining with Assumption~\ref{Assum_Bounded} we can reach the following inequality
		\begin{align*}
		&\|{\utwi{\lambda}}^{*,(k+1)} - {\utwi{\lambda}}^{*,(k)}\| \\
		\leq&\tilde{\sigma}_{{\bf A}}\sigma_{\bf x}+\tilde{\sigma}_{{\bf B}}\sigma_{\bf y}+\sigma_{\bf b}+\sigma_{\bf A}\mathcal{x}(\sigma_1) + \sigma_{\bf B}\mathcal{x}(\sigma_2)
		\triangleq \sigma_{{\utwi{\lambda}}},
	\end{align*}
which gives us
\begin{align}
	\|{\bf w}^{*,(k-1)} - {\bf w}^{*,(k)}\|_{\bf G} \leq \psi :=  \sqrt[]{\frac{\sigma_{\bf x}^2}{\alpha_1} +\frac{\sigma_{\bf y}^2}{\alpha_2} + 2\sigma_{{\utwi{\lambda}}}^2}
	\end{align}
and the desired result is obtained.

\subsection{Proof of Corollary \ref{linear}}

By recursively applying the  Theorem~\ref{main}, we have that
\begin{align*}
	\|{\bf w}^{(k)} \hspace{-1mm}-\hspace{-1mm} {\bf w}^{*,(k)}\|_{\bf G}
    & \leq r^k\|{\bf w}^{(0)}-{\bf w}^{*,(0)}\|_{\bf G} \\
    & +  \sum_{i=1}^{(k)}r^{k-i+1}\|{\bf w}^{*,(i-1)}-{\bf w}^{*,(i)}\|_{\bf G}
\end{align*}
where $\|{\bf w}^{*,(k-1)} - {\bf w}^{*,(k)}\|_{\bf G} \leq \psi$.  Taking $k\rightarrow+\infty$, we can derive
	\begin{align*}
	&\lim_{k\rightarrow\infty}\|{\bf w}^{(k)} - {\bf w}^{*,(k)}\|_{\bf G}  \nonumber\\
    \leq&\lim_{k\rightarrow\infty}\left( \frac{r(1-r^{(k)})}{1-r}\psi(\sigma) + r^k\|{\bf w}^{(0)}-{\bf w}^{*,(0)}\|_{\bf G}\right)\\
	\Rightarrow &\lim\sup_{k\rightarrow\infty}\|{\bf w}^{(k)} - {\bf w}^{*,(k)}\|_{\bf G} \leq \frac{r}{1-r}\psi = \frac{\psi}{\delta}.
	\end{align*}
	The desired result is then obtained.

\subsection{Remarks on Corollary 2}
The result is obtained by choosing the biggest step sizes;  that is:
\begin{align*}
    \alpha_1 = \frac{1}{(1+\beta\gamma)\tilde{\sigma^2}_{\bf A} + \frac{\tilde{L}_f^2}{\tilde{v}_f}}, ~\alpha_2 = \frac{1}{\frac{2\beta^2\max\tilde{\sigma}^4_{\bf B}}{\tilde{v}_g}+ \frac{\tilde{L}_g^2}{\tilde{v}_g}}
\end{align*}
Recall that:
\begin{align}
	    \delta &\leq \frac{\tilde{v}_f}{(1+\beta\gamma)\tilde{\sigma}^2_{\bf A} + \frac{\tilde{L}_f^2}{\tilde{v}_f}},~\delta \leq \frac{\tilde{v}_g}{\frac{4\beta^2\tilde{\sigma}^4_{\bf B}}{\tilde{v}_g} + \frac{2\tilde{L}_g^2}{\tilde{v}_g}}.\label{eqn:delta}
\end{align}
We already know that $\delta \leq \beta\gamma$; therefore, based on~\eqref{eqn:alpha1}--\eqref{eqn:alpha2}, one can pick $\delta$ as
\begin{align*}
\delta = \min \left(\frac{\tilde{v}_f}{(1+\beta\gamma)\tilde{\sigma}^2_{\bf A} + \frac{\tilde{L}_f^2}{\tilde{v}_f}},\frac{\tilde{v}_g}{\frac{4\beta^2\tilde{\sigma}^4_{\bf B}}{\tilde{v}_g} + \frac{2\tilde{L}_g^2}{\tilde{v}_g}}, \beta\gamma \right).
\end{align*}.

\ifCLASSOPTIONcaptionsoff
  \newpage
\fi

\bibliographystyle{IEEEtran}

\section{\@startsection {section}{1}{\z@}%
                                   {-1.5ex \@plus -1ex \@minus -.1ex}%
                                    {0.5ex \@plus.1ex}%
                                   {\normalfont\large\bfseries}}
\renewcommand\subsection{\@startsection{subsection}{2}{\z@}%
                                     {-1.0ex\@plus -1ex \@minus -.1ex}%
                                     {0.4ex \@plus .1ex}%
                                     {\normalfont\large\bfseries}}

\section{Introduction}
	\textcolor{black}{This paper considers  time-varying optimization problems for network systems, where objective and constraints evolve over time~\cite{popkov2005gradient,simonetto2014double,Shahrampour2018,hall2015online,bernstein2018online,simonetto2016class}}. The applicability of time-varying optimization problems is evident in a number of domains including power grids~\cite{dall2018optimal,tang2017real,hauswirth2016projected}, communication systems~\cite{Low1999Flow,Chen12}, and online methods in signal processing \cite{asif2014sparse}, just to name a few; see also the representative works~\cite{rahili2017distributed,fazlyab2016self,neely2017online} and~\cite{dall2019optimization} for additional time-varying models and application examples.

	\textcolor{black}{In particular, assume that the temporal domain is discretized as $\{k \tau, k \in \mathbb{N}\}$, with $\tau>0$ being a given interval time~\cite{simonetto2016class,dall2019optimization}. This paper focuses on time-varying optimization problems in the following form~\cite{popkov2005gradient,simonetto2014double}:
		\begin{subequations}
			\label{P1}
			\begin{align}
				\textrm{(P1)~~} \min_{{\bf x}\in\mathbb{R}^m,{\bf y}\in\mathbb{R}^n}~&f^{(k)}({\bf x}) + g^{(k)}({\bf y})\label{P1C}\\
				\text{s.t.}~&{\bf A}^{(k)}{\bf x}+{\bf B}^{(k)}{\bf y}={\bf b}^{(k)}\label{constraint_intro}\\
				&{\bf x}\in\mathcal{X}^{(k)}, {\bf y}\in\mathcal{Y}^{(k)},
			\end{align}
		\end{subequations}
		where $k$ is the index for the time steps $\{k \tau, k \in \mathbb{N}\}$;  $f^{(k)}(\cdot)$ and $g^{(k)}(\cdot)$ are strongly convex functions (for all times $k \tau$); ${\bf A}^{(k)}\in\mathbb{R}^{\ell\times m}, {\bf B}^{(k)}\in\mathbb{R}^{\ell\times n}$ are time-varying matrices; and, $\mathcal{X}^{(k)},\mathcal{Y}^{(k)}$ are convex constraint sets. For optimization and learning problems with streams of data, the interval $\tau$ coincides with the inter-arrival date of data points~\cite{dall2019optimization};} when problem \eqref{P1} is associated with a network, the functions $f^{(k)}(\cdot)$ and $g^{(k)}(\cdot)$ can capture performance objectives that evolve over time, whereas~\eqref{constraint_intro} can capture time-varying physical or logical interactions in the network. Denoting as ${\bf x}^{*,(k)}$ an optimal solution of~\eqref{P1} at time $k \tau$, the optimization model~\eqref{P1} leads to an optimal \emph{trajectory}. The problem addressed in this paper pertains to the development of algorithms that enable \emph{tracking of the optimal trajectory} $\{{\bf x}^{*, (k)}\}_{k \in \mathbb{N}}$. It is worth noticing that, for  problems with inequality constraints, one can always add slack variables to re-write them as equality constraints -- thus fitting  the formulation \eqref{P1}.

	Previous efforts\cite{kar2011gossip,olfati2007consensus} have addressed dynamic (consensus-type) problems as a series of static problems, and have assumed a time-scale separation between algorithms and variability of the problem,  so that convergence is reached for each time. However, this might not be  when the system parameters and the problem inputs change fast, at a time scale that is comparable with the execution of one (or a few) algorithmic steps; \cite{rahili2017distributed,fazlyab2016self, cherukuri2017saddle} successfully approach dynamic problem in continuous time, but only for isolated systems where time-varying exogenous inputs are available at a central processor.

	In recent years, an intensive research has focused on real-time implementation: \cite{rahili2017distributed} presents a control algorithm for real-time multi-agent systems with the ability to track optimal trajectory, however, only the cost function is time-varying; \cite{tang2017real} proposes an online algorithm for optimal power flow problem based on quasi-Newton method. It can be shown that proposed algorithm is able to provide suboptimal solution at a fast timescale. The tracking ability hinges on the accurate estimation of second order information. For the same application, \cite{dall2018optimal} leverages dual subgradient method and system feedback (or measurements) to design a tracking algorithm based on a double smoothing strategy. Regularization terms are added in both primal and dual subproblems to prove Q-linear convergence to a neighborhood of optimal solution for each time instance; \cite{bernstein2018online} further extend double smoothing algorithm to more general settings and provide a regret analysis. \textcolor{black}{The authors in~\cite{onlineSaddle} presents a saddle-point method for networked online convex optimization, and provide a regret analysis for problems with a Lipshitz-continuous function}. Lastly,~\cite{colombino2019online,zheng2019implicit} considered online optimization methods for output regulation problems in dynamical systems. \textcolor{black}{It is also worth mentioning that prediction-correction methods have been utilized to solve time-varying convex optimization problems; see e.g.,~\cite{simonetto2016class,fazlyab2017prediction} and the recent survey paper~\cite{simonetto2020time}.}

	This paper proposes the development of an online algorithm for time-varying convex problems based on the  alternating direction method of multipliers (ADMM) method\cite{boyd2011distributed}. Using a quadratic regularization term, ADMM can allow one to deal with nonsmooth terms,  and it exhibits improved convergence properties relative to dual (sub)gradient methods, especially for problems with ill-conditioned dual functions \cite{deng2016global}\cite{hong2017linear}. The paper present a new algorithm that has following characteristics: i) at each step the primal subproblems are solved via proximal gradient descent -- providing favorable scalability to large-scale problems and accommodating non-smooth objectives; ii) a dual perturbation method is utilized, where the dual variables are suitably perturbed at every iteration to gain in convergence rate. Related works along this line include the following: \cite{ling2014decentralized} leverages ADMM to solve a real-time multi-agent problem. But it differs from the present work because  it  considers only consensus constraints (a special case of our general formulation); \cite{cao2019dynamic} considers a dynamic sharing problem, and convergence to a neighborhood is provided under standard assumptions; however, the constraint is also a special case of our formulation.

	\textcolor{black}{We note that relative to existing online primal-dual methods~\cite{onlineSaddle,simonetto2014double,bernstein2018online}, the proposed method can handle non-differentiable costs. We also note that the choice of ADMM as opposed to, e.g. an Arrow-Hurwicz method, is due to the following two reasons: i) problem \eqref{P1} has a two-block structure (i.e., ${\bf x},{\bf y}$) which can be effectively handled by ADMM; and, ii) the objective function has nonsmooth terms, which would nevertheless require modifications of the Arrow-Hurwicz method.}

	Overall our work has the following contributions: i) An ADMM-based dynamic algorithm is proposed to solve a family of time-varying problems. ii) Tracking ability is rigorously proved under mild assumptions; nonsmooth terms are allowed in the objective and no full row-rank assumption is made for the constraints; see Table \ref{tab:conditions} for detailed comparison between the convergence of a few algorithms.

	\begin{table}
		\begin{center}
			\begin{tabular}{ |p{1.0cm}||p{1.3cm}|p{1.5cm}|p{1.7cm}| p{1.1cm}| }
				\hline
				& {\bf Strong Convexity} & {\bf Lipschitz Continuity} & {\bf Full Row Rank} & {\bf Optimality} \\
				\hline
				\hline
				\multirow{4}{4em}{\bf Classic ADMM} & $f^{(k)}$ & $\nabla f^{(k)}$&${\bf A}^{(k)},({\bf B}^{(k)})^T$ & \multirow{4}{4em}{Optimal solution} \\
				\cline{2-4}&$f^{(k)}, g^{(k)}$ & $ \nabla f^{(k)} $ & ${\bf A}^{(k)}$ &\\
				\cline{2-4}&$f^{(k)}$ & $\nabla f^{(k)}, \nabla g^{(k)}$&$({\bf B}^{(k)})^T$ &  \\
				\cline{2-4}&$f^{(k)}, g^{(k)}$ & $\nabla f^{(k)}, \nabla g^{(k)}$&\xmark &\\
				\hline
				{\bf Proposed Algorithm} &   $f^{(k)}, g^{(k)}$  & \xmark   & \xmark & Perturbed solution [cf.~\eqref{AKKT}] \\
				\hline
			\end{tabular}
		\end{center}
		\caption{Trade off between optimality and conditions for linear convergence \cite{deng2016global}. }\label{tab:conditions}
	\end{table}

	Our previous work \cite{zhang2017dynamic} focuses on ADMM-based online algorithms  to track a solution of a domain-specific linearized optimal power flow (OPF) problem in power grids. In this application domain, this work significantly extends \cite{zhang2017dynamic} in the following ways: 1) We extend the work to a more general OPF formulations; 2) we enable the solution of  OPF problems in a distributed setting, where the power system is divided into areas~\cite{autonomous_grids}; and, 3) a perturbed ADMM is proposed, and convergence is proved under milder conditions, which enables a wider range of constraints and objectives.

	The remainder of paper is organized as follows. Section \RomanNumeralCaps 2 will give the general time-varying problem formulation. Section \RomanNumeralCaps 3 will introduce our dynamic algorithm. Section \RomanNumeralCaps 4 will apply proposed algorithm to two applications, one is in power systems, the other one is route selection. Tracking ability is shown in \RomanNumeralCaps 5 and \RomanNumeralCaps 6 in the form of convergence analysis and simulation, respectively.

	\section{Problem Formulation}

	Consider the time-varying  problem \eqref{P1}\footnote{Throughout this paper, boldface characters denote vectors or matrices; characters with superscript ${(k)}$ denote time varying iterates and parameters; for a given vector ${\bf x}$ and matrix ${\bf G}$, $\|{\bf x}\|^2_{{\bf G}} := {\bf x}^T{\bf G}{\bf x}$; $<{\bf x}, {\bf y}>$ denotes the inner product between the vectors ${\bf x}$ and ${\bf y}$. Given a non-differentiable function $h$, the proximal operator is defined as
		$\text{prox}_{h}({\bf x}) = \arg\min_{\bf z}\|{\bf z}-{\bf x}\|^2 + h({\bf z}).$}.  At time $k$, if problem \eqref{P1} is solved to global optimality, then we say that the {\it perfect tracking} is achieved. However, in many applications~\cite{dall2019optimization} such perfect tracking may not be possible because before the problem at time $k$ is solved, it may have already evolved to a new problem. Specifically, iterative algorithms often involve multiple iterations of computing and communication, and by the time algorithms converge for time $k$, problem parameters such as  ${\bf A}^{(k)},{\bf B}^{(k)}, {\bf b}^{(k)}$ might have already changed. Therefore it is desirable to design algorithms with certain ``tracking ability", which means that the iterates can be continuously steered to stay close to the time-varying optimal solutions.

	Let us reformulate problem \eqref{P1} as follows. First, we rewrite the time-varying constraint sets $\mathcal{X}^{(k)},\mathcal{Y}^{(k)}$ into indicator functions in the objective; and then we separate objective into non-differential functions $f_0^{(k)}({\bf x}), g_0^{(k)}({\bf y})$ and differential functions $f_1^{(k)}({\bf x}), g_1^{(k)}({\bf y})$. At time $k$ we consider the following  time-varying problem:
	\begin{subequations}
		\label{P2}
		\begin{align}
			\textrm{(P2)~~}\min_{{\bf x}\in\mathbb{R}^m,{\bf y}\in\mathbb{R}^n}~&f^{(k)}({\bf x}) + g^{(k)}({\bf y}) \\
			\text{s.t.}~&{\bf A}^{(k)}{\bf x}+{\bf B}^{(k)}{\bf y}={\bf b}^{(k)}\label{P2C}
		\end{align}
	\end{subequations}
	where
	\begin{align}
		& f^{(k)}({\bf x}) \hspace{-1mm}:=\hspace{-1mm} f_0^{(k)}({\bf x})+f_1^{(k)}({\bf x}), g^{(k)}({\bf y})\hspace{-1mm} := \hspace{-1mm} g_0^{(k)}({\bf y})+g_1^{(k)}({\bf y}),\nonumber\\
		& {\bf A}^{(k)} := {\bf A}(t_k),{\bf B}^{(k)} := {\bf B}(t_k),{\bf b}^{(k)} := {\bf b}(t_k),\nonumber\\
		& \mathcal{X}^{(k)}:=\mathcal{X}(t_k),\mathcal{Y}^{(k)}:=\mathcal{Y}(t_k).\nonumber
	\end{align}
	Throughout the paper we will assume that the following assumption holds.\textcolor{black}{
		\begin{assumption}
			\label{Assum1}
			For each time $k$, \eqref{P1} is feasible and it satisfies the Slater's  condition; i.e., there exists a point $(\bar{{\bf x}}^{(k)}, \bar{{\bf y}}^{(k)})$, where $\bar{{\bf x}}^{(k)}\in \mbox{\rm int}(\mathcal{X}^{(k)}),\bar{{\bf y}}^{(k)}\in\mbox{\rm int}(\mathcal{Y}^{(k)})$, where $\mbox{\rm int}(\cdot)$  denotes the interior of the set.
		\end{assumption}
	}
	If Assumption \ref{Assum1} does not hold at a time  $k$, the problem formulation would not be well posed  since there is no solution trajectory to track. We introduce additional assumptions next.
	\begin{assumption}
		\label{Assum_uniform}
		For each time $k$, $f^{(k)}, g^{(k)}$ satisfy
		\vspace{-1mm}
		\begin{align}
			\hspace{-4mm}\langle\partial f^{(k)}({\bf x}_1)\hspace{-1mm}-\hspace{-1mm}\partial f^{(k)}({\bf x}_2), {\bf x}_1\hspace{-1mm}-\hspace{-1mm}{\bf x}_2\rangle&\hspace{-1mm}\geq\tilde{v}_f\|{\bf x}_1\hspace{-1mm}-\hspace{-1mm}{\bf x}_2\|, \forall {\bf x}_1,{\bf x}_2\label{eqn:convexity_f}\\
			\hspace{-3mm}\langle\partial g^{(k)}({\bf y}_1)\hspace{-1mm}-\hspace{-1mm}\partial g^{(k)}({\bf y}_2), {\bf y}_1\hspace{-1mm}-\hspace{-1mm}{\bf y}_2\rangle&\geq\tilde{v}_g\|{\bf y}_1\hspace{-1mm}-\hspace{-1mm}{\bf y}_2\|,\forall {\bf y}_1,{\bf y}_2\label{eqn:convexity_g}
		\end{align}
		where $\tilde{v}_f,\tilde{v}_g$ are uniform lower bounds of strongly convex constants for $f^{(k)}, g^{(k)}$. Functions $f_1^{(k)}, g_1^{(k)}$ have Lipschitz-continuous gradients, i.e.,
		\vspace{-1mm}
		\begin{align}
			\|\nabla f_1^{(k)}({\bf x}_1)\hspace{-1mm} -\hspace{-1mm} \nabla f_1^{(k)}({\bf x}_2)\|&\leq \tilde{L}_f\|{\bf x}_1\hspace{-1mm}-\hspace{-1mm}{\bf x}_2\|,\forall {\bf x}_1,{\bf x}_2 \label{eqn:Lp_f}\\
			\|\nabla g_1^{(k)}({\bf y}_1) \hspace{-1mm}- \hspace{-1mm}\nabla g_1^{(k)}({\bf y}_2)\|&\leq\tilde{L}_g\|{\bf y}_1\hspace{-1mm}-\hspace{-1mm}{\bf y}_2\|,\forall{\bf y}_1,{\bf y}_2\label{eqn:Lp_g}
		\end{align}
		where $\tilde{L}_f,\tilde{L}_g$ are uniform upper bounds of Lipschitz constants for $\nabla f_1^{(k)}, \nabla g_1^{(k)}$.
	\end{assumption}
	\begin{assumption}\label{Assum_slater}
		For each time $k$, the functions $f^{(k)}({\bf x}), g^{(k)}({\bf y})$ are coercive;  i.e.,
		\begin{align*}
			f^{(k)}({\bf x})\rightarrow\infty~\text{as}~\|{\bf x}\|\rightarrow\infty,~g^{(k)}({\bf y})\rightarrow\infty~\text{as}~\|{\bf y}\|\rightarrow\infty \, .
		\end{align*}
	\end{assumption}
	\textcolor{black}{
		Assumption \ref{Assum_slater} will be instrumental to ensure that the iterates are bounded. Since continuous coercive functions' level sets $\{{\bf x}| f({\bf x}) \leq \mu_1, \forall \mu_1\}, \{{\bf y}| g({\bf y}) \leq \mu_2, \forall \mu_2\} $ are always compact, the optimal solutions to problem \eqref{P2}, defined as ${\bf x}^{\text{opt},(k)}, {\bf y}^{\text{opt},(k)}$, are bounded, i.e.,
		\begin{align*}
			&\|{\bf x}^{\text{opt},(k)}\|\leq \sigma_1, ~\|{\bf y}^{\text{opt},(k)}\|\leq \sigma_2
		\end{align*}
		for some positive constants  $\sigma_1, \sigma_2$.}

	\section{Online Proximal-ADMM using Perturbations}

	This section presents an ADMM-based algorithm to track an optimal solution trajectory of the time-varying problem \eqref{P1}. As summarized in Table I, the proposed algorithm exhibits linear convergence guarantees under less stringent  conditions relative to existing ADMM-based methods (even for static problems). In fact, although classic ADMM is conceptually simple and easy to implement, the conditions under which it is convergent is shown to be quite restrictive \cite{deng2016global}. We propose a new algorithm by leveraging the idea of dual perturbation \cite{hajinezhad2019perturbed,koshal2011multiuser} and gradient steps; this will provide a way to demonstrate convergence for a larger family of problems. However, a linear convergence rate at milder conditions comes at the cost of ensuring tracking of an approximate Karush-Kuhn-Tucker (KKT) point~\cite{hajinezhad2019perturbed,koshal2011multiuser,bernstein2018online}.

	Accordingly, we propose to add a small perturbation to the dual variable ${\mbox{\boldmath $\lambda$}}$ in the form of $1-\beta\gamma$, where $\gamma>0$ is the perturbation parameter and $\beta \gamma \in(0,1)$. The perturbed augmented Lagrangian function is then defined as
	\begin{align}
		&\mathcal{L}^{(k)}({\bf x},{\bf y};{\mbox{\boldmath $\lambda$}}) =\mathcal{L}_1^{(k)}({\bf x},{\bf y};{\mbox{\boldmath $\lambda$}})+  f_0^{(k)}({\bf x})+g_0^{(k)}({\bf y}) \label{AL2}
	\end{align}
	where
	\begin{align}
		&\mathcal{L}_1^{(k)}({\bf x},{\bf y};{\mbox{\boldmath $\lambda$}})\hspace{-1mm} = \hspace{-1mm}f^{(k)}_1({\bf x})\hspace{-1mm}+\hspace{-1mm}g^{(k)}_1({\bf y}) +\frac{\beta}{2}\|{\bf A}^{(k)}{\bf x}\hspace{-1mm}+\hspace{-1mm}{\bf B}^{(k)}{\bf y}\hspace{-1mm}-\hspace{-1mm}{\bf b}^{(k)}\|^2\nonumber\\
		&- (1-\beta\gamma){\mbox{\boldmath $\lambda$}}^T\hspace{-1mm}({\bf A}^{(k)}{\bf x}+{\bf B}^{(k)}{\bf y}\hspace{-1mm}-\hspace{-1mm}{\bf b}^{(k)}) \nonumber.
	\end{align}
	Mirroring~\cite{chambolle2011first},  to update ${\bf x}$ and ${\bf y}$, one can performs the following steps in an online fashion (where $k$ is here the time index):
	\begin{align*}
		{\bf y}^{(k+1)} \hspace{-1mm}&=\hspace{-1mm}\arg\min_{{\bf y}}\;\left\langle\frac{\partial\mathcal{L}_1^{(k+1)}({\bf x}^{(k)},{\bf y}^{(k)};{\mbox{\boldmath $\lambda$}}^{(k)})}{\partial{\bf y}^{(k)}},{\bf y}-{\bf y}^{(k)}\right\rangle\\
		&+ g_0^{(k+1)}({\bf y})  + \frac{1}{2\alpha_2}\|{\bf y}-{\bf y}^{(k)}\|^2,\\
		{\bf x}^{(k+1)} \hspace{-1mm}&=\hspace{-1mm}\arg\min_{{\bf x}}\; \left\langle\frac{\partial\mathcal{L}_1^{(k+1)}({\bf x}^{(k)},{\bf y}^{(k+1)};{\mbox{\boldmath $\lambda$}}^{(k)})}{\partial{\bf x}^{(k)}},{\bf x}-{\bf x}^{(k)}\right\rangle\\
		&+ f_0^{(k+1)} ({\bf x}) + \frac{1}{2\alpha_1}\|{\bf x}-{\bf x}^{(k)}\|^2,
	\end{align*}
	where $\alpha_1,\alpha_2$ are step sizes. We are now ready to outline the \emph{online  proximal-ADMM algorithm} with perturbations, whose steps are the following (we stress again that $k$ is  the time index):
	\begin{subequations}
		\label{perturb}
		\begin{align}
			\hspace{-1cm}&\hspace{-2mm}{\bf y}^{(k+1)} \hspace{-1mm}= \hspace{-1mm}\underset{g_0^{(k+1)}}{\textrm{prox}}\hspace{-1mm}\left({\bf y}^{(k)} \hspace{-1mm}-\hspace{-1mm} \alpha_2\frac{\partial\mathcal{L}_1^{(k+1)}({\bf x}^{(k)},{\bf y}^{(k)};{\mbox{\boldmath $\lambda$}}^{(k)})}{\partial {\bf y}^{(k)}}\right),\label{itr:y}\\
			&\hspace{-2mm}{\bf x}^{(k+1)} \hspace{-2mm}= \hspace{-1mm}\underset{f_0^{(k+1)}}{\textrm{prox}}\hspace{-1mm}\left({\bf x}^{(k)} \hspace{-2mm}- \hspace{-1mm}\alpha_1\frac{\partial\mathcal{L}_1^{(k+1)}({\bf x}^{(k)},{\bf y}^{(k+1)};{\mbox{\boldmath $\lambda$}}^{(k)})}{\partial {\bf x}^{(k)}}\hspace{-1mm}\right),\label{itr:x}\\
			&\hspace{-2mm}{\mbox{\boldmath $\lambda$}}\hspace{-0.5mm}^{(k+1)} \hspace{-2mm}= \hspace{-1mm}(1\hspace{-1mm}-\hspace{-1mm}\beta\gamma){\mbox{\boldmath $\lambda$}}\hspace{-0.5mm}^{(k)}\hspace{-1mm} - \hspace{-1mm}\beta\hspace{-1mm}\left({\bf A}\hspace{-0.5mm}^{(k+1)}{\bf x}\hspace{-0.5mm}^{(k+1)}\hspace{-1.5mm}+\hspace{-1mm}{\bf B}\hspace{-0.5mm}^{(k+1)}{\bf y}\hspace{-0.5mm}^{(k+1)}\hspace{-1.5mm}-\hspace{-1mm}{\bf b}\hspace{-0.5mm}^{(k+1)}\hspace{-1mm}\right)\label{perturb_dual}.
		\end{align}
	\end{subequations}
	The underlying assumption here is that for each time $k$, the interval $\tau$ is sufficient  to run at least one iteration of \eqref{perturb}.

	Compared to classical ADMM-based algorithms (for both static and time-varying optimization),  key differences here are in the proximal gradient steps in the primal update and the perturbation added to ${\utwi{\lambda}}$. The proximal gradient steps may provide favorable computational gains when applied to a large-scale problem; it also facilitate ones to develop measurement-based algorithms as discussed in, e.g.,~\cite{dall2019optimization}. \textcolor{black}{The perturbation added to ${\utwi{\lambda}}$ emerges when considering a regularized Lagrangian function of the form $\mathcal{L}^{(k)}({\bf x},{\bf y};{\utwi{\lambda}}) - \frac{\gamma}{2}\|{\utwi{\lambda}}\|^2$ as in e.g.,~\cite{bernstein2018online,koshal2011multiuser,hajinezhad2019perturbed}. This additional term renders the regularized Lagrangian strongly concave in ${\utwi{\lambda}}$. Adding a (small) perturbation in dual variable (or, equivalently, considering a regularized Lagrangian) } is a very useful technique to ensure convergence of the ADMM, and even obtain a linear convergence behavior as explained shortly. To gain   intuition,  let us consider a toy example as follows:
	\begin{align}
		\min_{{\bf x}} \quad 0, \quad \mbox{s.t.} \quad {\bf A} {\bf x} =0
	\end{align}
	where ${\bf A}$ is some fixed matrix, not necessarily positive semidefinite.
	The optimality condition for the above problem can be written down as the following saddle point problem
	\begin{align}
		\min_{{\bf x}}\max_{{\mbox{\boldmath $\lambda$}}}~ {\bf x}^T{\bf A}{\mbox{\boldmath $\lambda$}}\label{perturb_prob}.
	\end{align}
	One can apply the alternating gradient descent/ascent method for solving problem \eqref{perturb_prob}, whose steps are similar as \eqref{perturb} and are given below
	\begin{align*}
		{\bf x}^{(k+1)} &= {\bf x}^{(k)} - \alpha({\bf A}{\mbox{\boldmath $\lambda$}}^{(k)}),\\
		{\mbox{\boldmath $\lambda$}}^{(k+1)} &= {\mbox{\boldmath $\lambda$}}^{(k)} + \beta({\bf A}^T{\bf x}^{(k)}).
	\end{align*}
	In Figure \ref{fig:perturb}, we plot ${\bf x}^T{\bf A}{\utwi{\lambda}}$ using a random matrix ${\bf A}$. An interesting observation  is that the algorithm will diverge if no perturbation is added to ${\bf y}$ as shown in Figure \ref{c1}; also see \cite{daskalakis2018limit} for a formal proof. However, once a small perturbation is added to ${\bf y}$ in both primal and dual updates, i.e.
	\begin{align*}
		{\bf x}^{(k+1)} &= {\bf x}^{(k)} - \alpha({\bf A}{\utwi{\lambda}}^{(k)}(1-\gamma\beta)),\\
		{\mbox{\boldmath $\lambda$}}^{(k+1)} &= {\mbox{\boldmath $\lambda$}}^k(1-\gamma\beta) + \beta({\bf A}^T{\bf x}^{(k)}),
	\end{align*}
	where $\gamma>0$ is a small number, the algorithm will converge as shown in Figure \ref{c2}. This example serves as a motivation to use the perturbation technique.

	\begin{figure}
		\centering
		\begin{subfigure}[b]{0.40\textwidth}
			\centering
			\includegraphics[width=\textwidth]{no_perturb-eps-converted-to.pdf}
			\caption[]%
			{{\small Performance without perturbation}}
			\label{c1}
		\end{subfigure}
		\hfill
		\begin{subfigure}[b]{0.40\textwidth}
			\centering
			\includegraphics[width=\textwidth]{with_perturb-eps-converted-to.pdf}
			\caption[]%
			{{\small Performance with perturbation.}}
			\label{c2}
		\end{subfigure}
		{\small }
		\caption{Example of trends of the objective value of \eqref{perturb_prob} for methods with and without perturbation.}
		\label{fig:perturb}
	\end{figure}

	\section{Convergence Analysis}

	In this section we provide analytical results for convergence and tracking ability of the proposed algorithm.

	From \cite{deng2016global}, it is known that existing ADMM has relatively strict conditions for linear convergence and these conditions may not hold true in some applications; for example, the coefficient matrices in constraints
	\eqref{P2C} might not have full row rank (this is the case for the application presented later in the paper). Further, in some applications, the objective function of \eqref{P2} may also contain non-smooth terms, which can jeopardize the Lipschitz continuity property.
	In contrast, the proposed algorithm could be utilized in a wider range of time-varying optimization problems. It is also worth pointing out that \cite{deng2016global} deals with static optimization problems; here, the focus is on time-varying settings. Next,  we analyze the convergence of the algorithm. To proceed, we concatenate primal and dual optimizer as $\{{\bf w}^*\} = \{{\bf x}^*;{\bf y}^*;{\mbox{\boldmath $\lambda$}}^*\}$ (for static case) as the optimizer of $\max\limits_{{\utwi{\lambda}}}\min\limits_{{\bf x},{\bf y}}\mathcal{L}^{(k)}$ at time $k$. For notation simplicity we neglect superscript $k$ for static case and we have:
	\begin{subequations}\label{AKKT}
		\begin{align}
			{\bf A}^T{\mbox{\boldmath $\lambda$}}^*-\nabla f_1({\bf x}^*)\in\partial f_0({\bf x}^*)\label{x_opt}\\
			{\bf B}^T{\mbox{\boldmath $\lambda$}}^*- \nabla g_1({\bf y}^*) \in \partial g_0({\bf y}^*) \label{y_opt}\\
			{\bf A}{\bf x}^*+{\bf B}{\bf y}^*-{\bf b}+\gamma{\mbox{\boldmath $\lambda$}}^* = 0\label{dual_opt}.
		\end{align}
	\end{subequations}
	Condition \eqref{AKKT} is a perturbed version of KKT conditions, related to approximate KKT (AKKT)\cite{dutta2013approximate,andreani2011sequential}. \textcolor{black}{Basically, optimizer ${\bf w}^*$ is not necessarily the KKT point of original problem \eqref{P1}, but rather an approximate solution. }
	\textcolor{black}{Let $\{{\bf x}^{\text{opt},(k)}, {\bf y}^{\text{opt},(k)}, {\utwi{\lambda}}^{\text{opt},(k)}\}$ be a KKT point of problem \eqref{P1}, which is also the solution to \eqref{AKKT} when $\gamma = 0$. From Slater's condition, we know that any optimal dual solution ${\utwi{\lambda}}^{\text{opt},(k)}$ is bounded \cite[Lemma 1]{nedic2009approximate}. Let ${\bf v}^{\text{opt},(k)} = \{{\bf x}^{\text{opt},(k)}, {\bf y}^{\text{opt},(k)}\}, {\bf v}^{*,(k)} = \{{\bf x}^{*,(k)}, {\bf y}^{*,(k)}\}$; we can then directly use the result in \cite[Proposition 3.1]{koshal2011multiuser} to show that the distance between ${\bf v}^{*,(k)}$ and ${\bf v}^{\text{opt},(k)}$ is bounded. that is,\footnote{The original result is for smooth strongly convex function, but it is easy to check that it still holds true to our problem.},
		\begin{align}\label{nedic}
			\kappa\|{\bf v}^{*,(k)} - {\bf v}^{\text{opt},(k)}\|^2 + \frac{\gamma}{2}\|{\utwi{\lambda}}^{*,(k)}\|^2 \leq \frac{\gamma}{2}\|{\utwi{\lambda}}^{\text{opt},(k)}\|^2,
		\end{align}
		where $\kappa>0$ is a constant. We can further derive that
		\begin{align}\label{nedic2}
			\|{\bf v}^{*,(k)} - {\bf v}^{\text{opt},(k)} \|\leq \sqrt{\gamma}\max\|{\utwi{\lambda}}^{\text{opt},(k)}\|\cdot\textrm{c}, \; \forall~k
		\end{align}
		where $c>0$ is some constant; see \cite[Proposition 3.1]{koshal2011multiuser} for a detailed discussion and the  proof of the result. Since every optimal dual variable ${\utwi{\lambda}}^{\text{opt},(k)}$ to the original problem \eqref{P1} is bounded \cite[Lemma 1]{nedic2009approximate}, it follows that the distance between the AKKT point and a KKT point is bounded too, and such a bound depends on the choice of $\gamma$ (the smaller the $\gamma$, the smaller the distance since $\|{\utwi{\lambda}}^{\text{opt},(k)}\|$ is independent of $\gamma$).
		In the following analysis, we focus on bounding the distance between iterates generated by \eqref{perturb} and AKKT points.
		To proceed, we first state the assumptions on ${\bf x}^{*,(k)}, {\bf y}^{*,(k)}$ and on problem parameters as follows.
		\begin{assumption}
			\label{Assum_Bounded}
			The successive difference between AKKT points is bounded:
			\begin{align}
				\|{\bf x}^{*,(k+1)} - &{\bf x}^{*,(k)}\|\leq \sigma_{{\bf x}}, \,\,\,\, \|{\bf y}^{*,(k+1)} - {\bf y}^{*,(k)}\|\leq \sigma_{{\bf y}},\label{diff}
			\end{align}
			where ${\bf x}^{*,(k)},{\bf y}^{*,(k)}$ are primal optimizer of \eqref{AKKT} at time $k$; $\sigma_{\bf x}>0, \sigma_{\bf y} >0$ are some constants. Also, the variation of the problem parameters is bounded as:
			\begin{align}
				&\|{\bf A}^{(k+1)} - {\bf A}^{(k)}\| \leq \sigma_{\bf A}, ~\|{\bf B}^{(k+1)} - {\bf B}^{(k)}\| \leq \sigma_{\bf B} \label{diff_ab}\\
				&\|{\bf A}^{(k)}\|\leq\tilde{\sigma}_{{\bf A}},~\|{\bf B}^{(k)}\|\leq\tilde{\sigma}_{{\bf B}},~\|{\bf b}^{(k+1)}-{\bf b}^{(k)}\|\leq\sigma_{\bf b}\label{parameter}
			\end{align}
			where $\sigma_{\bf A},\sigma_{\bf B},\tilde{\sigma}_{{\bf A}},\tilde{\sigma}_{{\bf B}},\sigma_{\bf b}$ are some given positive constants.
		\end{assumption}
	}
	Assumption \ref{Assum_Bounded} is  common in time-varying optimization~\cite{dall2019optimization, simonetto2014double, fazlyab2016self,bernstein2018online,Shahrampour2018,hall2015online,madden2020bounds}; worst-case bounds for \eqref{diff} can be obtained assuming that  the  sets $\mathcal{X}^{(k)},\mathcal{Y}^{(k)}$ are compact uniformly in time. Another approach is to measure the  distance based on the \textit{optimal drift}, without assuming a specific bound; see, e.g.,~\cite{ling2014decentralized,cao2019dynamic}. The parameters  $\sigma_{{\bf x}}$ and $\sigma_{{\bf y}}$ quantify the maximum variation of the optimal solutions over two consecutive time steps. Since the paper deals with a tracking problem, conventional wisdom would suggest that better tracking performance can be achieved when  \eqref{P1} is not changing rapidly; this will be confirmed in the convergence results presented later (see also ~\cite{dall2019optimization}).

	The main result of the paper is stated next. For notation simplicity, let $\{{\bf w}^{(k)}\}=\{{\bf x}^{(k)};{\bf y}^{(k)};{\utwi{\lambda}}^{(k)}\}$ be the iterates generated by~\eqref{perturb}, and  let  $\{{\bf w}^{*,(k)}\}=\{{\bf x}^{*,(k)},{\bf y}^{*,(k)},{\utwi{\lambda}}^{*,(k)}\}$ be an optimizer of $\max\limits_{{\utwi{\lambda}}}\min\limits_{{\bf x},{\bf y}}\mathcal{L}^{(k)}$.

	\begin{theorem}\label{main}
		Suppose that Assumptions~\ref{Assum1}--\ref{Assum_Bounded} hold for each time $k$.  Let ${\bf G} =
		\textrm{diag}(\frac{1}{\alpha_1}{\bf I}, \frac{1}{\alpha_2}{\bf I}, \frac{1}{\beta}{\bf I})$ be a  positive definite matrix. Assume that the  step size $\beta$ and the perturbation constant $\gamma$ satisfy: $\beta\gamma + \beta \leq 1, \beta \leq 1$. Finally, assume that the step sizes satisfy the  following:
		\begin{align*}
			0 & < \alpha_1 \leq \left((1+\beta\gamma)\tilde{\sigma}^2_{\bf A} + \frac{\tilde{L}_f^2}{\tilde{v}_f} \right)^{-1}  \nonumber \\
			0 & < \alpha_2 \leq \left(\frac{2\beta^2\max\tilde{\sigma}^4_{\bf B}}{\tilde{v}_g}+ \frac{\tilde{L}_g^2}{\tilde{v}_g} \right)^{-1} \, .
		\end{align*}

		Then, at every time $k$, the tracking error of the algorithm~\eqref{perturb} evolves as:
		\begin{align}
			\|{\bf w}^{(k)} \hspace{-1mm}-\hspace{-1mm} {\bf w}^{*,(k)}\|_{\bf G} & \leq
			\frac{1}{1+\delta} \|{\bf w}^{(k-1)}-{\bf w}^{*,(k-1)}\|_{\bf G} \nonumber \\
			& \,\, + \frac{1}{1+\delta} \left(\frac{\sigma_{\bf x}^2}{\alpha_1} +\frac{\sigma_{\bf y}^2}{\alpha_2} + 2\sigma_{{\utwi{\lambda}}}^2\right)^{\frac{1}{2}}
			\label{eq:linear_conv}
		\end{align}
		where $ \sigma_{{\utwi{\lambda}}} :=\tilde{\sigma}_{{\bf A}}\sigma_{\bf x}+\tilde{\sigma}_{{\bf B}}\sigma_{\bf y}+\sigma_{\bf b}+\sigma_{\bf A}\mathcal{x}(\sigma_1) + \sigma_{\bf B}\mathcal{x}(\sigma_2)$, and  $\delta$ satisfies the condition:
		\begin{align}
			0 < \delta \leq \min\left(\frac{\tilde{v}_f}{(1+\beta\gamma)\tilde{\sigma}^2_{\bf A} + \frac{\tilde{L}_f^2}{\tilde{v}_f}},\frac{\tilde{v}^2_g}{4\beta^2\tilde{\sigma}^4_{\bf B} + 2\tilde{L}_g^2}, \beta\gamma\right).
		\end{align}
	\end{theorem}

	\begin{corollary}
		\label{linear}
		Under the assumptions of Theorem~\ref{main}, one has the following asymptotic behavior for the tracking error:
		\begin{align}
			\label{eq:result_theorem}
			\limsup_{k\rightarrow\infty} \|{\bf w}^{(k)} - {\bf w}^{*,(k)}\|_{\bf G} \leq
			\frac{1}{\delta} \left(\frac{\sigma_{\bf x}^2}{\alpha_1} +\frac{\sigma_{\bf y}^2}{\alpha_2} + 2\sigma_{{\utwi{\lambda}}}^2\right)^{\frac{1}{2}} .
		\end{align}
	\end{corollary}

	The proofs of Theorem~\ref{main} and Corollary~\ref{linear} are provided in the Appendix.  \textcolor{black}{We also provide a comment on the boundedness of dual iterates; specifically, the dual iterates ${\utwi{\lambda}}^{(k)}$ can be bounded as follows:
		\begin{align*}
			\|{\utwi{\lambda}}^{(k)}\| &= \|{\utwi{\lambda}}^{(k)} - {\utwi{\lambda}}^{*,(k)} + {\utwi{\lambda}}^{*,(k)}\| \\&\leq \|{\utwi{\lambda}}^{(k)} - {\utwi{\lambda}}^{*,(k)}\| + \|{\utwi{\lambda}}^{*,(k)}\| \, .
		\end{align*}
		From Corollary 1, we know that $\|{\utwi{\lambda}}^{(k)} - {\utwi{\lambda}}^{*,(k)}\|$ is bounded. Also we know from \cite[Proposition 3.1]{koshal2011multiuser} and \eqref{nedic} that $\|{\utwi{\lambda}}^{*,(k)}\|\leq \|{\utwi{\lambda}}^{\text{opt},(k)}\|$, where ${\utwi{\lambda}}^{\text{opt},(k)}$ is the optimal dual solution (without perturbation), which is bounded because of Slater's condition \cite[Lemma 1]{nedic2009approximate}. Therefore, $\|{\utwi{\lambda}}^{(k)}\|$ is bounded too.}

	The result~\eqref{eq:linear_conv} asserts that the proposed algorithm exhibits linear convergence with a contraction coefficient of $1/(1-\delta)$~\cite{dall2019optimization}; the evolution of the tracking error depends on the temporal variability of the optimal trajectory, which is bounded by the second term on the right-hand-side of~\eqref{eq:linear_conv}. It is worth pointing out that, if the problem~\eqref{P1} is static, then~\eqref{eq:linear_conv} boils down to
	$$
	\|{\bf w}^{(k)} \hspace{-1mm}-\hspace{-1mm} {\bf w}^{*}\|_{\bf G} \leq
	\frac{1}{1+\delta}\|{\bf w}^{(k-1)}-{\bf w}^{*}\|_{\bf G}
	$$
	showing linear convergence of the proximal-ADMM method in batch optimization (${\bf w}^{*}$ is in this case the solution of the static problem).

	The asymptotic result~\eqref{eq:result_theorem} matches existing results in online methods for time-varying optimization~\cite{bernstein2018online,dall2019optimization,dixit2019online}. In particular, the bound depends on the $\delta$ (which affects the contraction coefficient) and the maximum variation of the optimal trajectory over two consecutive time steps, and it shows  how the variation of the problem parameters and optimal solutions can affect the  tracking performance.

	Although~\eqref{eq:result_theorem} asserts that the maximum tracking error is bounded, its tightness is to be investigated on a case-by-case basis (i.e., based on the particular evolution of the solution)\footnote{For example, ~\cite{madden2020bounds} showed that a bound of the form~\eqref{eq:result_theorem} is actually tightly met for online gradient and proximal-gradient descent  for a particular sequence of adversarial cost functions.}.

	The following corollary presented  shows how to maximize $\delta$ (and, hence, how to minimize the worst-case tracking bound).

	\begin{corollary}
		\label{stepsize}
		If the step sizes are selected as
		\begin{align*}
			\alpha_1 =\hspace{-1mm} \left(\hspace{-1mm}(1+\beta\gamma)\tilde{\sigma}^2_{\bf A} + \frac{\tilde{L}_f^2}{\tilde{v}_f} \right)^{-1}\hspace{-2mm},
			\alpha_2 =\hspace{-1mm} \left(\frac{2\beta^2\max\tilde{\sigma}^4_{\bf B}}{\tilde{v}_g}+ \frac{\tilde{L}_g^2}{\tilde{v}_g} \right)^{-1}
		\end{align*}
		then one has that the contraction coefficient $\frac{1}{1+ \delta}$ can be computed using the following expression:
		\begin{align}
			\delta = \min\left(\frac{\tilde{v}_f}{(1+\beta\gamma)\tilde{\sigma}^2_{\bf A} + \frac{\tilde{L}_f^2}{\tilde{v}_f}},\frac{\tilde{v}^2_g}{4\beta^2\tilde{\sigma}^4_{\bf B} + 2\tilde{L}_g^2}, \beta\gamma\right).
		\end{align}
	\end{corollary}

	Notice that, for example, if we specify $\beta=0.5$ and $ \gamma=1$, then $\alpha_1,\alpha_2, \delta$ depend only  on the problem itself; i.e,
	\vspace{-2mm}
	\begin{align}
		\delta &= \min\left(\frac{\tilde{v}_f}{\frac{3}{2}\tilde{\sigma}^2_{\bf A} + \frac{\tilde{L}_f^2}{\tilde{v}_f}},\frac{\tilde{v}^2_g}{\tilde{\sigma}^4_{\bf B} + 2\tilde{L}_g^2}, \frac{1}{2}\right),\label{eqn:delta1}\\
		\alpha_1 &= \frac{1}{\frac{3}{2}\tilde{\sigma}^2_{\bf A} + \frac{\tilde{L}_f^2}{\tilde{v}_f}}, ~\alpha_2 = \frac{1}{\frac{\max\tilde{\sigma}^4_{\bf B}}{2\tilde{v}_g}+ \frac{\tilde{L}_g^2}{\tilde{v}_g}}.\label{eqn:a12}
	\end{align}
	As long as one picks $\delta$ as in \eqref{eqn:delta1}, there exist suitable $\alpha_1,\alpha_2$ to ensure convergence (see \eqref{eqn:alpha1}--\eqref{eqn:alpha2} in the proof).

	\section{Example of Motivating Applications}\label{simulation}

	\subsection{Multi-area power grid optimization}\label{power}

	In this section, we briefly outline an example in power grids. We consider a distribution network featuring distributed energy resources (DERs), and we apply the proposed methodology to drive the DER output powers to the solution of an optimization problem encapsulating voltage constraints and given performance objectives. We demonstrate that the proposed methodology is amenable to settings where the distribution system is partitioned in areas; each area is autonomously controlled, and it ``trades'' power with adjacent areas based on given economic objectives~\cite{autonomous_grids}. In contrast,  previous works in the context of real-time optimal power flow involve centralized algorithms~\cite{tang2017real,hauswirth2016projected} or algorithms with a gather-and-broadcast architecture~\cite{dall2018optimal}.

	Similar to~\cite{autonomous_grids}, consider partitioning a power distribution network into $C$ clusters, and denote as ${\mbox{$\mathcal{C}$}}_{i}$ the set of electrical nodes within cluster $i = 1, \ldots, C$. Two clusters $i$ and $j$ are adjacent if there is at least an electrical node $i$ such that $i \in {\mbox{$\mathcal{C}$}}_{i}$ and $i \in {\mbox{$\mathcal{C}$}}_{j}$. Let ${\cal B}_{i,j} := {\mbox{$\mathcal{C}$}}_{i} \cap {\mbox{$\mathcal{C}$}}_{j}$ be the set of \emph{boundary nodes} connecting cluster $i$ to cluster $j$, and define ${\cal B}_{i} := \cup_{j\neq i} {\cal B}_{i,j}$. Further, let ${\cal I}_{i} := {\mbox{$\mathcal{C}$}}_{i} \backslash {\cal B}_{i}$ be the set of \emph{internal nodes} for cluster $i$. For future developments, let $N_{i} := |{\cal I}_{i}|$ be the number of internal nodes if cluster $i$, and let ${\cal N}_{i} \subset \{1, \ldots, C\}$ be the set of neighboring clusters of the $i$th one (i.e., cluster connected to the $i$th one).

	Let ${\bf x}_{j}^{i} := [P_{j}^{i}, Q_{j}^{i}]^\textsf{T} \in \mathbb{R}^2$ collect the net active and reactive powers injected by DERs at the node $j \in {\cal I}_{i}$ of cluster $i$. Particularly, ${\bf x}_{j}^{i} $ can represent the powers injected by \emph{one} DER located at node $j$, or the aggregate net power injections of \emph{a group of} DERs located at node $j$ (e.g., a household with multiple controllable devices) and we stack the setpoints $\{{\bf x}_{j}^{i}\}_{ j \in {\cal I}_{i}}$ in the vector ${\bf x}^{i} \in \mathbb{R}^{2 N_{i}}$. If no controllable DERs are present at a given location, the corresponding vector  ${\bf x}_{j}^{i}$ is  set to $\mathbf{0}$\footnote{For notation simplicity, the model is outlined for balanced systems and for the case where one household/building with DERs is located at a node. However, the model can be trivially extended to multiphase networks~\cite{linModels} and for the case where multiple households/buildings with DERs are located at a node (at the cost of increasing the complexity of the notation).}. On the other hand, ${\utwi{\ell}}_{j}^{i} \in \mathbb{R}^2$ denotes the net non-controllable loads at node $j \in {\cal I}_{i}$, and ${\utwi{\ell}}^{i} \in \mathbb{R}^{2 N_{i}}$ stacks the loads $\{{\utwi{\ell}}_{j}^{i}\}_{ j \in {\cal I}_{i}}$. It is assumed that no DERs and no non-controllable loads are located at the boundary nodes ${\cal B}_{i,j}$.

	Let $V_{j}^{i} \in \mathbb{C}$ denote the complex line-to-ground voltage phasor at node $j$ of cluster $i$, and let ${\bf v}^{i} := [\{|V_{j}^{i}|, j \in {\cal I}_{i}\}]^\textsf{T}$ be the vector of voltage  magnitudes of the internal nodes ${\cal I}_{i}$. For each pair of neighboring clusters $(i,j)$, let ${\bf x}_{n}^{j \rightarrow i} := [P_{n}^{j \rightarrow i}, Q_{n}^{j \rightarrow i}]^\textsf{T} \in \mathbb{R}^2$ represent the active and reactive powers flowing into area $i$ from area $j$ through node $n \in {\cal B}_{i,j} $; on the other hand, ${\bf x}_{n}^{i \rightarrow j} \in \mathbb{R}^2$ contains the active and reactive powers flowing into area $j$ from area $i$ through node $n \in {\cal B}_{i,j} $. From Kirchhoff's Law, it holds that ${\bf x}_{n}^{j \rightarrow i} + {\bf x}_{n}^{i \rightarrow j} = \mathbf{0}$. To facilitate the syntheses of computationally-affordable algorithms, we leverage the following approximate linear relationship between net injected power and voltage magnitude (see e.g.,~\cite{dall2018optimal} and references therein):
	\begin{subequations} \label{relation1}
		\begin{align}
			\tilde{{\bf v}}^{i}  &:=  \hspace{-2mm} \sum_{j \in {\cal I}_{i}} {\bf A}_{j}^{i} ({\bf x}^{i}_j - {\utwi{\ell}}^{i}_j) +\hspace{-2mm}  \sum_{j \in {\cal N}_{i}} \sum_{n \in {\cal B}_{i,j}}\hspace{-2mm} {\bf A}_{n}^{j \rightarrow i} {\bf x}_{n}^{j \rightarrow i} \hspace{-1mm}+\hspace{-1mm} {\bf a}, \label{eqn:lin_v} \\
			&=  {\bf A}^{i} ({\bf x}^{i} - {\utwi{\ell}}^{i}) +  \sum_{j \in {\cal N}^{i}}  {\bf A}^{j \rightarrow i} {\bf x}^{j \rightarrow i} + {\bf a}, \label{eqn:lin_v2}
		\end{align}
	\end{subequations}
	\textcolor{black}{where ${\bf A}^i = [{\bf A}_j^i]_{j\in{\cal I}_i}, {\bf A}^{j\rightarrow i} = [{\bf A}^{j\rightarrow i}_n]_{n\in{\cal B}_{i,j}}, {\bf a}$ are time-varying problem parameters derived from linearized power flow equation. }
	Another linear relationship between net injected power and power between clusters is captured in the following equation:
	\begin{subequations}\label{relation2}
		\begin{align}
			{\bf x}^{j \rightarrow i} &:=  \sum_{k \in {\cal I}_{i}} {\bf M}^{j \rightarrow i}_{k}  ({\bf x}^{i}_k - {\utwi{\ell}}^{i}_k)  + \mathbf{m}^{j \rightarrow i}, \label{eqn:lin_s0} \nonumber\\  + & \sum_{k \in {\cal N}_{i} \backslash \{j\}} \sum_{n \in {\cal B}_{i,k}} {\bf M}^{k, j \rightarrow i}_{n}{\bf x}_{n}^{k \rightarrow i}\\
			&=   {\bf M}^{j \rightarrow i}  ({\bf x}^{i}\hspace{-1mm} - \hspace{-1mm}{\utwi{\ell}}^{i}) + \mathbf{m}^{j \rightarrow i}+ \hspace{-4mm} \sum_{k \in {\cal N}_{i} \backslash \{j\}} \hspace{-3mm}{\bf M}^{k, j \rightarrow i} {\bf x}^{j \rightarrow i}, \label{eqn:lin_s01}
		\end{align}
	\end{subequations}
	where ${\bf M}^{j\rightarrow i} \hspace{-1mm}=\hspace{-1mm} [{\bf M}^{j\rightarrow i}_k]_{k\in{\cal I}_i}, {\bf M}^{k,j\rightarrow i} \hspace{-1mm}=\hspace{-1mm} [{\bf M}^{k,j\rightarrow i}_{n}]_{n\in{\cal B}_{i,k}}$,${\bf m}^{j\rightarrow i}$ are also time-varying problem parameters depending on the actual network physics.
	All model parameters in \eqref{relation1}--\eqref{relation2} can be obtained as shown in~\cite{linModels}.
	Now we are ready to state our real-time OPF problem as follows:
	\begin{subequations}
		\begin{align}
			&\hspace{-5mm}\min_{\substack{\{{\bf x}^{i}\}, \{{\bf x}_{n}^{j \rightarrow i}, {\bf x}_{n}^{i \rightarrow j}\}}}  \hspace{.2cm} \sum_{i=1}^C  [f^{i}({\bf x}^{i}) + g^i(\{{\bf x}^{j \rightarrow i}\})] \tag{P3} \label{obj} \\
			\text{s.t.} ~&{\bf x}_{j}^{i} \in {\cal Y}_{j}^{i}, \forall \,\, j \in {\cal I}^{i}, \,\, i = 1, \ldots C \label{eqn:constr_Y_exact}  \\
			& v^{\text{min}} \mathbf{1} \leq \tilde{{\bf v}}^{i} \leq v^{\text{max}} \mathbf{1}, \forall \,\, i = 1, \ldots C \label{cons}\\
			& {\bf x}^{j \rightarrow i} =   {\bf M}^{j \rightarrow i}  ({\bf x}^{i} - {\utwi{\ell}}^{i}) + \mathbf{m}^{j \rightarrow i}+\hspace{-5mm}  \sum_{k \in {\cal N}_{i} \backslash \{j\}} \hspace{-5mm}{\bf M}^{k, j \rightarrow i} {\bf x}^{j \rightarrow i} \nonumber\\
			&\hspace{3cm}, \, \forall \, j \in {\cal N}_{i}, \, \, i = 1, \ldots, C \label{between}\\
			& {\bf x}^{j \rightarrow i} + {\bf x}^{i \rightarrow j} = \mathbf{0} , \,\, \forall \,\, \textrm{neighboring~areas~} (i,j)
		\end{align}
	\end{subequations}
	where the time-varying objective function models the amount of real power curtailed and the amount of reactive power injected or absorbed (which leads to non-smooth term in the objective, e.g. $\ell_1$ term). For notation simplicity, we write objective function in \eqref{obj} as $\Psi({\bf x})$. Consider ${\bf M}^{j\rightarrow i}$ consists of $1,0$, with $1$ for real power, $0$ for reactive power. Putting \eqref{eqn:lin_s01} back to \eqref{eqn:lin_v2}, adding slack variables ${\utwi{\gamma}}^i,{\utwi{\beta}}^i$ to \eqref{cons} formulate equality constraints, and adding strongly convex term w.r.t ${\utwi{\gamma}} = \{{\utwi{\gamma}}^i\},{\utwi{\beta}} = \{{\utwi{\beta}}^i\}$ we have the following formulation:
	\begin{subequations}
		\begin{align}
			&\hspace{-0.5cm}\min_{\substack{\{\{{\bf x}^{i}\}, \{{\bf x}_{n}^{j \rightarrow i}, {\bf x}_{n}^{i \rightarrow j}\}\},\{{\utwi{\gamma}}^{i},{\utwi{\beta}}^{i}\geq 0\}\} } } \Psi({\bf x})+ a\|{\utwi{\gamma}}\|^2 + b\|{\utwi{\beta}}\|^2\tag{P4} \label{obj2} \\
			\text{s.t.~} & {\bf x}_{j}^{i} \in {\cal Y}_{j}^{i},  \forall \,\, j \in {\cal I}_{i}, \,\, i = 1, \ldots C \label{eqn:constr_Y_exact}  \\
			& v^{\text{min}}\mathbf{1} -\tilde{{\bf v}}^{i} + {\utwi{\gamma}}^{i}= \mathbf{0},  \forall \,\, i = 1, \ldots C \label{cons2}\\
			&\tilde{{\bf v}}^{i} + {\utwi{\beta}}^{i} - v^{\text{max}}\mathbf{1}=\mathbf{0} ,  \forall \,\, i = 1, \ldots C \label{cons3}\\
			& {\bf x}^{j \rightarrow i} =   {\bf M}^{j \rightarrow i}  ({\bf x}^{i} - {\utwi{\ell}}^{i}) + \mathbf{m}^{j \rightarrow i}+\hspace{-2mm}  \sum_{k \in {\cal N}_{i} \backslash \{j\}}\hspace{-2mm} {\bf M}^{k, j \rightarrow i} {\bf x}^{j \rightarrow i} \nonumber\\
			&\hspace{3cm}, \, \forall \, j \in {\cal N}_{i}, \, \, i = 1, \ldots, C \label{eqn:between}\\
			& {\bf x}^{j \rightarrow i} + {\bf x}^{i \rightarrow j} = \mathbf{0} , \,\, \forall \,\, \textrm{neighboring~areas~} (i,j).\label{eqn:consensus}
		\end{align}
	\end{subequations}
	We can now clearly see a mapping from \eqref{obj2} to \eqref{P1}: objective functions
	are $\Psi({\bf x})$ and $a\|{\utwi{\gamma}}\|^2 + b\|{\utwi{\beta}}\|^2$(where $a,b > 0$ are small); two blocks of variables are $\{{\bf x}^{i},{\bf x}^{j\rightarrow i}\}$ and $\{{\utwi{\gamma}}^{i},{\utwi{\beta}}^{i}\}$; constraints are all linear and separable w.r.t each network node. \textcolor{black}{Problem \eqref{obj2} is time varying in both objective function and constraint parameters.} In order to better illustrate how the proposed algorithm can be applied, we use a 4-cluster network (see Figure \ref{fig:net}) as an example. First, we substitute ${\bf x}^{j\rightarrow i}$ in \eqref{eqn:lin_v2} with \eqref{eqn:between}; then, we substitute $\tilde{{\bf v}}^i$ in \eqref{cons2}--\eqref{cons3} with \eqref{eqn:lin_v2};
	last, we define the corresponding augmented Lagrangian function as follows:
	\vspace{-2mm}
	\begin{align*}
		&\mathcal{L}({\bf x},{\utwi{\gamma}},{\utwi{\beta}},{\utwi{\lambda}}) = \Psi({\bf x}) + a\|{\utwi{\gamma}}\|^2 + b\|{\utwi{\beta}}\|^2 \nonumber\\
		+ &\sum_{k}\hspace{-1mm}\sum_{i\in C_k}\frac{\rho}{2}\|\textcolor{black}{\sum_{j\in C_k}({\bf A}^{i}_j+{\bf A}_j^{3}){\bf x}^{i}_j  \hspace{-1mm}+\hspace{-1mm} {\bf a}_k }\hspace{-1mm}+\hspace{-1mm}{\utwi{\beta}}^{i} \hspace{-1mm}-\hspace{-1mm} v^{\text{max}} \hspace{-1mm}+ \hspace{-1mm}\frac{{\utwi{\lambda}}_1(1-\gamma)}{\rho}\|^2\\
		+&\sum_{k}\hspace{-1mm}\sum_{i\in C_k} \frac{\rho}{2}\|v^{\text{min}} \hspace{-1mm}\textcolor{black}{-\hspace{-1mm}\sum_{j\in C_k} ({\bf A}^{i}_j\hspace{-1mm}+\hspace{-1mm}{\bf A}^{3}_j){\bf x}^{i}_j \hspace{-1mm}- \hspace{-1mm}{\bf a}_k} \hspace{-1mm}+\hspace{-1mm} {\utwi{\gamma}}^{i}\hspace{-1mm}+\hspace{-1mm}\frac{{\utwi{\lambda}}_2(1-\gamma)}{\rho}\|^2\\
		+& \sum_{k}\frac{\rho}{2}\|\sum_{i}{\bf x}^{i\rightarrow j} - \textcolor{black}{\sum_{j\in C_k} {\bf x}^{j}} + \frac{{\utwi{\lambda}}_4(1-\gamma)}{\rho}\|^2\\
		+&\frac{\rho}{2}\|\sum_{i\neq j}{\bf x}^{i\rightarrow j}+\frac{{\utwi{\lambda}}_3(1-\gamma)}{\rho}\|^2
	\end{align*}

	The detailed updates follow the same way as \eqref{perturb} and from Table I we know that linear convergence to AKKT is guaranteed. To further improve our algorithm for this particular application. We incorporate system measurements in both primal and dual updates in the following way:
	\begin{align}
		\sum_{j\in C_k}({\bf A}^{i}_j+{\bf A}_j^{3}){\bf x}^{i}_j +{\bf a}_k &\rightarrow \phi({\bf x}),~ \sum_{j\in C_k}{\bf x}^{j}\rightarrow \psi({\bf x}),
	\end{align}
	where $\phi({\bf x}), \psi({\bf x})$ are measurements. This is beneficial in that: i) A natural distributed computing scheme is achieved while without feedback it is not clear whether the algorithm can be implemented in a distributed way; ii) feedback terms are much less than uncontrollable terms, which essentially shrinks the measuring time; iii) it is easier to satisfy power flow equations with the help of system measurements.

	\begin{figure}
		\centering
		\includegraphics[width=8cm]{net}
		\caption{power distribution network with 4 clusters, node 3 is a boundary node that belongs to all 4 clusters.}
		\label{fig:net}
	\end{figure}

	\subsection{Simulations}
	In this section, we test our algorithm using the same power systems settings. We consider a similar system as in \cite{dall2018optimal}, where a modified IEEE 37-node test feeder is utilized. The network is obtained by considering a single phase equivalent, and by replacing the loads on phase ``c" specified in the original dataset with real load data measured from feeders in a neighborhood called Anatolia in California during a week in August 2012. It is assumed that the aggregations of photovoltaic systems are located at nodes 4, 7, 10, 13, 17, 20, 22, 23, 26, 28, 29, 30, 31, 32, 33, 34, 35, and 36. The rating of  these inverters are 300kVA for $i=3$, 350kVA for $i=15, 16$ and 200kVA for the remaining ones. The objective is set to be $f^i({\bf x}^i) = c_p(P_{\text{av},i} - P_i)^2 + c_q(Q_i)^2 + \bar{c}_q|Q_i|, g^i({\bf x}^{j\rightarrow i}) = 0$ where $P_{\text{av},i}$ is the maximum real power available from the PV system $i$, and $c_p = 3, c_q = 1, \bar{c}_q = 0.1$. The voltage limits are set to be $V^{\text{min}} = 0.95$pu, $V^{\text{max}} = 1.05$pu. The generation profiles are simulated based on real solar irradiance data and have a granularity of 1 second.
	\begin{figure}
		\centering
		\includegraphics[width=8cm]{a2_p-eps-converted-to.pdf}
		\caption{Real power at feeder head during 12:00-12:30.}
		\vspace{-.2cm}
		\label{f2}
	\end{figure}
	First we specify a given trajectory for the power at the common coupling, which is color-coded in red in Fig. \ref{f2} (negative power indicates reverse power flows). It can be seen that our algorithm is able to regulate $P_0^k$ close to $P_{0,\text{set}}^k$ in real time. Figure \ref{f3} illustrates the voltage profiles for selected nodes. From 10:00 to 12:00 we observe a few flickers, which is caused by rapid variations of the solar irradiance. Other than that, it can be seen that voltage regulation is enforced and a flat voltage profile is obtained. Note that even there are some relatively large jumps from around 12:00 to 14:00, our algorithm is still able to track the optimal trajectory. A comparison with double smoothing algorithm \cite{dall2018optimal} is presented in Figure \ref{f4}. The proposed strategy has potentially better voltage regulation ability, especially for extreme cases e.g. the two spikes from 10:00 to 12:00.
	\begin{figure}
		\centering
		\includegraphics[width=8cm]{a2_v-eps-converted-to.pdf}
		\caption{Voltage  profile achieved (only some nodes are considered for illustration purposes). }
		\vspace{-.2cm}
		\label{f3}
	\end{figure}
	\begin{figure}
		\centering
		\includegraphics[width=8cm]{voltage_sum_measurements-eps-converted-to.pdf}
		\caption{Index for the overall voltage violation across the system $\sum\limits_{n\in\mathcal{N}}\left(\max(|V_n^k| - v^{\text{max}} ,0)+\max(v^{\text{min}} - |V_n^k|,0)\right)$}
		\vspace{-.2cm}
		\label{f4}
	\end{figure}
	\begin{figure}[t]
		\centering
		\includegraphics[width=8cm]{cons_0916_c1-eps-converted-to.pdf}
		\caption{Power violation of Cluster 1: power violation for each cluster is defined as $\|\sum\limits_{i\in\mathcal{N}^{(j)}}{\bf x}^{i\rightarrow j} - \sum\limits_{j\in C}{\bf x}^{j}\|^2$}
		\label{power_vio}
	\end{figure}%
	\begin{figure}[t]
		\centering
		\includegraphics[width=8cm]{cons_0916_cv-eps-converted-to.pdf}
		\caption{Consensus violation: $\|{\bf x}^{(i\rightarrow j)}+{\bf x}^{(j\rightarrow i)}\|^2$}
		\label{cons_power}
	\end{figure}

	We proceed to test in the same setting except we are adding consensus constraints. In Figure \ref{power_vio} we can see that for all 4 clusters, power violation decreases dramatically in first a few minutes and remains at a  low level of $10^{-10}$. The power consensus violation is shown in Figure \ref{cons_power}, where a steep drop at the begging and flat low line after that are observed.

	\section{Conclusion}
	This paper gives a general online optimization problem formulation and proposes a dynamic algorithm based on alternating direction method of multipliers that can continuously track optimal solution in real time. The steps of ADMM are proximal gradient steps with modification of adding perturbation to dual variable and incorporating system feedback for certain applications. The resulting algorithm is proved to converge to a neighborhood of optimal solution for each time instance. Numerical results for power systems applications also demonstrate the practicality of the proposed algorithm. Our future research will focus on general online nonconvex optimization problems.

	\appendices

	\subsection{Proof of Theorem \ref{main}}

	To show the result of Theorem~\ref{main}, we start from the following lemma.

	\begin{lemma}
		\label{lem:linear}
		Under the assumptions of Theorem~\ref{main}, it holds that:
		\begin{align}
			\label{eqn:lin}
			\|{\bf w}^{(k)}-{\bf w}^{*,(k)}\|^2_{\bf G} \geq (1+\delta)\|{\bf w}^{(k+1)}-{\bf w}^{*,(k)}\|^2_{\bf G} .
		\end{align}

	\end{lemma}

	\begin{proof}
		Since the lemma focuses on a particular time instant $k$, we replace ${\bf w}^{*,(k)}$ with ${\bf w}^{*}$ for notation simplicity. Also, to better fit the equations in the columns, we replace the subscript $^{(k)}$ and $^{(k+1)}$ with $^{k}$ and $^{k+1}$, respectively, in the proof.

		\textcolor{black}{From the optimality condition of subproblem \eqref{itr:y}, one has that:
			\begin{align}
				&{\bf A}^T{\mbox{\boldmath $\lambda$}}^k(1-\beta\gamma) - \beta {\bf A}^T({\bf A}{\bf x}^{k} + {\bf B}{\bf y}^{k+1}-{\bf b})\nonumber\\-&\nabla f_1({\bf x}^k) + \frac{1}{\alpha_1}{\bf I}({\bf x}^k\hspace{-1mm}-{\bf x}^{k+1})\in\partial f_0({\bf x}^{k+1}).\nonumber
			\end{align}
			Rearrange terms we get
			\begin{align}
				& {\bf A}^{T}{\mbox{\boldmath $\lambda$}}^{k+1}\hspace{-1mm}+\hspace{-1mm}\beta{\bf A}^T{\bf A}({\bf x}^{k+1}\hspace{-1mm} - \hspace{-1mm}{\bf x}^{k})\hspace{-1mm} +\hspace{-1mm} \nabla f_1({\bf x}^{k+1})\hspace{-1mm} -\hspace{-1mm} \nabla f_1({\bf x}^k) \nonumber\\
				+& \frac{1}{\alpha_1}{\bf I}({\bf x}^k-{\bf x}^{k+1})\in\partial f_0({\bf x}^{k+1}) + \nabla f_1({\bf x}^{k+1}).\label{x_k}
			\end{align}
			Also from optimality condition of subproblem and \eqref{itr:x}, one has that
			\begin{align*}
				&{\bf B}^T{\mbox{\boldmath $\lambda$}}^k(1\hspace{-1mm}-\hspace{-1mm}\beta\gamma)\hspace{-1mm} - \hspace{-1mm} \beta {\bf B}^{T}({\bf A}{\bf x}^{k+1}\hspace{-1mm}+\hspace{-1mm} {\bf B}{\bf y}^{k+1}\hspace{-1mm} -\hspace{-1mm} {\bf b})\hspace{-1mm}\nonumber\\
				+&\beta{\bf B}^T{\bf A}({\bf x}^{k+1}-{\bf x}^k) + \beta{\bf B}^T{\bf B}({\bf y}^{k+1}-{\bf y}^k)\nonumber\\
				+& \frac{1}{\alpha_2}{\bf I}({\bf y}^k\hspace{-1mm} -\hspace{-1mm} {\bf y}^{k+1}) -\nabla g_1({\bf y}^k)\in\partial g_0({\bf y}^{k+1}) .\nonumber
			\end{align*}
			Rearrange terms we get
			\begin{align}
				& {\bf B}^T({\mbox{\boldmath $\lambda$}}^{k+1} + \beta {\bf A}({\bf x}^{k+1}-{\bf x}^k) +\beta{\bf B}({\bf y}^{k+1}-{\bf y}^k))\nonumber\\
				&+ \nabla g_1({\bf y}^{k+1}) - \nabla g_1({\bf y}^k)+ \frac{1}{\alpha_2}{\bf I}({\bf y}^k-{\bf y}^{k+1}) \nonumber\\
				\in &\partial g_0({\bf y}^{k+1}) + \nabla g_1({\bf y}^{k+1})\label{y_k}.
			\end{align}
			Furthermore, from the dual update \eqref{perturb_dual}, one can obtain:
			\begin{align}
				\frac{1}{\beta}({\mbox{\boldmath $\lambda$}}^k - {\mbox{\boldmath $\lambda$}}^{k+1}) = \gamma{\mbox{\boldmath $\lambda$}}^k + ({\bf A}{\bf x}^{k+1} + {\bf B}{\bf y}^{k+1} - {\bf b}),\label{dual}
			\end{align}
		}
		and, together with optimality condition, one obtains:
		\begin{align}\label{dual_opt2}
			\frac{1}{\beta}({\mbox{\boldmath $\lambda$}}^k\hspace{-1mm} - \hspace{-1mm}{\mbox{\boldmath $\lambda$}}^{k+1})\hspace{-1mm} =\hspace{-1mm} \gamma({\mbox{\boldmath $\lambda$}}^k\hspace{-1mm}-\hspace{-1mm}{\mbox{\boldmath $\lambda$}}^*) \hspace{-1mm}+\hspace{-1mm} {\bf A}({\bf x}^{k+1}\hspace{-2mm}-\hspace{-1mm}{\bf x}^{*})\hspace{-1mm} + \hspace{-1mm}{\bf B}({\bf y}^{k+1}\hspace{-2mm}-\hspace{-1mm}{\bf y}^{*})
		\end{align}
		Since the functions $f = f_0 + f_1$ and $g = g_0+g_1$ are strongly convex, we leverage \eqref{eqn:convexity_f} and \eqref{eqn:convexity_g}
		and, by plugging the optimality condition \eqref{x_k} and \eqref{y_k}, we have
		\begin{align}
			&\langle{\bf A}^{T}({\mbox{\boldmath $\lambda$}}^{k+1}\hspace{-1mm}-\hspace{-1mm}{\utwi{\lambda}}^{*})+\beta{\bf A}^T{\bf A}({\bf x}^{k+1}\hspace{-1mm} -\hspace{-1mm} {\bf x}^{k})\hspace{-1mm} +\hspace{-1mm} \nabla f_1({\bf x}^{k+1}) \hspace{-1mm}- \hspace{-1mm}\nabla f_1({\bf x}^k)\nonumber\\
			&+ \frac{1}{\alpha_1}{\bf I}({\bf x}^k-{\bf x}^{k+1}),{\bf x}^{k+1}-{\bf x}^*\rangle\geq v_f\|{\bf x}^{k+1}-{\bf x}^*\|^2\label{x}\\
			&\langle{\bf B}^T({\mbox{\boldmath $\lambda$}}^{k+1}-{\utwi{\lambda}}^* + \beta {\bf A}({\bf x}^{k+1}-{\bf x}^k) +\beta{\bf B}({\bf y}^{k+1}-{\bf y}^k)) \nonumber\\
			&+ \nabla g_1({\bf y}^{k+1}) - \nabla g_1({\bf y}^k)+\frac{1}{\alpha_2}{\bf I}({\bf y}^k-{\bf y}^{k+1}),{\bf y}^{k+1}-{\bf y}^*\rangle \nonumber\\
			&\geq v_g\|{\bf y}^{k+1}-{\bf y}^*\|^2\label{y}
		\end{align}
		Next, add \eqref{x} and \eqref{y} together, define $\Phi = v_f\|{\bf x}^{k+1}-{\bf x}^*\|^2+v_g\|{\bf y}^{k+1}-{\bf y}^*\|^2$ and plug in \eqref{dual_opt2} to obtain:

			\begin{align}
				&\langle\beta{\bf A}({\bf x}^{k+1}\hspace{-1mm}-\hspace{-1mm}{\bf x}^k)\hspace{-1mm}+\hspace{-1mm}{\utwi{\lambda}}^{k+1}\hspace{-1mm}-\hspace{-1mm}{\utwi{\lambda}}^*,\frac{1}{\beta}({\utwi{\lambda}}^k\hspace{-1mm}-\hspace{-1mm}{\utwi{\lambda}}^{k+1}))\hspace{-1mm}-\hspace{-1mm}\gamma({\utwi{\lambda}}^k\hspace{-1mm}-\hspace{-1mm}{\utwi{\lambda}}^*) \rangle \nonumber\\
				+ &\langle\frac{1}{\alpha_1}({\bf x}^k\hspace{-1mm}-\hspace{-1mm}{\bf x}^{k+1})+\nabla f_1({\bf x}^{k+1})-\nabla f_1({\bf x}^k), {\bf x}^{k+1}-{\bf x}^*\rangle \nonumber\\
				+& \langle\nabla g_1({\bf y}^{k+1})\hspace{-1mm}-\hspace{-1mm}\nabla g_1({\bf y}^k)+\beta{\bf B}^T{\bf B}({\bf y}^{k+1}\hspace{-1mm}-\hspace{-1mm}{\bf y}^k), {\bf y}^{k+1}\hspace{-1mm}-\hspace{-1mm}{\bf y}^* \rangle \nonumber\\
				+&\langle \frac{1}{\alpha_2}({\bf y}^k-{\bf y}^{k+1}), {\bf y}^{k+1}-{\bf y}^*\rangle
				\geq \Phi. \label{cross_all}
			\end{align}
			We then proceed to split the cross terms and one can notice that there are similar terms for ${\bf x},{\bf y},{\utwi{\lambda}}$ in the following form:
			\begin{align*}
				&\langle{\utwi{\lambda}}^{k+1}-{\utwi{\lambda}}^*,\frac{1}{\beta}({\utwi{\lambda}}^k-{\utwi{\lambda}}^{k+1})\rangle + \langle \frac{1}{\alpha_1}({\bf x}^k\hspace{-1mm}-\hspace{-1mm}{\bf x}^{k+1}), {\bf x}^{k+1}\hspace{-1mm}-\hspace{-1mm}{\bf x}^*\rangle\\
				&+\frac{1}{\alpha_2}({\bf y}^k-{\bf y}^{k+1}), {\bf y}^{k+1}-{\bf y}^*\rangle
			\end{align*}
			We group them together and define the following quantities:
			$${\bf G} = \left(\begin{matrix}
			\frac{1}{\alpha_1}{\bf I}&0&0\\0&\frac{1}{\alpha_2}{\bf I}&0\\0&0&\frac{1}{\beta}{\bf I}
			\end{matrix}\right),{\bf w}^k = \left(\begin{matrix}
			{\bf x}^k\\{\bf y}^k\\ {\mbox{\boldmath $\lambda$}}^k
			\end{matrix}\right) , \,\, {\bf w}^* = \left(\begin{matrix}
			{\bf x}^*\\{\bf y}^*\\ {\mbox{\boldmath $\lambda$}}^*
			\end{matrix}\right),$$

			so that one can rewrite the inequality as follows:
			\begin{align}
				&({\bf w}^{k+1} - {\bf w}^*)^T{\bf G}({\bf w}^{k} - {\bf w}^{k+1}) + \gamma\langle {\mbox{\boldmath $\lambda$}}^*-{\mbox{\boldmath $\lambda$}}^k, {\mbox{\boldmath $\lambda$}}^{k+1} - {\mbox{\boldmath $\lambda$}}^*\rangle\nonumber\\
				+&\langle {\mbox{\boldmath $\lambda$}}^{k}\hspace{-1mm} - \hspace{-1mm}{\mbox{\boldmath $\lambda$}}^{k+1}, {\bf A}({\bf x}^{k+1}-{\bf x}^{k})\rangle + \beta\gamma\langle {\mbox{\boldmath $\lambda$}}^*\hspace{-1mm}-\hspace{-1mm}{\mbox{\boldmath $\lambda$}}^k, {\bf A}({\bf x}^{k+1}-{\bf x}^{k})\rangle \nonumber\\
				+&\langle\beta{\bf B}^T{\bf B}({\bf y}^{k+1} - {\bf y}^k),{\bf y}^{k+1}-{\bf y}^* \rangle\nonumber\\
				+&\langle{\bf x}^{k+1}\hspace{-1mm}-\hspace{-1mm}{\bf x}^*,\nabla f_1(x^{k+1})\hspace{-1mm} -\hspace{-1mm} \nabla f_1({\bf x}^k)\rangle\nonumber\\
				+&\langle{\bf y}^{k+1}\hspace{-1mm}-\hspace{-1mm}{\bf y}^*,\nabla g_1(y^{k+1})\hspace{-1mm} -\hspace{-1mm} \nabla g_1({\bf y}^k)\rangle\hspace{-1mm}\geq \Phi \, . \label{eqn:crossx}
			\end{align}
			To tackle the cross terms related only to ${\utwi{\lambda}}$, we consider the following equality
			\begin{align}
				\|{\bf a}\hspace{-1mm}-\hspace{-1mm}{\bf c}\|^2_{\bf G}\hspace{-1mm} - \hspace{-1mm}\|{\bf b}\hspace{-1mm}-\hspace{-1mm}{\bf c}\|^2_{\bf G}
				\hspace{-1mm}=\hspace{-1mm} 2({\bf a}\hspace{-1mm}-\hspace{-1mm}{\bf c})^T{\bf G}({\bf a}\hspace{-1mm}-\hspace{-1mm}{\bf b}) - \|{\bf a}\hspace{-1mm}-\hspace{-1mm}{\bf b}\|^2_{\bf G},\label{eqn:fact}
			\end{align}
			and use it in \eqref{eqn:crossx} to arrive at the following inequality:
			\begin{align}
				&({\bf w}^{k+1} - {\bf w}^*)^T{\bf G}({\bf w}^{k} - {\bf w}^{k+1}) \geq \frac{\gamma}{2}\|{\utwi{\lambda}}^{k+1}-{\utwi{\lambda}}^*\|^2\nonumber\\
				-& \frac{\gamma}{2}\|{\utwi{\lambda}}^{k}\hspace{-1mm}-\hspace{-1mm}{\utwi{\lambda}}^{k+1}\|^2\hspace{-1mm} + \hspace{-1mm}\frac{\gamma}{2}\|{\utwi{\lambda}}^{k}\hspace{-1mm}-\hspace{-1mm}{\utwi{\lambda}}^*\|^2
				\hspace{-1mm}+\hspace{-1mm}\langle {\mbox{\boldmath $\lambda$}}^{k+1}\hspace{-1mm} - \hspace{-1mm}{\mbox{\boldmath $\lambda$}}^{k}, {\bf A}({\bf x}^{k+1}\hspace{-1mm}-\hspace{-1mm}{\bf x}^{k})\rangle\nonumber\\
				+& \beta\gamma\langle {\mbox{\boldmath $\lambda$}}^k\hspace{-1mm}-\hspace{-1mm}{\mbox{\boldmath $\lambda$}}^*\hspace{-1mm}, {\bf A}({\bf x}^{k+1}\hspace{-1mm}-\hspace{-1mm}{\bf x}^{k})\rangle\hspace{-1mm}+\hspace{-1mm}\langle{\bf x}^{*}\hspace{-1mm}-\hspace{-1mm}{\bf x}^{k+1}\hspace{-1mm},\nabla\hspace{-1mm} f_1(x^{k+1})\hspace{-1mm} -\hspace{-1mm} \nabla\hspace{-1mm} f_1({\bf x}^k)\rangle\hspace{-1mm} \nonumber\\
				+&\langle{\bf y}^{*}\hspace{-1mm}-\hspace{-1mm}{\bf y}^{k+1},\nabla g_1(y^{k+1})\hspace{-1mm} -\hspace{-1mm} \nabla g_1({\bf y}^k)\rangle\hspace{-1mm}\nonumber\\
				+&\langle\beta{\bf B}^T{\bf B}({\bf y}^{k+1} - {\bf y}^k),{\bf y}^{k+1}-{\bf y}^* \rangle+\Phi.\label{eqn:crossy}
			\end{align}
			Then, we utilize the Cauchy-Schwarz inequality to bound the following term:
			\begin{align}
				&\langle {\bf A}({\bf x}^{k+1} -{\bf x}^k),{\mbox{\boldmath $\lambda$}}^{k+1}-{\mbox{\boldmath $\lambda$}}^k\rangle \nonumber\\
				\geq& -\frac{1}{2\rho_1}\|{\bf A}({\bf x}^{k+1}\hspace{-1mm}-\hspace{-1mm}{\bf x}^k)\|^2 \hspace{-1mm}-\hspace{-1mm}\frac{\rho_1}{2}\|{\mbox{\boldmath $\lambda$}}^{k+1} \hspace{-1mm}-\hspace{-1mm}{\mbox{\boldmath $\lambda$}}^k\|^2, \forall \rho_1 >0\label{cross3}
			\end{align}
			We implement same process for the rest of cross terms:
			\begin{align*}
				&\beta\gamma\langle {\bf A}({\bf x}^{k+1} -{\bf x}^k),{\mbox{\boldmath $\lambda$}}^k-{\mbox{\boldmath $\lambda$}}^*\rangle \nonumber\\
				\geq& -\frac{\beta\gamma}{2\rho_2}\|{\bf A}({\bf x}^{k+1}\hspace{-1mm}-\hspace{-1mm}{\bf x}^k)\|^2\hspace{-1mm} -\hspace{-1mm}\frac{\beta\gamma\rho_2}{2}\|{\mbox{\boldmath $\lambda$}}^{k} \hspace{-1mm}-\hspace{-1mm} {\mbox{\boldmath $\lambda$}}^*\|^2, \forall \rho_2 >0\label{cross4}\\
				&\langle\beta{\bf B}^T{\bf B}({\bf y}^k- {\bf y}^{k+1}),{\bf y}^{k+1}-{\bf y}^* \rangle\nonumber\\
				\geq&-\frac{\beta}{\rho_3}\|{\bf B}({\bf y}^{k+1} - {\bf y}^{k})\|^2 - \beta\rho_3\|{\bf B}({\bf y}^{k+1}-{\bf y}^*)\|^2\nonumber\\
				\geq& -\frac{\beta\sigma_{\max}^2({\bf B})}{\rho_3}\|{\bf y}^{k+1}\hspace{-1mm} -\hspace{-1mm} {\bf y}^{k}\|^2-\beta\rho_3\sigma_{\max}^2({\bf B})\|{\bf y}^{k+1}\hspace{-1mm}-\hspace{-1mm}{\bf y}^*\|^2 \, ,
			\end{align*}
			where $\sigma_{\max}({\bf B})$ is the maximal singular value of ${\bf B}$. The remaining terms in~\eqref{eqn:crossy} are related only to gradients of $f_1, g_1$, which can be bounded as follows:
			\begin{align}
				&\langle{\bf x}^{*}-{\bf x}^{k+1},\nabla f_1(x^{k+1}) - \nabla f_1({\bf x}^k)\rangle
				\nonumber\\
				&+\langle{\bf y}^{*}-{\bf y}^{k+1},\nabla g_1(y^{k+1}) - \nabla g_1({\bf y}^k)\rangle\nonumber \\
				\geq& -\frac{L_f^2}{2\rho_4}\|{\bf x}^{k+1}-{\bf x}^k\|^2 - \frac{\rho_4}{2}\|{\bf x}^{k+1}-{\bf x}^*\|^2 \nonumber\\
				&-\frac{L_g^2}{2\rho_5}\|{\bf y}^{k+1}-{\bf y}^k\|^2 - \frac{\rho_5}{2}\|{\bf y}^{k+1}-{\bf y}^*\|^2,
			\end{align}
			where we have used the Cauchy-Schwarz inequality and we leveraged the Lipschitz continuity of $f_1,g_1$. Also, from \eqref{eqn:fact}, it can be noticed that:
			\begin{align*}
				&\|{\bf w}^k-{\bf w}^*\|^2_{\bf G} - \|{\bf w}^{k+1}-{\bf w}^*\|^2_{\bf G} \\
				= ~&2({\bf w}^k-{\bf w}^*)^TG({\bf w}^k-{\bf w}^{k+1}) - \|{\bf w}^{k} - {\bf w}^{k+1}\|^2_{\bf G} \, .
			\end{align*}
			It therefore follows that:
			\begin{align*}
				&\|{\bf w}^k\hspace{-1mm}-\hspace{-1mm}{\bf w}^*\|^2_{\bf G} \hspace{-1mm}-\hspace{-1mm} \|{\bf w}^{k+1}\hspace{-1mm}-\hspace{-1mm}{\bf w}^*\|^2_{\bf G}  \nonumber\\
				\geq& \|{\bf w}^k\hspace{-1mm}-\hspace{-1mm}{\bf w}^{k+1}\|^2_{\bf G} + \gamma\|{\utwi{\lambda}}^{k+1}-{\utwi{\lambda}}^*\|^2 \nonumber\\
				-& \gamma\|{\utwi{\lambda}}^{k}-{\utwi{\lambda}}^{k+1}\|^2 + \gamma\|{\utwi{\lambda}}^{k}-{\utwi{\lambda}}^*\|^2\nonumber\\
				-&\frac{\sigma_{\max}^2({\bf A})}{\rho_1}\|{\bf x}^{k+1}\hspace{-1mm}-\hspace{-1mm}{\bf x}^k\|^2 \hspace{-1mm}-\hspace{-1mm}\rho_1\|{\mbox{\boldmath $\lambda$}}^{k+1} \hspace{-1mm}-\hspace{-1mm}{\mbox{\boldmath $\lambda$}}^k\|^2\nonumber\\
				-&\frac{\beta\gamma\sigma_{\max}^2({\bf A})}{\rho_2}\|{\bf x}^{k+1}\hspace{-1mm}-\hspace{-1mm}{\bf x}^k\|^2\hspace{-1mm} -\hspace{-1mm}\beta\gamma\rho_2\|{\mbox{\boldmath $\lambda$}}^{k} \hspace{-1mm}-\hspace{-1mm} {\mbox{\boldmath $\lambda$}}^*\|^2\nonumber\\
				-&\frac{L_f^2}{\rho_4}\|{\bf x}^{k+1}-{\bf x}^k\|^2 - \rho_4\|{\bf x}^{k+1}-{\bf x}^*\|^2 \nonumber\\
				-&\frac{L_g^2}{\rho_5}\|{\bf y}^{k+1}-{\bf y}^k\|^2 - \rho_5\|{\bf y}^{k+1}-{\bf y}^*\|^2+ \Phi\nonumber \\
				-&\frac{\beta\sigma_{\max}^2({\bf B})}{\rho_3}\|{\bf y}^{k+1}\hspace{-1mm} -\hspace{-1mm} {\bf y}^{k}\|^2\hspace{-1mm}-\hspace{-1mm}\beta\rho_3\sigma_{\max}^2({\bf B})\|{\bf y}^{k+1}\hspace{-1mm}-\hspace{-1mm}{\bf y}^*\|^2
			\end{align*}
			Rearranging the terms in a suitable way, we arrive at the following inequality:
			\begin{align}
				&\|{\bf w}^k\hspace{-1mm}-\hspace{-1mm}{\bf w}^*\|^2_{\bf G} \hspace{-1mm}-\hspace{-1mm} \|{\bf w}^{k+1}\hspace{-1mm}-\hspace{-1mm}{\bf w}^*\|^2_{\bf G}
				\nonumber\\
				\geq& (\frac{1}{\alpha_1}\hspace{-1mm}-\hspace{-1mm} \frac{\sigma_{\max}^2({\bf A})}{\rho_1} \hspace{-1mm}-\hspace{-1mm} \frac{\beta\gamma\sigma_{\max}^2({\bf A})}{\rho_2} - \frac{L_f^2}{\rho_4})\|{\bf x}^k-{\bf x}^{k+1}\|^2\nonumber\\
				+&(\frac{1}{\alpha_2} - \beta\frac{\sigma_{\max}^2({\bf B})}{\rho_3} - \frac{L_g^2}{\rho_5})\|{\bf y}^k-{\bf y}^{k+1}\|^2\nonumber\\
				+&(\frac{1}{\beta}-\gamma-\rho_1)\|{\utwi{\lambda}}^k-{\utwi{\lambda}}^{k+1}\|^2+(2v_f-\rho_4)\|{\bf x}^{k+1}-{\bf x}^*\|^2\nonumber\\
				+&(2v_g-\rho_5-\beta\rho_3\sigma_{\max}^2({\bf B}))\|{\bf y}^{k+1}-{\bf y}^*\|^2 \nonumber\\
				+&\gamma\|{\utwi{\lambda}}^{k+1}-{\utwi{\lambda}}^*\|^2 + (\gamma-\beta\gamma\rho_2)\|{\utwi{\lambda}}^k-{\utwi{\lambda}}^*\|^2\label{cross2}.
			\end{align}
		Recall that the goal is to prove the following inequality:
		\begin{align}
			\|{\bf w}^k-{\bf w}^*\|^2_{\bf G} \geq (1+\delta)\|{\bf w}^{k+1}-{\bf w}^*\|^2_{\bf G},\label{linconv1}
		\end{align}
		where $\delta>0$ is a constant. For brevity, denote the right-hand-side of \eqref{cross2} as $C$; then it is sufficient to prove that:
		\begin{align*}
			&C\geq \delta\|{\bf w}^{k+1}\hspace{-1mm}-\hspace{-1mm}{\bf w}^*\|^2_{\bf G} \nonumber\\
			&\hspace{-3mm}=\hspace{-1mm}\frac{\delta}{\alpha_1}\|{\bf x}^{k+1}\hspace{-1mm}-\hspace{-1mm}{\bf x}^*\|^2 \hspace{-1mm}+\hspace{-1mm} \frac{\delta}{\alpha_2}\|{\bf y}^{k+1}\hspace{-1mm}-\hspace{-1mm}{\bf y}^*\|^2\hspace{-1mm}+\hspace{-1mm}\frac{\delta}{\beta}\|{\mbox{\boldmath $\lambda$}}^{k+1}\hspace{-1mm}-\hspace{-1mm}{\mbox{\boldmath $\lambda$}}^*\|^2 \,
		\end{align*}
		which  requires  the following to hold true:
		\begin{align*}
			&\frac{1}{\alpha_1} - \frac{\max\sigma^2({\bf A})}{\rho_1} - \frac{\beta\gamma\max\sigma^2({\bf A})}{\rho_2} - \frac{L_f^2}{\rho_4} \geq 0\\
			&\frac{1}{\alpha_2} - \frac{\beta\max\sigma^2({\bf B})}{\rho_3} - \frac{L_g^2}{\rho_5}\geq 0,~ \frac{1}{\beta}-\gamma-\rho_1 \geq 0 \\
			&2v_f-\rho_4-\frac{\delta}{\alpha_1} \geq 0,~ 2v_g - \rho_5-\beta\rho_3\max\sigma^2({\bf B}) -\frac{\delta}{\alpha_2} \geq 0\\
			&\gamma-\frac{\delta}{\beta} \geq 0,~ \gamma-\beta\gamma\rho_2\geq 0.
		\end{align*}
		From the inequalities above, one can notice  that the constant $\alpha_1,\alpha_2$ is closely related to various constants as well as the  singular values of ${\bf A}$ and ${\bf B}$, denoted as $\sigma({\bf A})$ and $\sigma({\bf B})$, respectively. Specifically, this leads to the following conditions for the step sizes:
		\begin{align}
			&\frac{\delta}{2v_f\hspace{-1mm}-\hspace{-1mm}\rho_4} \leq \alpha_1 \leq \frac{1}{(\frac{1}{\rho_1}\hspace{-1mm}+\hspace{-1mm}\frac{\beta\gamma}{\rho_2})\max\sigma^2({\bf A})\hspace{-1mm} + \hspace{-1mm}\frac{L_f^2}{\rho_4}},  \label{eqn:alpha1}\\
			&\frac{\delta}{2v_g\hspace{-1mm}-\hspace{-1mm}\rho_5\hspace{-1mm}-\hspace{-1mm}\beta\rho_3\max\sigma^2({\bf B})}\leq \alpha_2\leq \frac{1}{\frac{\beta\max\sigma^2({\bf B})}{\rho_3} \hspace{-1mm}+\hspace{-1mm} \frac{L_g^2}{\rho_5}} \, . \label{eqn:alpha2}
		\end{align}
		To ensure that there exists step sizes $\alpha_1,\alpha_2$ that satisfy the condition above, one can choose $\rho_1 = 1, \rho_2 = 1, \rho_3 = \frac{v_g}{2\beta\max\sigma^2({\bf B})}, \rho_4 = v_f, \rho_5 = v_g$, and from Assumption \ref{Assum_Bounded},\ref{Assum_uniform} we know we have all problem dependent parameters uniform bounded.

		As for other parameters, one has that:
		\begin{align*}
			\frac{1}{\beta}-\gamma-1 \geq 0,~ \gamma-\beta\gamma\geq 0\Rightarrow~ \beta\gamma+\beta \leq 1,~ \beta\leq 1.
		\end{align*}
		This completes the proof.
	\end{proof}
	Using  Lemma~\ref{lem:linear}, we can now proven the result of  Theorem~\ref{main}.  For notation simplicity, we define $r = \frac{1}{1+\delta}$. Using~\eqref{eqn:lin} and the triangle inequality, we have that
	\begin{align*}
		&\|{\bf w}^{(k)} \hspace{-1mm}-\hspace{-1mm} {\bf w}^{*,(k)}\|_{\bf G} \leq r\|{\bf w}^{(k-1)}\hspace{-1mm} - \hspace{-1mm} {\bf w}^{*,(k-1)} \hspace{-1mm}+\hspace{-1mm} {\bf w}^{*,(k-1)} \hspace{-1mm}-\hspace{-1mm} {\bf w}^{*,(k)}\|_{\bf G} \\
		& \leq r\|{\bf w}^{(k-1)}\hspace{-1mm} - \hspace{-1mm} {\bf w}^{*,(k-1)}\|_{\bf G} +  r\|{\bf w}^{*,(k-1)} - {\bf w}^{*,(k)}\|_{\bf G} \, .
	\end{align*}
	We now find a bound on $\|{\bf w}^{*,(k-1)} \hspace{-1mm}-\hspace{-1mm} {\bf w}^{*,(k)}\|_{\bf G}$.
	From \eqref{diff}, \eqref{dual_opt} and the fact that we have chosen $\gamma=1$, we know that
	\textcolor{black}{
		\begin{align*}
			&{\bf A}^{(k+1)}{\bf x}^{*,(k+1)}\hspace{-1mm}-\hspace{-1mm}{\bf A}^{(k)}{\bf x}^{*,(k)}\hspace{-1mm} +\hspace{-1mm} {\bf B}^{(k+1)}{\bf y}^{*,(k+1)}\hspace{-1mm}-\hspace{-1mm}{\bf B}^{(k)}{\bf y}^{*,(k)}\nonumber\\
			+&{\bf b}^{(k)}-{\bf b}^{(k+1)}+ {\utwi{\lambda}}^{*,(k+1)} -{\utwi{\lambda}}^{*,(k)} = 0.
		\end{align*}
		Move ${\utwi{\lambda}}$ terms to the other side of the equation and take norm for both hand sides, it follows:
		\begin{align*}
			&\|{\utwi{\lambda}}^{*,(k+1)} - {\utwi{\lambda}}^{*,(k)}\| \\
			\leq& \|{\bf A}^{(k+1)}({\bf x}^{*,(k+1)} \hspace{-1mm}-\hspace{-1mm} {\bf x}^{*,(k)}) + ({\bf A}^{(k+1)}-{\bf A}^{(k)}){\bf x}^{*,(k)}\| \\
			&+\|{\bf B}^{(k+1)}({\bf y}^{*,(k+1)} \hspace{-1mm}-\hspace{-1mm} {\bf y}^{*,(k)}) + ({\bf B}^{(k+1)}-{\bf B}^{(k)}){\bf y}^{*,(k)}\|\\
			&+\|{\bf b}^{(k)}-{\bf b}^{(k+1)}\|.
		\end{align*}
		From triangle inequality of norm we know
		\begin{align*}
			&\|{\utwi{\lambda}}^{*,(k+1)} - {\utwi{\lambda}}^{*,(k)}\| \\
			\leq& \|{\bf A}^{(k+1)}\|\|{\bf x}^{*,(k+1)} \hspace{-1mm}-\hspace{-1mm} {\bf x}^{*,(k)}\|\hspace{-1mm} +\hspace{-1mm} \|{\bf B}^{(k+1)}\|\|{\bf y}^{*,(k+1)}\hspace{-1mm} -\hspace{-1mm} {\bf y}^{*,(k)}\| \\ &+ \|{\bf b}^{(k)}-{\bf b}^{(k+1)}\|\hspace{-1mm}+\hspace{-1mm} \|{\bf A}^{(k+1)}-{\bf A}^{(k)}\|\|{\bf x}^{*,(k)}\|\\
			&+\|{\bf B}^{(k+1)}-{\bf B}^{(k)}\|\|{\bf y}^{*,(k)}\|
		\end{align*}
		From \eqref{nedic2} we know that $\|{\bf x}^{*,(k)}\|\leq \sigma_1 + \sqrt{\gamma}\max\|{\utwi{\lambda}}^{\text{opt}, (k)}\|\cdot c := \mathcal{x}(\sigma_1), \|{\bf y}^{*,(k)}\|\leq \sigma_2 + \sqrt{\gamma}\max\|{\utwi{\lambda}}^{\text{opt}, (k)}\|\cdot c:= \mathcal{x}(\sigma_2)$
		Combining with Assumption~\ref{Assum_Bounded} we can reach the following inequality
		\begin{align*}
			&\|{\utwi{\lambda}}^{*,(k+1)} - {\utwi{\lambda}}^{*,(k)}\| \\
			\leq&\tilde{\sigma}_{{\bf A}}\sigma_{\bf x}+\tilde{\sigma}_{{\bf B}}\sigma_{\bf y}+\sigma_{\bf b}+\sigma_{\bf A}\mathcal{x}(\sigma_1) + \sigma_{\bf B}\mathcal{x}(\sigma_2)
			\triangleq \sigma_{{\utwi{\lambda}}},
		\end{align*}
		which gives us
		\begin{align}
			\|{\bf w}^{*,(k-1)} - {\bf w}^{*,(k)}\|_{\bf G} \leq \psi :=  \sqrt[]{\frac{\sigma_{\bf x}^2}{\alpha_1} +\frac{\sigma_{\bf y}^2}{\alpha_2} + 2\sigma_{{\utwi{\lambda}}}^2}
		\end{align}
		and the desired result is obtained. }

	\subsection{Proof of Corollary \ref{linear}}

	By recursively applying the  Theorem~\ref{main}, we have that
	\begin{align*}
		\|{\bf w}^{(k)} \hspace{-1mm}-\hspace{-1mm} {\bf w}^{*,(k)}\|_{\bf G}
		& \leq r^k\|{\bf w}^{(0)}-{\bf w}^{*,(0)}\|_{\bf G} \\
		& +  \sum_{i=1}^{(k)}r^{k-i+1}\|{\bf w}^{*,(i-1)}-{\bf w}^{*,(i)}\|_{\bf G}
	\end{align*}
	where $\|{\bf w}^{*,(k-1)} - {\bf w}^{*,(k)}\|_{\bf G} \leq \psi$.  Taking $k\rightarrow+\infty$, we can derive
	\begin{align*}
		&\lim_{k\rightarrow\infty}\|{\bf w}^{(k)} - {\bf w}^{*,(k)}\|_{\bf G}  \nonumber\\
		\leq&\lim_{k\rightarrow\infty}\left( \frac{r(1-r^{(k)})}{1-r}\psi(\sigma) + r^k\|{\bf w}^{(0)}-{\bf w}^{*,(0)}\|_{\bf G}\right)\\
		\Rightarrow &\lim\sup_{k\rightarrow\infty}\|{\bf w}^{(k)} - {\bf w}^{*,(k)}\|_{\bf G} \leq \frac{r}{1-r}\psi = \frac{\psi}{\delta}.
	\end{align*}
	The desired result is then obtained.

	\subsection{Remarks on Corollary 2}
	The result is obtained by choosing the biggest step sizes;  that is:
	\begin{align*}
		\alpha_1 = \frac{1}{(1+\beta\gamma)\tilde{\sigma^2}_{\bf A} + \frac{\tilde{L}_f^2}{\tilde{v}_f}}, ~\alpha_2 = \frac{1}{\frac{2\beta^2\max\tilde{\sigma}^4_{\bf B}}{\tilde{v}_g}+ \frac{\tilde{L}_g^2}{\tilde{v}_g}}
	\end{align*}
	Recall that:
	\begin{align}
		\delta &\leq \frac{\tilde{v}_f}{(1+\beta\gamma)\tilde{\sigma}^2_{\bf A} + \frac{\tilde{L}_f^2}{\tilde{v}_f}},~\delta \leq \frac{\tilde{v}_g}{\frac{4\beta^2\tilde{\sigma}^4_{\bf B}}{\tilde{v}_g} + \frac{2\tilde{L}_g^2}{\tilde{v}_g}}.\label{eqn:delta}
	\end{align}
	We already know that $\delta \leq \beta\gamma$; therefore, based on~\eqref{eqn:alpha1}--\eqref{eqn:alpha2}, one can pick $\delta$ as
	\begin{align*}
		\delta = \min \left(\frac{\tilde{v}_f}{(1+\beta\gamma)\tilde{\sigma}^2_{\bf A} + \frac{\tilde{L}_f^2}{\tilde{v}_f}},\frac{\tilde{v}_g}{\frac{4\beta^2\tilde{\sigma}^4_{\bf B}}{\tilde{v}_g} + \frac{2\tilde{L}_g^2}{\tilde{v}_g}}, \beta\gamma \right).
	\end{align*}.

	%
	\ifCLASSOPTIONcaptionsoff
	\newpage
	\fi

	\bibliographystyle{IEEEtran}

\section{\@startsection {section}{1}{\z@}%
                                   {-1.5ex \@plus -1ex \@minus -.1ex}%
                                    {0.5ex \@plus.1ex}%
                                   {\normalfont\large\bfseries}}
\renewcommand\subsection{\@startsection{subsection}{2}{\z@}%
                                     {-1.0ex\@plus -1ex \@minus -.1ex}%
                                     {0.4ex \@plus .1ex}%
                                     {\normalfont\large\bfseries}}

redpajama_set_name: RedPajamaArXiv