# A Globally Convergent Distributed Jacobi Scheme for Block-Structured
Nonconvex Constrained Optimization Problems

Anirudh Subramanyam, Youngdae Kim, Michel Schanen, François Pacaud, and Mihai
Anitescu The authors are with the Mathematics and Computer Science Division,
Argonne National Laboratory, Lemont IL 60439. This work was supported by the
U.S. Department of Energy, Office of Science, under contract DE-
AC02-06CH11357.

###### Abstract

Motivated by the increasing availability of high-performance parallel
computing, we design a distributed parallel algorithm for linearly-coupled
block-structured nonconvex constrained optimization problems. Our algorithm
performs Jacobi-type proximal updates of the augmented Lagrangian function,
requiring only local solutions of separable block nonlinear programming (NLP)
problems. We provide a cheap and explicitly computable Lyapunov function that
allows us to establish global and local sublinear convergence of our
algorithm, its iteration complexity, as well as simple, practical and
theoretically convergent rules for automatically tuning its parameters. This
in contrast to existing algorithms for nonconvex constrained optimization
based on the alternating direction method of multipliers that rely on at least
one of the following: Gauss-Seidel or sequential updates, global solutions of
NLP problems, non-computable Lyapunov functions, and hand-tuning of
parameters. Numerical experiments showcase its advantages for large-scale
problems, including the multi-period optimization of a 9000-bus AC optimal
power flow test case over 168 time periods, solved on the Summit supercomputer
using an open-source Julia code.

###### Index Terms:

distributed optimization, augmented Lagrangian, nonconvex optimization

##  I Introduction

Block-structured nonlinear optimization models are ubiquitous in science and
engineering applications. Some examples include nonlinear model predictive
control, multi-stage stochastic programming, supervised machine learning,
optimization with differential-algebraic equations, as well as network control
with geographically distributed agents. These models are used in a wide
variety of areas including power systems, telecommunications, sensor networks,
smart manufacturing, and chemical process systems, to name but a few. The
block structure in several of these models have the following form, which is
the subject of the present paper.

$\displaystyle\mathop{\text{minimize}}_{x_{1},\ldots,x_{T}}$
$\displaystyle\sum_{t=1}^{T}f_{t}(x_{t})$ (1) subject to $\displaystyle
x_{t}\in X_{t},\;\;t\in\\{1,2,\ldots,T\\},$
$\displaystyle\sum_{t=1}^{T}A_{t}x_{t}=b,$

where $X_{t}$ are compact (possibly nonconvex) sets,
$A_{t}\in\mathbb{R}^{m\times n_{t}}$ are matrices, and
$f_{t}:\mathbb{R}^{n_{t}}\mapsto\mathbb{R}$ are continuously differentiable
(possibly nonconvex) functions.

In formulation (1), the decision variables are grouped into $T$ blocks
$(x_{1},x_{2},\ldots,x_{T})$ coupled only via the linear constraints
$\sum_{t=1}^{T}A_{t}x_{t}=b$. Such linearly coupled block-structured models
arise in a wide variety of applications. For example, consider a discrete-time
dynamical system:

$x_{t+1}=\phi(x_{t},u_{t})\;,$

where $x_{t}\in\mathbb{R}^{n_{x}}$ and $u_{t}\in\mathbb{R}^{n_{u}}$ are the
vectors of states and inputs, respectively, and
$\phi:\mathbb{R}^{n_{x}}\times\mathbb{R}^{n_{u}}\mapsto\mathbb{R}^{n_{x}}$ is
some smooth state transition function. In such systems, the optimal control
problem that must be solved in the context of model predictive control [1] can
be brought into the form of (1) by introducing additional variables
$\tilde{x}_{t}\in\mathbb{R}^{n_{x}}$ and additional constraints:

$x_{t+1}=\tilde{x}_{t},\;\;\tilde{x}_{t}=\phi(x_{t},u_{t}),$

and by defining the $t^{\mathrm{th}}$ block to be
$X_{t}=\\{(x_{t},u_{t},\tilde{x}_{t})\in\mathbb{R}^{n_{x}}\times\mathbb{R}^{n_{u}}\times\mathbb{R}^{n_{x}}:\tilde{x}_{t}=\phi(x_{t},u_{t})\\}$.
The $T$ blocks are then coupled via the linear equations
$x_{t+1}=\tilde{x}_{t}$, $t=1,2,\ldots,T-1$.

The goal of this paper is to devise a parallel algorithm that can solve
instances of formulation (1) when either the number of blocks $T$ or the
number of variables/constraints in individual blocks $X_{t}$ is large. In such
cases, memory or storage requirements may prohibit the direct solution of (1)
using conventional nonlinear programming (NLP) solvers. Also, structure-
agnostic solvers may not be able to exploit any available distributed parallel
computing capabilities, resulting in excessively high computational times that
can be detrimental in several applications including model predictive control.

###  I-A Literature review

A popular approach to solving large-scale instances of (1) is to decompose the
problem iteratively into individual block subproblems that are easier to
handle numerically. This decomposition can be done either at the linear
algebra level of an interior point method [2, 3, 4] or by utilizing
(augmented) Lagrangian functions to solve block-separable dual problems that
are derived using (local) convex duality [5, 6, 7]. The alternating direction
method of multipliers (ADMM) belongs to the latter class of methods and has
seen a recent surge in popularity because of its suitability for distributed
computation; e.g., see [8]. However, the vast majority of provably convergent
ADMM approaches for solving (1) either exploit convexity in the objective
function $f_{t}$ or constraints $X_{t}$ or they are tailored for specific
instances of (1); e.g., see [9, 10, 11, 12, 13, 14].

A recent body of literature [15, 16, 17, 18, 19, 20, 21, 22, 23, 24] has
analyzed the global and local convergence of ADMM approaches for non-convex
instances of formulation (1). It is well-understood by now (e.g., see [17])
that for $T>2$, convergence requires the existence of a block variable that is
constrained only by the linear coupling equation, and whose objective function
is globally Lipschitz differentiable. Moreover, it must be coupled in such a
way that it can “control” the iterates of the coupling constraints’ dual
variables (see Remark 1 later in the paper).

One way to achieve this is to first introduce an additional block of slack
variables $z\in\mathbb{R}^{m}$ that relax the linear coupling constraints, and
then drive $z$ to $0$ by adding a smooth quadratic penalty term in the
objective function:

$\displaystyle\mathop{\text{minimize}}_{x_{1},\ldots,x_{T},z}$
$\displaystyle\sum_{t=1}^{T}f_{t}(x_{t})+\frac{\theta}{2}\left\lVert
z\right\rVert^{2}$ (2) subject to $\displaystyle
z\in\mathbb{R}^{m},\;\;x_{t}\in X_{t},\;\;t\in\\{1,2,\ldots,T\\},$
$\displaystyle\sum_{t=1}^{T}A_{t}x_{t}+z=b,$

where the penalty coefficient $\theta$ is a function of the target tolerance
$\epsilon>0$ for satisfying the linear coupling constraints. In [18], this
idea is exploited to design algorithms that globally converge to
$\epsilon$-stationary points of (1) in $O(\epsilon^{-6})$ iterations.
Moreover, the algorithms require the global solution of NLP problems [18,
Remark 3.8], which can be numerically demanding when the local subproblems are
nonconvex. A similar quadratic penalty idea is adopted in [21] although it
requires the augmented Lagrangian function (with respect to all nonlinear and
coupling constraints) to satisfy the Kurdyka-Łojasiewicz property [15];
furthermore, that analysis does not reconcile the penalty formulation (2) with
the original formulation (1). Each iteration linearizes the augmented
Lagrangian function around the previous iterate that is then alternately
minimized with respect to the $T$ blocks. Other approaches using local linear
and convex approximations of the nonlinear constraints or of the augmented
Lagrangian function, and methods establishing convergence under additional
assumptions such as coercivity of the objective function, have been proposed
in [25, 16, 26, 23].

Instead of using a quadratic penalty, [19, 20] replace the objective in (2) by
the augmented Lagrangian function

$\sum_{t=1}^{T}f_{t}(x_{t})+\beta^{\top}z+\frac{\theta}{2}\left\lVert
z\right\rVert^{2},$ (3)

and drive $z$ to $0$ by iteratively updating the Lagrange multipliers $\beta$,
where each so-called outer loop iteration uses another inner-level iterative
ADMM scheme that minimizes (3) subject to the constraints of (2). An
alternative two-level decomposition for (possibly nonlinear) coupling
constraints is proposed in [22] but it depends intimately on randomized block
updating of the augmented Lagrangian function which in turn requires locally
tight upper bounds of the objective function with respect to each variable
block, while fixing the other blocks.

All of the aforementioned algorithms rely on a Gauss-Seidel updating scheme,
where the variable blocks are iteratively optimized in a particular sequence.
Specifically, when updating a given variable block (say $x_{t}$), the values
of all other blocks that appear earlier in the sequence (e.g.,
$x_{1},\ldots,x_{t-1}$) must be fixed to their newest values. The overall
convergence relies critically on the sequential nature of this update.
Unfortunately, this makes these algorithms unsuitable for distributed parallel
computations. A common workaround (e.g., see [6, 27]) to enable parallel
computation is to equivalently reformulate (1), (2) by introducing additional
variables $y_{t}$ as follows:

$\begin{array}[]{r@{\;\;}l}\displaystyle\mathop{\text{minimize}}_{\begin{subarray}{c}x_{1},\ldots,x_{T}\\\
y_{1},\ldots,y_{T}\end{subarray}}&\displaystyle\sum_{t=1}^{T}f_{t}(x_{t})\\\
\text{subject to}&\displaystyle x_{t}\in
X_{t},\;\;A_{t}x_{t}=y_{t},\;\;t\in\\{1,2,\ldots,T\\},\\\
&\displaystyle\sum_{t=1}^{T}y_{t}=b,\end{array}$ (4)

The augmented Lagrangian function with respect to the last equation is
completely decomposable in the $x$-variables:

$\sum_{t=1}^{T}f_{t}(x_{t})+\beta^{\top}\Big{(}\sum_{t=1}^{T}y_{t}-b\Big{)}+\frac{\theta}{2}\Big{\|}\sum_{t=1}^{T}y_{t}-b\Big{\|}^{2}\;.$
(5)

Therefore, the reformulation (4) can be interpreted as a two-block model where
the $x$\- and $y$-variables constitute respectively the first and second
blocks. Since (5) admits a closed-form minimization with respect to $y$ (for
fixed values of $x$), one can use any of the aforementioned Gauss-Seidel
algorithms to enable parallel implementation. However, this reformulation
substantially increases the problem dimension in terms of both the number of
variables and constraints. This can slow down convergence and we demonstrate
this empirically when we compare it with our proposed method.

Recently, [16] suggested an alternative strategy to enable parallel
computation. Similar to the prox-linear method [28], it proposes to linearize
around the previous iterate the quadratic penalty term in the augmented
Lagrangian function, which makes the latter decomposable with respect to the
$x$ blocks. Although the method is shown to asymptotically converge for (2)
(no convergence rate type is provided), it can be slow to converge in practice
for models with highly nonconvex constraints. Indeed, that analysis suggests
an augmented Lagrangian penalty parameter that scales as $O(\theta^{4})$
posing numerical challenges. Empirically, that method failed to converge for
the smallest of our test instances with $T=3$.

###  I-B Contributions

We propose a novel distributed Jacobi scheme for solving the block-structured
optimization problem (1) with smooth nonconvex objective functions $f_{t}$ and
constraint sets $X_{t}$ that are algebraically described by smooth nonconvex
functions $c_{t}$. The scheme only requires local solutions of individual
block NLP problems. In contrast to solving the higher-dimensional
reformulation (4), our algorithm directly performs Jacobi updates of the
augmented Lagrangian function with respect to the linear coupling constraints
in (2), in which optimizing a single variable block does not require the
newest values of the other blocks. This makes it particularly suitable for
parallel computation, where the cost of each iteration can be reduced by a
factor of $T$ compared to existing Gauss-Seidel schemes.

We show that the algorithm converges globally to an $\epsilon$-approximate
stationary point of (1) in no more than $O(\epsilon^{-4})$ iterations and
locally converges to an approximate local minimizer at a sublinear rate under
mild assumptions. The proof uses an easy-to-compute Lyapunov function that
does not require a (typically unknown) local minimizer. We supplement the
analysis by showing (empirically) that convergence fails if the proximal
weights in the algorithm are not chosen appropriately. The algorithm can be
interpreted as a nonconvex constrained extension of the Jacobi algorithms
proposed in [29, 30, 27, 14], and as a Jacobi extension of the Gauss-Seidel
algorithm for nonconvex problems proposed in [18]. In contrast to the majority
of existing approaches for nonconvex problems, we provide a practical and
automatic parameter tuning scheme, an open-source Julia implementation (that
can be downloaded from https://github.com/exanauts/ProxAL.jl), and an
empirical demonstration of convergence on a large-scale multi-period optimal
power flow test case.

###  I-C Notation, Assumptions and Preliminaries

For any integer $T$, we use $[T]$ to denote the index set
$\\{1,2,\ldots,T\\}$. We use $x$ (without subscript) as shorthand for the
entire vector of decisions $(x_{1},\ldots,x_{T})\in\mathbb{R}^{n}$, and $X$ as
the corresponding shorthand for the feasible set $X_{1}\times\ldots\times
X_{T}$, where we define $n\coloneqq\sum_{t=1}^{T}n_{t}$. For an arbitrary
vector $w=(w_{1},\ldots,w_{T})\in\mathbb{R}^{n}$, we define
$Aw\coloneqq\sum_{t=1}^{T}A_{t}w_{t}$, and for any $t\in[T]$, we define
$A_{\neq t}w_{\neq t}\coloneqq\sum_{s\in[T]\setminus\\{t\\}}A_{s}w_{s}$. We
let $D\coloneqq\mathrm{diag}(A_{1},\ldots,A_{T})\in\mathbb{R}^{Tm\times n}$
denote the block-diagonal matrix with $A_{1}$, $\ldots$, $A_{T}$ along its
diagonal. Iterates at the $k^{\text{th}}$ iteration are denoted with
superscript $k$. The general normal cone [31, Definition 6.3] to a set $X$ at
a point $x\in X$ is denoted as $N_{X}(x)$. For a vector $z$ and matrix $M$, we
use $\left\lVert z\right\rVert$ and $\left\lVert M\right\rVert$ to denote
their Euclidean and spectral norms, respectively. When $M$ is positive
(semi-)definite, we use $\left\lVert z\right\rVert_{M}$ to denote the
(semi-)norm $\sqrt{z^{\top}Mz}$. For a vector $z\in\mathbb{R}^{N}$ and set
$S\subseteq\mathbb{R}^{N}$, we define
$\mathop{\mathrm{dist}}\left(z,S\right)\coloneqq\min_{w\in S}\left\lVert
z-w\right\rVert$. For any $\epsilon>0$, we let $B_{\epsilon}(z)$ denote the
open Euclidean ball in $\mathbb{R}^{N}$ with center $z$ and radius $\epsilon$.
We use $\mathrm{I}$ and $\mathrm{0}$ to denote the identity and zero matrices,
respectively; unless indicated otherwise, their dimensions should be clear
from the context. For two square symmetric matrices $M_{1},M_{2}$, we use
$M_{1}\succ M_{2}$ ($M_{1}\succeq M_{2}$) or $M_{2}\prec M_{1}$ ($M_{2}\preceq
M_{2}$) to indicate that $M_{1}-M_{2}$ is positive definite (positive
semidefinite).

Throughout the paper, we make the following assumptions.

  1. (A1)

$X_{t}$ is non-empty and compact for all $t\in[T]$.

  2. (A2)

$f_{t}:\mathbb{R}^{n_{t}}\mapsto\mathbb{R}$ is $C^{2}$ for all $t\in[T]$.

  3. (A3)

The matrix $A\coloneqq[A_{1}\ldots A_{T}]\in\mathbb{R}^{m\times n}$ has full
row rank.

  4. (A4)

Problem (1) has a feasible solution.

We say $x^{*}\in X$ is a stationary point of (1), if there exist Lagrange
multipliers $\lambda^{*}\in\mathbb{R}^{m}$ satisfying:

$\displaystyle Ax^{*}=b,$ (6a) $\displaystyle\nabla
f_{t}(x_{t}^{*})+A_{t}^{\top}\lambda^{*}\in-N_{X_{t}}(x_{t}^{*}),\;\;t\in[T].$
(6b)

By introducing the primal and dual residual functions,
$\pi:X\mapsto\mathbb{R}$ and
$\delta_{t}:X_{t}\times\mathbb{R}^{m}\mapsto\mathbb{R}$, respectively,
conditions (6) can be stated as $\pi(x^{*})=0$ and
$\delta_{t}(x_{t}^{*},\lambda^{*})=0$, $t\in[T]$, where

$\displaystyle\pi(x)$ $\displaystyle=\left\lVert Ax-b\right\rVert,$ (7a)
$\displaystyle\delta_{t}(x_{t},\lambda)$
$\displaystyle=\mathop{\mathrm{dist}}\left(\nabla
f_{t}(x_{t})+A_{t}^{\top}\lambda,-N_{X_{t}}(x_{t})\right),t\in[T]$ (7b)

Under an appropriate constraint qualification, the above are equivalent to the
Karush-Kuhn-Tucker (KKT) conditions, and they must be necessarily satisfied if
$x^{*}$ is a local solution of (1); see [32, Lemma 12.9] for a discussion of
the linear independence constraint qualification (LICQ).

The augmented Lagrangian function of (2) with respect to its coupling
constraints is parameterized with penalty parameters $\rho,\theta>0$ and is
defined as follows (with domain $X\times\mathbb{R}^{m}\times\mathbb{R}^{m}$):

$\displaystyle\mathcal{L}(x,z,\lambda)$
$\displaystyle=\sum_{t=1}^{T}f_{t}(x_{t})+\frac{\theta}{2}\left\lVert
z\right\rVert^{2}$ (8)
$\displaystyle\phantom{=}+\lambda^{\top}\left[Ax+z-b\right]+\frac{\rho}{2}\left\lVert
Ax+z-b\right\rVert^{2}.$

For convenience, we let $\mathcal{L}(x_{t};\bar{x}_{\neq
t},\bar{z},\bar{\lambda})$ denote the above function with domain $X_{t}$ that
is obtained by fixing all $x_{s}$ variables to $\bar{x}_{s}$,
$s\in[T]\setminus\\{t\\}$, $z$ to $\bar{z}$, $\lambda$ to $\bar{\lambda}$, and
ignoring any constants:

$\displaystyle\mathcal{L}(x_{t};\bar{x}_{\neq t},\bar{z},\bar{\lambda})$
$\displaystyle=f_{t}(x_{t})+\bar{\lambda}^{\top}A_{t}x_{t}$ (9)
$\displaystyle\phantom{=}+\frac{\rho}{2}\left\lVert A_{t}x_{t}+A_{\neq
t}\bar{x}_{\neq t}+\bar{z}-b\right\rVert^{2}.$

##  II Algorithm

Algorithm 1 outlines the basic distributed scheme for solving formulation (1).
Here, $x^{0},z^{0}$ and $\lambda^{0}$ denote the initial guesses of the primal
and dual variables of the penalty formulation (2), whereas $\rho$, $\theta$
and $\tau_{x},\tau_{z}$ are scalars denoting the penalty parameters in (2) and
its augmented Lagrangian function, and proximal weights corresponding to the
$x$ and $z$ variables, respectively.

Algorithm 1 Basic distributed proximal Jacobi scheme

1:$x^{0}\in X$, $z^{0}\in\mathbb{R}^{m}$, $\lambda^{0}\in\mathbb{R}^{m}$,
scalars $\rho,\theta,\tau_{x},\tau_{z}>0$.

2:for $k=1,2,\ldots$

3: Update $x$: compute a local minimizer $x_{t}^{k}$ of (10) by warm-starting
with $x_{t}^{k-1}$ and solve in parallel for $t\in[T]$:  $\min_{x_{t}\in
X_{t}}\mathcal{L}(x_{t};x_{\neq
t}^{k-1},z^{k-1},\lambda^{k-1})+\frac{\tau_{x}}{2}\left\lVert
x_{t}-x_{t}^{k-1}\right\rVert_{A_{t}^{\top}A_{t}}^{2}$ (10)

4: Update $z$:
$z^{k}=\frac{\tau_{z}z^{k-1}-\rho\left[Ax^{k}-b\right]-\lambda^{k-1}}{\tau_{z}+\rho+\theta}$
(11)

5: Update $\lambda$:
$\lambda^{k}=\lambda^{k-1}+\rho\left[Ax^{k}+z^{k}-b\right]$ (12)

6:end for

We note a few important points about Algorithm 1. The optimization problems
(10) in line 3 that compute new values of the $x$ variables can be solved
completely in parallel using any local NLP solver. These problems minimize the
sum of the augmented Lagrangian function (8) and the proximal terms, over the
$x_{t}$ variables alone. Similarly, the update step (11) is the closed-form
solution of the quadratic problem that is defined by minimizing the sum of the
augmented Lagrangian function (8) and the proximal terms, over the $z$
variables alone, and by fixing the $x$ and $\lambda$ variables to $x^{k}$ and
$\lambda^{k-1}$, respectively. Even though this update occurs sequentially
after solving (10), our analysis and results remain unchanged if we use a
variant of the algorithm where the former is performed in parallel using only
information about $x^{k-1}$. We omit this variant for ease of exposition.
Indeed, since the update step (11) is relatively cheap compared to solving
problem (10), we can compute the former along with the latter locally on each
node of a parallel computing architecture. Finally, the dual variables in (12)
are updated as per the standard augmented Lagrangian method.

###  II-A Global convergence

We first establish that under an appropriate choice of the scalar parameters
$\rho,\theta,\tau_{x},\tau_{z}$, the sequence $\\{x^{k}\\}$ generated by
Algorithm 1 converges to a point satisfying the first-order stationarity
conditions (6). The key idea is to establish that the sum of the augmented
Lagrangian function and proximal terms is a candidate Lyapunov function:

$\displaystyle\Phi(x,z,\lambda,\hat{x},\hat{z})=\mathcal{L}(x,z,\lambda)$
$\displaystyle+\frac{\tau_{z}}{4}\left\lVert z-\hat{z}\right\rVert^{2}$ (13)
$\displaystyle+\sum_{t=1}^{T}\frac{\tau_{x}}{4}\left\lVert
x_{t}-\hat{x}_{t}\right\rVert^{2}_{A_{t}^{\top}A_{t}}.$

Specifically, we show that the sequence $\\{\Phi^{k}\\}$, which consists of
$\Phi$ evaluated at the iterates generated by the algorithm, decreases
monotonically and is bounded from below by $\hat{\Phi}$:

$\displaystyle\Phi^{k}$
$\displaystyle\coloneqq\Phi(x^{k},z^{k},\lambda^{k},x^{k-1},z^{k-1}),\;k\geq
1,$ (14) $\displaystyle\hat{\Phi}$ $\displaystyle\coloneqq\min_{x\in
X}\sum_{t=1}^{T}f_{t}(x_{t}).$ (15)

In the remainder of the paper, for $k\geq 1$, we define $\Delta z^{k}\coloneqq
z^{k}-z^{k-1}$ and similarly define $\Delta x^{k}$, $\Delta\lambda^{k}$, and
$\Delta\Phi^{k}$. Also, let $\Delta
z^{0}\coloneqq-\tau_{z}^{-1}(\lambda^{0}+\theta z^{0})$, $\Delta
x^{0}\coloneqq 0$ and
$\Phi^{0}\coloneqq\mathcal{L}(x^{0},z^{0},\lambda^{0})+\frac{\tau_{z}}{4}\|\Delta
z^{0}\|^{2}$. Finally, we let $\pi^{k}\coloneqq\pi(x^{k})$ and
$\delta^{k}_{t}\coloneqq\delta_{t}(x_{t}^{k},\lambda^{k})$.

######  Theorem 1 (Global convergence).

Suppose that assumptions (A1)–(A4) hold, and the sequence
$\\{(x^{k},z^{k},\lambda^{k})\\}$ is generated by Algorithm 1 with parameters
$\rho,\theta,\tau_{x},\tau_{z}>0$, such that

$\eta_{x}\coloneqq\frac{\tau_{x}}{4}-\frac{(T-1)\rho}{2}>0,\;\eta_{z}\coloneqq\frac{\tau_{z}}{4}-\frac{2(\theta+\tau_{z})^{2}}{\rho}>0.$
(16)

Then, for all $K\geq 1$, there exists $j\in[K]$ such that

$\displaystyle\pi^{j}\leq\sqrt{\frac{2\big{(}\Phi^{1}-\hat{\Phi}\big{)}}{\theta}\left(1+\frac{2(\theta+\tau_{z})^{2}}{K\eta_{z}\rho}\right)}$
$\displaystyle\delta_{t}^{j}\leq(\rho+\tau_{x})\left\lVert
A_{t}\right\rVert\sqrt{\frac{2(T+1)\big{(}\Phi^{1}-\Phi^{K}\big{)}}{K\min\\{\eta_{x},\eta_{z}\\}}},\;t\in[T].$

In particular, for any $\epsilon\in(0,1)$, if we choose

$\theta=\frac{1}{\epsilon^{2}},\quad\rho=\frac{64}{\epsilon^{2}},\quad\tau_{x}=\frac{256\,(T-1)}{\epsilon^{2}},\quad\tau_{z}=\frac{2}{\epsilon^{2}},\quad$

then after $K=O(\epsilon^{-4})$ iterations, the iterates
$\\{(x^{k},z^{k},\lambda^{k})\\}$ converge to an $\epsilon$-stationary point
of problem (1) satisfying
$\min_{j\in[K]}\max\\{\pi^{j},\max_{t\in[T]}\delta_{t}^{j}\\}=O(\epsilon)$.

###### Proof.

See Appendix A. ∎

######  Remark 1.

The iteration complexity can be improved to $O(\epsilon^{-2})$, whenever the
original problem (1) has the same form as problem (2); see also [18, 19, 17]
for related discussions. Specifically, we can avoid introducing variables $z$
with penalty weights $O(\epsilon^{-2})$, whenever formulation (1) contains a
variable block, say $x_{T}$, that is constrained only via the linear coupling
constraints; i.e., $X_{T}=\mathbb{R}^{n_{T}}$, and its coupling matrix
satisfies $A_{T}=\mathrm{I}$ (or more generally, its image
$\mathop{\mathrm{Im}}\left(A_{T}\right)\supseteq\mathop{\mathrm{Im}}\left([b,A_{1},\ldots,A_{T-1}]\right)$),
and its objective satisfies $\nabla^{2}f_{T}(x_{T})\preceq M\mathrm{I}$ for
all $x\in X$, for some fixed $M>0$.

###  II-B Local convergence

We now establish that the sequence $\\{x^{k}\\}$ generated by Algorithm 1
converges to an approximate local minimizer of the original problem (1), under
some additional assumptions. In contrast to the previous section, where we
showed global convergence to stationary points of the original problem (1),
the proof proceeds by showing convergence to local minimizers of the quadratic
penalty formulation (2). To relate local minimizers of the former to those of
the latter, we introduce the (full) Lagrangian functions of problems (1) and
(2), $\Lambda$ and $\Lambda_{\theta}$, respectively. Under assumption (A5)
stated below, these functions have domains
$\mathbb{R}^{n}\times\mathbb{R}^{r}\times\mathbb{R}^{m}$ and
$\mathbb{R}^{n}\times\mathbb{R}^{r}\times\mathbb{R}^{m}\times\mathbb{R}^{m}$,
respectively, where $r\coloneqq\sum_{t=1}^{T}r_{t}$, and are given by:

$\displaystyle\Lambda(x,\mu,\lambda)=\sum_{t=1}^{T}[f_{t}(x_{t})+\mu_{t}^{\top}c_{t}(x_{t})]+\lambda^{\top}[Ax-b]$
(17)
$\displaystyle\Lambda_{\theta}(x,z,\mu,\lambda)=\Lambda(x,\mu,\lambda)+\frac{\theta}{2}\left\lVert
z\right\rVert^{2}+\lambda^{\top}z$ (18)

where $\mu_{t}\in\mathbb{R}^{r_{t}}$ and $\mu$ denotes the vector
$(\mu_{1},\ldots,\mu_{T})\in\mathbb{R}^{r}$. We make the following additional
assumptions.

  1. (A5)

$X_{t}=\\{x_{t}\in\mathbb{R}^{n_{t}}:c_{t}(x_{t})=0\\}$ can be algebraically
described using only equality constraints, where
$c_{t}:\mathbb{R}^{n_{t}}\mapsto\mathbb{R}^{r_{t}}$ is $C^{2}$ for all
$t\in[T]$.

  2. (A6)

For all $\theta>0$,111In what follows, we suppress dependence of
$x^{*},z^{*},\mu^{*},\lambda^{*}$ on $\theta$ for ease of exposition. there
exist $x^{*}\in X$ and $z^{*}\in\mathbb{R}^{m}$, such that $(x^{*},z^{*})$ is
a local solution of problem (2). Also, the constraint Jacobian $J(x^{*})$ of
the original problem (1) at $x^{*}$ has full row rank, where we define:

$\small J(x)=\begin{bmatrix}\nabla c(x)^{\top}\\\ A\end{bmatrix},\nabla
c(x)=\begin{bmatrix}\nabla c_{1}(x_{1})&\ldots&0\\\ \vdots&\ddots&\vdots\\\
0&\ldots&\nabla c_{T}(x_{T})\end{bmatrix}.$

  3. (A7)

For all $\theta>0$, there exist $\mu^{*}\in\mathbb{R}^{r}$ and
$\lambda^{*}\in\mathbb{R}^{m}$, such that for the same $x^{*},z^{*}$ in
assumption (A6), we have

$\displaystyle
w_{x}^{\top}\nabla^{2}_{xx}\Lambda_{\theta}(x^{*},z^{*},\mu^{*},\lambda^{*})w_{x}+\theta
w_{z}^{\top}w_{z}>0,$

for all $(w_{x},w_{z})\in\mathbb{R}^{n}\times\mathbb{R}^{m}$ satisfying
$\nabla c(x^{*})^{\top}w_{x}=0$, $Aw_{x}+w_{z}=0$ and $\left\lVert
w_{x}\right\rVert+\left\lVert w_{z}\right\rVert>0$.

  4. (A8)

For all $k\geq 1$, the constraint Jacobian of problem (10) at its optimal
solution, $\nabla c_{t}(x_{t}^{k})^{\top}$, has full row rank.

  5. (A9)

For all $k\geq 1$, the $x$-update step in line 3 computes a local minimizer of
problem (10) that is closest to either $x_{t}^{k-1}$ if $A_{t}$ has full
column rank or to $x_{t}^{*}$ otherwise.

We make some remarks about these assumptions. Assumption (A5) is without loss
of generality since inequalities can be reformulated as equality constraints
by adding squares of additional slack variables, and it is also used in the
convergence analyses of existing augmented Lagrangian methods (e.g., see [5,
19, 24]). Assumptions (A6) and (A7) state that the linear independence
constraint qualification (LICQ) and second-order sufficient conditions (SOSC)
are satisfied, respectively, at the local solution $(x^{*},z^{*})$ of the
quadratic penalty formulation (2) for all (sufficiently large) $\theta>0$.
Assumption (A8) states that LICQ is also satisfied at local solutions of
problem (10) that correspond to the $x$-updates. As before, these are fairly
standard to establish local convergence of augmented Lagrangian methods.
Finally, Assumption (A9) is satisfied by any well-behaved NLP solver if
$A_{t}$ has full column rank since all variables have an associated proximal
term in problem (10) in this case; when $A_{t}$ does not have full column
rank, the assumption is necessary for convergence and it is implicit in
classical analyses of the standard augmented Lagrangian method [5, see eq.(6)
in Proposition 2.4].

######  Remark 2.

Assumption (A9) can be relaxed if the proximal term in problem (10) is
replaced with $\frac{\tau_{x}}{2}\|x-x^{k-1}\|^{2}_{P_{t}}$, where
$P_{t}\succeq A_{t}^{\top}A_{t}$ has full column rank. It can be shown that
local (and global) convergence continues to hold in this case (with minor
adjustments to the statement of Theorem 2) if we relax (A9) to require only
that problem (10) computes a local minimizer that is closest to $x_{t}^{k-1}$.
We do not present this generalization for ease of exposition.

The key idea is to show that the iterates $(x^{k},z^{k})$ generated by
Algorithm 1 correspond to local minimizers of a perturbed variant of problem
(2), that is parameterized by $p\in\mathbb{R}^{m}$ and
$d=(d_{1},\ldots,d_{T},d_{z})\in\mathbb{R}^{n+m}$. These parameters correspond
to the primal and dual residuals of the quadratic penalty formulation (2) at
the iterates of Algorithm 1. Specifically, for any $k\geq 1$, we define these
as follows:

$\displaystyle p^{k}$ $\displaystyle\coloneqq Ax^{k}+z^{k}-b$ (19a)
$\displaystyle d_{t}^{k}$ $\displaystyle\coloneqq\rho A_{t}^{\top}A_{\neq
t}\Delta x_{\neq t}^{k}-\rho A_{t}^{\top}\Delta
z^{k}-\tau_{x}A_{t}^{\top}A_{t}\Delta x_{t}^{k}$ (19b) $\displaystyle
d_{z}^{k}$ $\displaystyle\coloneqq-\tau_{z}\Delta z^{k}$ (19c)

In the following theorem, recall that
$D=\mathop{\mathrm{diag}}(A_{1},\ldots,A_{T})$ so that
$Dx=(A_{1}x_{1},\ldots,A_{T}x_{T})\in\mathbb{R}^{Tm}$.

######  Theorem 2 (Local convergence).

Suppose that assumptions (A1)–(A9) hold, $\epsilon\in(0,1)$ is a fixed
constant, and parameters $\rho,\theta,\tau_{x},\tau_{z}>0$ are chosen such
that $\theta=\epsilon^{-2}$ and (16) holds. Then, there exist constants
$\bar{\rho},\epsilon^{\prime}>0$ such that if $\rho>\bar{\rho}$ and the
sequence $\\{(x^{k},z^{k},\lambda^{k})\\}$ is generated by Algorithm 1 with
input $(x^{0},z^{0},\lambda^{0})$ satisfying
$\left\lVert(Dx^{0},z^{0},\lambda^{0})-(Dx^{*},z^{*},\lambda^{*})\right\rVert<\epsilon^{\prime}$
and $\lambda^{0}=-\theta z^{0}$, then the sequence $\\{x^{k}\\}$ converges at
a sublinear rate to the $\epsilon$-approximate local minimizer $x^{*}$ of
problem (1) satisfying:

$\displaystyle\max\left\\{\pi(x^{*}),\max_{t\in[T]}\delta_{t}(x_{t}^{*},\lambda^{*})\right\\}=O(\epsilon),$
$\displaystyle
w^{\top}\nabla^{2}_{xx}\Lambda(x^{*},\mu^{*},\lambda^{*})w>0,\;\forall
w:J(x^{*})w=0,w\neq 0.$

###### Proof.

See Appendix B. ∎

###  II-C Adaptive parameter tuning

In practice, the values for the parameters $\rho,\theta,\tau_{x},\tau_{z}$
that are suggested in Theorems 1 and 2 are quite conservative and can be
significantly larger than what is required for convergence. Therefore,
Algorithm 2 presents a practical strategy for adaptively tuning these
parameter values.

Algorithm 2 Adaptive distributed proximal Jacobi scheme

1:$\epsilon\in(0,1)$, $\rho_{0},\omega,\kappa_{x},\kappa_{z},\zeta,\Psi>0$,
$\nu_{x},\nu_{\rho},\nu_{\theta},\chi>1$

2:Initialize $\theta=\epsilon^{-2}$, $\rho=\rho_{0}$,
$\tau_{x}=\kappa_{x}\rho$, $\tau_{z}=\kappa_{z}\rho$, $\psi=0$.

3:for $k=1,2,\ldots$

4: Update $x$, $z$, $\lambda$ as per Algorithm 1

5: if $\Phi^{k}-\Phi^{k-1}>\zeta|\Phi^{k}|$

6: $\tau_{x}\leftarrow\min\\{\nu_{x}\tau_{x},(2T-1)\rho\\}$

7: end if

8: if
$\max\left\\{\|{p^{k}}\|_{\infty},\|{d^{k}}\|_{\infty}\right\\}\leq\epsilon$
and $\|{Ax^{k}-b}\|_{\infty}>\epsilon$

9: $\theta\leftarrow\nu_{\theta}\theta$

10: end if

11: if $\|{p^{k}}\|_{\infty}>\chi\|{d^{k}}\|_{\infty}$ and $\rho<\omega\theta$

12: $\rho\leftarrow\min\\{\nu_{\rho}\rho,\omega\theta\\}$,
$\tau_{x}\leftarrow\kappa_{x}\rho$, $\tau_{z}\leftarrow\kappa_{z}\rho$

13: else if $\|{d^{k}}\|_{\infty}>\chi\|{p^{k}}\|_{\infty}$ and $\psi<\Psi$

14: $\rho\leftarrow\rho/\nu_{\rho}$, $\tau_{x}\leftarrow\kappa_{x}\rho$,
$\tau_{z}\leftarrow\kappa_{z}\rho$, $\psi\leftarrow\psi+1$

15: end if

16: if $\|{Ax^{k}-b}\|_{\infty}\leq\epsilon$ stop end if

17:end for

The parameter tuning rules are motivated from the proofs of Theorems 1 and 2.
In particular, the (cheaply computable) difference in Lyapunov function
values, $\Phi^{k}-\Phi^{k-1}$, is guaranteed to be negative, whenever the
penalty parameters $\rho,\theta$ and proximal weights $\tau_{x},\tau_{z}$ are
large enough to ensure $\eta_{x},\eta_{z}>0$. Therefore, Algorithm 2 starts
with relatively small values of these parameters and gradually adjusts them.
In particular, $\tau_{x}$ is increased in line 6 by a factor $\nu_{x}>1$
whenever $\Delta\Phi^{k}$ (appropriately scaled by $|\Phi^{k}|$) is larger
than some $\zeta>0$; however, the increase is limited to $(2T-1)\rho$ based on
the definition (16) of $\eta_{x}$. The penalty parameter $\theta$ is increased
by a factor $\nu_{\theta}>1$ in line 9 whenever the primal and dual residuals
with respect to quadratic penalty formulation (2), $p^{k}$ and $d^{k}$, are
smaller than the tolerance $\epsilon$, but the true primal residual $Ax^{k}-b$
continues to be larger than $\epsilon$.

The parameter $\rho$ is increased or decreased by a factor $\nu_{\rho}>1$ to
ensure that the magnitudes of $p^{k}$ and $d^{k}$ are within a factor $\chi>1$
of each other, similar to [8]. However, the number of times $\rho$ is allowed
to decrease is bounded by some integer $\Psi>0$, and moreover, it is not
allowed to increase after $\rho=\omega\theta$. Therefore, $\rho$ will not be
adjusted an infinite number of times. The proximal weight $\tau_{z}$ is fixed
at $\kappa_{z}\rho$, where $\kappa_{z}$ is selected to ensure that
$\eta_{z}>0$ whenever $\rho$ is also sufficiently large. In particular, the
proof of Theorem 1 shows that whenever $\rho/\theta=\omega>32$, then
$\kappa_{z}=1/32$ suffices. Finally, since Theorem 1 ensures that the
termination condition in line 16 will be met after a finite number of
iterations (e.g., when $\eta_{x},\eta_{z}>0$), the parameter $\theta$ will not
be increased an infinite number of times in line 9. Therefore, Algorithm 2 is
guaranteed to converge to an $\epsilon$-approximate stationary point of the
original problem (1) after $O(\epsilon^{-4})$ iterations.

##  III Numerical Experiments

We demonstrate the computational performance of the proximal Jacobi scheme on
multi-period AC Optimal Power Flow (ACOPF) problems. The ACOPF is used to
determine optimal dispatch levels of generators to satisfy electrical loads in
power systems [33]. Its multi-period variant entails the solution of multiple
single-period ACOPF problems that are coupled over a long time horizon. Such
extended horizons are essential to model operational constraints such as
ramping limits of thermal generators, as well as planning decisions involving
energy storage units [34] or production costing with accurate representations
of system voltages [35].

We consider a $T$-period ACOPF problem, where the variables
$x_{t}=(p^{g}_{t},q^{g}_{t},V_{t},\vartheta_{t})$ in time period $t$ are the
real ($p^{g}_{t}$) and reactive ($q^{g}_{t}$) power generation levels and
nodal voltage magnitudes ($V_{t}$) and angles ($\vartheta_{t}$), respectively.
The constraint set $X_{t}$ captures AC power balances with respect to the real
($p^{d}_{t}$) and reactive loads ($q^{d}_{t}$) in time period $t$, whereas the
objective function $f_{t}(p^{g}_{t})$ minimizes the total production cost in
that time period. For simplicity, we ignore constraints on transmission line
flow limits. The single-period ACOPF problems are coupled via inter-temporal
constraints that capture the physical ramping limits of generators; that is,
generator $i$ cannot change its real power output between periods $t$ and
$t+1$ by more than $r_{i}\Delta_{t}$ units, where $\Delta_{t}$ denotes the
length of period $t$. Formally, if $G$ and $B$ denote the sets of generators
and buses, respectively, then this problem can be formulated as follows.

$\displaystyle\mathop{\text{minimize}}_{x_{1},\ldots,x_{T}}$
$\displaystyle\sum_{t\in[T]}f_{t}(p^{g}_{t})$ subject to $\displaystyle
x_{t}=(p^{g}_{t},q^{g}_{t},V_{t},\vartheta_{t})\in[\underline{x},\overline{x}],\;t\in[T],$
$\displaystyle\sum_{j\in
G_{i}}p^{g}_{jt}-p^{d}_{it}=c_{i}^{\mathrm{re}}(V_{t},\vartheta_{t}),\;i\in
B,\;t\in[T],$ $\displaystyle\sum_{j\in
G_{i}}q^{g}_{jt}-q^{d}_{it}=c_{i}^{\mathrm{im}}(V_{t},\vartheta_{t}),\;i\in
B,\;t\in[T],$ $\displaystyle\left\lvert
p^{g}_{i,t+1}-p^{g}_{it}\right\rvert\leq r_{i}\Delta_{t},\;i\in
G,\;t\in[T-1],$

where $[\underline{x},\overline{x}]$ represent variable bounds, $G_{i}$ is the
set of generators connected to bus $i$, and $c_{i}^{\mathrm{re}}$ and
$c_{i}^{\mathrm{im}}$ are smooth nonconvex functions of $V_{t}$ and
$\vartheta_{t}$. We denote by $Y=Y^{\mathrm{re}}+\sqrt{-1}Y^{\mathrm{im}}$ the
$|B|\times|B|$ admittance matrix (with $Y_{ii}=y_{ii}-\sum_{j\neq i}Y_{ij}$,
$y_{ii}=y^{\mathrm{re}}_{ii}+\sqrt{-1}y^{\mathrm{im}}_{ii}$, where
$y^{\mathrm{re}}_{ii}$ and $y^{\mathrm{im}}_{ii}$ are the shunt conductance
and susceptance at bus $i\in B$, respectively). Using the polar formulation,
the functions $c_{i}^{\mathrm{re},\mathrm{im}}$ can be written, for all buses
$i\in B$, as [36]

$\displaystyle c_{i}^{\mathrm{re}}(V,\vartheta)$
$\displaystyle=Y^{\mathrm{re}}_{ii}V_{i}^{2}$ $\displaystyle+\sum_{j\in
N_{i}}V_{i}V_{j}\big{[}Y^{\mathrm{re}}_{ij}\cos(\vartheta_{i}-\vartheta_{j})+Y^{\mathrm{im}}_{ij}\sin(\vartheta_{i}-\vartheta_{j})\big{]},$
$\displaystyle c_{i}^{\mathrm{im}}(V,\vartheta)$
$\displaystyle=-Y^{\mathrm{im}}_{ii}V_{i}^{2}$ $\displaystyle+\sum_{j\in
N_{i}}V_{i}V_{j}\big{[}Y^{\mathrm{re}}_{ij}\sin(\vartheta_{i}-\vartheta_{j})-Y^{\mathrm{im}}_{ij}\cos(\vartheta_{i}-\vartheta_{j})\big{]},$

where $N_{i}\subseteq B$ is the set of neighboring buses of $i$.

The multi-period ACOPF can be formulated as an instance of problem (1) by
introducing extra variables $s^{g}_{i,t+1}\in[0,2r_{i}\Delta_{t}]$ and
reformulating the inter-temporal ramping limits as equality constraints:
$p^{g}_{i,t+1}-p^{g}_{it}+s^{g}_{i,t+1}=r_{i}\Delta_{t}$. Note that although
the objective function is linear or convex quadratic, the nonconvexity of
$c^{\mathrm{re}}$ and $c^{\mathrm{im}}$ makes the problem nonconvex.

We consider the standard IEEE test cases ‘118’, ‘1354pegase’, and ‘9241pegase’
available from MATPOWER [37]. Since the original data only provides a single
vector of loads, we generate a load profile over multiple periods as follows.
We first obtain hourly load data from New England ISO [38] over a typical week
(7 days) of operations. We then use this $T=24\times 7=168$ period profile as
a multiplier for the load at each bus, downscaling if required to ensure
single-period ACOPF feasibility. Thus, the load distribution (across buses) in
each period remains the same as the original MATPOWER case, and only the total
load (across time periods) follows the imposed profile. We consider fairly
stringent ramp limits $r_{i}\in[0.33\%,0.50\%]\times p^{\max}_{i}$ (per
minute, which is the standard unit for this application class), where
$p^{\max}_{i}$ is the maximum rated output of generator $i$ (note that
$\Delta_{t}=60$ min). For each test case, the objective function is scaled by
$10^{-3}$ (roughly $T^{-1}$) to ensure $\rho=O(1)$. All runs are initialized
with $x^{0}=(\underline{x}+\overline{x})/2$, $z^{0}=0$, $\lambda^{0}=0$ and
unless mentioned otherwise, following parameter values are used in Algorithm
2: $\rho_{0}=10^{-3},\kappa_{x}=2$ for case 118 and
$\rho_{0}=10^{-5},\kappa_{x}=2.5$ otherwise, and $\omega=32$,
$\kappa_{z}=1/32$, $\zeta=10^{-4}$, $\nu_{x}=2$, $\nu_{\rho}=2$,
$\nu_{\theta}=10$, $\chi=10$, $\Psi=100$.

Our algorithm is implemented in Julia and is available at
https://github.com/exanauts/ProxAL.jl. It uses the Message Passing Interface
(MPI) which allows the use of distributed parallel computing resources. All
NLP subproblems were solved using the JuMP modeling interface [39] and Ipopt
[40] (with default options) as the NLP solver. The computational times
reported in Section III-C were obtained on the Summit supercomputer at Oak
Ridge National Laboratory [41].

###  III-A Role of proximal terms

The proof of Theorem 1 shows that for fixed values of $\theta$ and $\rho$, the
Lyapunov $\\{\Phi^{k}\\}$ must be monotonically decreasing, whenever the
proximal weights $\tau_{x}$ are sufficiently large. To verify this
empirically, we fix $\theta=10^{6}$, $\rho=1$, $\tau_{z}=\kappa_{z}\rho$ and
then consider $\tau_{x}\in\\{0,1,2\\}$; we also disable their automatic
updates in Algorithm 2. Figure 1 shows the values of the corresponding
Lyapunov sequences
$\\{\big{(}\Phi^{k}-\underline{\Phi}\big{)}/{\underline{\Phi}}\\}$, where
$\underline{\Phi}$ is a normalizing constant equal to the smallest observed
value of $\Phi^{k}$ (across 100 iterations). We find that when the proximal
terms are absent ($\tau_{x}=0$) or when they are present but not sufficiently
large ($\tau_{x}=1$), the Lyapunov sequence (and hence also
$\\{(x^{k},z^{k})\\}$) diverges. This is contrast to classical Gauss-Seidel
based algorithms [19, 18], which are not as critically dependent on proximal
terms. Indeed, the Jacobi nature of the $x$-update in Algorithm 1 requires
that the new iterate $x^{k}$ be sufficiently close to its old value $x^{k-1}$
for convergence.

Figure 1: Lyapunov sequence for increasing values of the proximal weight
$\tau_{x}$ and for fixed values of $\theta$, $\rho$ and $\tau_{z}$.

###  III-B Comparison with existing method

We compare the performance of our proposed scheme with the two-level ADMM
algorithm proposed in [19]. Since the latter uses Gauss-Seidel or sequential
updating rules, which are not amenable to distributed parallelization, the two
algorithms are not directly comparable. Therefore, to solve the multi-period
ACOPF with the two-level algorithm, we reformulate it using the variable
splitting technique described in [19, Section 1.2]. Specifically, for each
$t\in[T]$, we introduce a so-called global copy $\bar{p}^{g}_{t}$ of
$p^{g}_{t}$, and for each $t\in[T-1]$, we introduce a so-called local copy
$\hat{p}^{g}_{t+1}$ (local to block $t+1$) of $p^{g}_{t}$. The ramping limit
for period $t\in[T-1]$ is then expressed as:
$p^{g}_{i,t+1}-\hat{p}^{g}_{i,t+1}+s^{g}_{i,t+1}=r_{i}\Delta_{t}$, which is
added as an additional constraint to $X_{t+1}$. Finally, the coupling
constraints linking the various blocks can be expressed as:
$p^{g}_{t}=\bar{p}^{g}_{t}$ for $t\in[T]$, and
$\hat{p}^{g}_{t+1}=\bar{p}^{g}_{t}$ for $t\in[T-1]$. For brevity, we consider
only the 118-bus case with a ramp limit of $0.33\%$.

As suggested in [19, Section 6], our implementation of the two-level ADMM uses
$\omega=0.75$, $\gamma=1.5$,
$[\underline{\lambda},\bar{\lambda}]=[-10^{6},10^{6}]$, $\beta^{1}=1000$,
$\rho^{k}=c_{1}\beta^{k}$, and the $k^{\mathrm{th}}$ level inner-level loop is
terminated based on [19, eq.14c] with
$\epsilon_{3}^{k}=\sqrt{2n_{g}(T-1)}/(c_{2}k\rho^{k})$, where $n_{g}$ is the
number of generators and all other symbols refer to the notation in [19].
After trying $(c_{1},c_{2})\in\\{2\times
10^{f}:f=-6,-4,-2,0,2\\}\times\\{1,1000\\}$, we set $c_{1}=2\times 10^{-4}$
and $c_{2}=1000$, since it produced the smallest infinity norm of the combined
primal and dual residual vector [19, eq.14a-14c] after 50 cumulative
iterations. We set the tolerance $\epsilon=10^{-8}$ in Algorithm 2.

Figure 2 shows the performance of our proximal Jacobi scheme and the two-level
algorithm. All residuals correspond to infinity norms; in our scheme, the
plotted primal and dual residuals are $\|Ax^{k}-b\|_{\infty}$ and
$\|d^{k}\|_{\infty}$; in the two-level scheme, they are
$\max\\{\|p^{g,k}-\bar{p}^{g,k}\|_{\infty},\|\hat{p}^{g,k}-\bar{p}^{g,k}\|_{\infty}\\}$
and the infinity norm of [19, eq.14a-14b] (plotted as a function of the number
of cumulative iterations), respectively. We observe that for the first few
iterations (up to $25$) and for small target tolerances (up to $10^{-3}$), the
two algorithms have similar behavior in terms of reducing the optimality
errors. However, our proposed scheme is able to decrease the residuals at an
almost linear rate even for tighter tolerances, whereas the decrease for the
two-level scheme seems to become worse than linear. This is likely because
Algorithm 2 allows $\rho$ to decrease; indeed, we observe in Figure 2 that the
primal and dual residuals are roughly within a factor of $\chi=10$ of each
other at all iterations. Finally, note that the run times of both schemes are
expected to be similar, since the per-iteration cost in either scheme is
dominated by the solution of (single-period) ACOPF NLP problems. The actual
run time is a function of the number of parallel processes available which we
study in the next section.

Figure 2: Computational comparison of the proposed proximal Jacobi scheme with
the two-level algorithm proposed in [19].

###  III-C Scalability

One of the main features of the proposed algorithm is its ability to make use
of distributed parallel computing resources. Therefore, we study its
scalability with respect to both the problem size (for a fixed number of
parallel processes) as well as the number of processes (for a fixed problem
size). Table I summarizes the computational performance for the different test
cases using $T=168$ parallel processes and a target tolerance of
$\epsilon=10^{-3}$ in each case. We observe from Table I that the number of
iterations and run times increase with increasing network size and with
decreasing values of the ramping limit (which makes the linear coupling
constraints harder to satisfy). Nevertheless, the total run time remains less
than 100 minutes for the most difficult test case consisting of more than 3
million variables and constraints. Also, Figure 3 (which is plotted for the
ramping limit of $0.33\%$) shows that although the initial primal residuals
are quite large, they are reduced by more than 4 orders of magnitude over the
course of roughly 50 iterations. Finally, Figure 4 shows the empirically
observed run times as a function of the number of parallel processes (MPI
ranks) for the 118-bus test case with $r=0.33\%$. We find that our
implementation exhibits a near-linear scaling in this particular instance.

TABLE I: Summary of computational performance using $T=168$ parallel processes for convergence to $\epsilon=10^{-3}$ tolerance. Case | # Vars. | # Cons. | Ramp % | # Iters. | Time (s)
---|---|---|---|---|---
118 | 66,810 | 48,666 | 0.33 | 24 | 5.2
0.50 | 13 | 3.9
1354pegase | 1,181,779 | 696,259 | 0.33 | 60 | 137.3
0.50 | 53 | 119.7
9241pegase | 3,235,756 | 3,148,396 | 0.33 | 67 | 4,755.5
0.50 | 59 | 3,511.9
Figure 3: Convergence of the proposed scheme to a tolerance
$\epsilon=10^{-3}$. Figure 4: Scalability as a function of the number of
parallel processes.

##  IV Conclusions

This paper proposed a distributed parallel decomposition algorithm for solving
linearly coupled block-structured nonconvex constrained optimization problems.
Such structures naturally arise in several applications, including nonlinear
model predictive control and stochastic optimization. The algorithm performs
Jacobi updates of the augmented Lagrangian function with appropriately chosen
proximal terms, requiring only local solutions of the individual block NLP
problems. By constructing an easily computable Lyapunov function, we showed
that the algorithm converges globally to $\epsilon$-approximate stationary
points in $O(\epsilon^{-4})$ iterations as well as locally to
$\epsilon$-approximate local minimizers of the original NLP at a sublinear
rate. We also provided a simple, practical and theoretically convergent
variant of the algorithm where its parameters are adaptively tuned during the
iterations. Our numerical experiments show that it can outperform existing
algorithms based on ADMM and that optimizing large-scale nonconvex AC optimal
power flow problems over more than 100 time steps can be done in the order of
roughly one hour. Future work includes the study of acceleration techniques
such as inexact NLP solutions and momentum methods, integration with the two-
level algorithm of [19], as well as the incorporation of second-order
information such as in [6].

##  Appendix A Global convergence proof

The proof of Theorem 1 requires several intermediate results. Throughout, we
suppose that the conditions in Theorem 1 hold; that is, assumptions (A1)–(A4)
are satisfied and $\rho,\theta,\eta_{x},\eta_{z}>0$ are chosen such that (16)
is satisfied.

First, Lemma 2 establishes that the Lyapunov sequence $\\{\Phi^{k}\\}$ defined
in (14) decreases by positive multiples of $\left\lVert\Delta
x^{k}_{t}\right\rVert^{2}_{A_{t}^{\top}A_{t}}$ and $\left\lVert\Delta
z^{k}\right\rVert^{2}$. Second, Lemma 3 shows that the sequence is bounded
below by $\hat{\Phi}$ defined in (15). Lemma 4 uses these results to bound
$\left\lVert\Delta x^{k}_{t}\right\rVert^{2}_{A_{t}^{\top}A_{t}}$ and
$\left\lVert\Delta z^{k}\right\rVert^{2}$ as a function of the iteration index
$k$. Finally, Propositions 1 and 2 establish bounds on the primal and dual
residuals, respectively. The proof of Theorem 1 follows directly from these.
We use the following intermediate result to simplify exposition.

######  Lemma 1.

For all $k\geq 0$, we have:

$\lambda^{k}=-\theta z^{k}-\tau_{z}\Delta z^{k}.$ (20)

###### Proof.

For $k=0$, equation (20) is equivalent to the definition of $\Delta z^{0}$.
For $k\geq 1$, the $\lambda$-update formula (12) implies that
$\lambda^{k-1}=\lambda^{k}-\rho\big{[}\sum_{t=1}^{T}A_{t}x_{t}^{k}+z^{k}-b\big{]}.$
Substituting this in the $z$-update formula (11) yields equation (20). ∎

######  Lemma 2.

For all $k\geq 1$, we have:

$\displaystyle\Delta\Phi^{k}\leq$
$\displaystyle-\eta_{x}\sum_{t=1}^{T}\left(\left\lVert\Delta
x_{t}^{k}\right\rVert^{2}_{A_{t}^{\top}A_{t}}+\left\lVert\Delta
x_{t}^{k-1}\right\rVert^{2}_{A_{t}^{\top}A_{t}}\right)$ (21)
$\displaystyle-\eta_{z}\left(\left\lVert\Delta
z^{k}\right\rVert^{2}+\left\lVert\Delta z^{k-1}\right\rVert^{2}\right).$

###### Proof.

Using the definitions (13) and (14), we have:

$\displaystyle\Delta\Phi^{k}$
$\displaystyle=\mathcal{L}(x^{k},z^{k},\lambda^{k})-\mathcal{L}(x^{k-1},z^{k-1},\lambda^{k-1})$
(22)
$\displaystyle\phantom{=}+\sum_{t=1}^{T}\frac{\tau_{x}}{4}\left[\left\lVert\Delta
x_{t}^{k}\right\rVert_{A_{t}^{\top}A_{t}}^{2}-\left\lVert\Delta
x_{t}^{k-1}\right\rVert_{A_{t}^{\top}A_{t}}^{2}\right]$
$\displaystyle\phantom{=}+\frac{\tau_{z}}{4}\left[\left\lVert\Delta
z^{k}\right\rVert^{2}-\left\lVert\Delta z^{k-1}\right\rVert^{2}\right].$

The first term in the above expression can be written as:

$\displaystyle\mathcal{L}(x^{k},z^{k},\lambda^{k})-\mathcal{L}(x^{k-1},z^{k-1},\lambda^{k-1})=(a)+(b)$
(23a)
$\displaystyle(a)=\big{[}\mathcal{L}(x^{k},z^{k},\lambda^{k})-\mathcal{L}(x^{k},z^{k},\lambda^{k-1})\big{]}$
(23b)
$\displaystyle(b)=\big{[}\mathcal{L}(x^{k},z^{k},\lambda^{k-1})-\mathcal{L}(x^{k-1},z^{k-1},\lambda^{k-1})\big{]}$
(23c)

First, we examine $(a)$. From the definitions of (8) and (12), it follows that
$(a)=\frac{1}{\rho}\left\lVert\Delta\lambda^{k}\right\rVert^{2}.$ Now,
replacing $k$ with $k-1$ in (20), we get $\lambda^{k-1}=-\theta
z^{k-1}-\tau_{z}\Delta z^{k-1}$ and subtracting this equation from (20), we
obtain:

$\Delta\lambda^{k}=-(\theta+\tau_{z})\Delta z^{k}+\tau_{z}\Delta z^{k-1}.$
(24)

Combining these and using the triangle inequality, we obtain:

$\displaystyle(a)$
$\displaystyle=\frac{1}{\rho}\left\lVert\Delta\lambda^{k}\right\rVert^{2}$
(25) $\displaystyle\leq\frac{1}{\rho}\left((\theta+\tau_{z})\left\lVert\Delta
z^{k}\right\rVert+\tau_{z}\left\lVert\Delta z^{k-1}\right\rVert\right)^{2}$
$\displaystyle\leq\frac{(\theta+\tau_{z})^{2}}{\rho}\left(\left\lVert\Delta
z^{k}\right\rVert+\left\lVert\Delta z^{k-1}\right\rVert\right)^{2}$
$\displaystyle\leq\frac{2(\theta+\tau_{z})^{2}}{\rho}\left(\left\lVert\Delta
z^{k}\right\rVert^{2}+\left\lVert\Delta z^{k-1}\right\rVert^{2}\right),$

where the second and third inequalities follow from $\theta>0$ and
$(a_{1}+a_{2})^{2}\leq 2(a_{1}^{2}+a_{2}^{2})$ for arbitrary real
$a_{1},a_{2}$, respectively.

We now examine $(b)$. For arbitrary $t\in[T]$ and $j,k,l\geq 0$, let
$(x_{<t}^{j},x_{t}^{k},x_{>t}^{l})$ denote the vector obtained by stacking all
$x_{s}^{j}$ for $s<t$, $x_{s}^{k}$ for $s=t$, and $x_{s}^{l}$ for $s>t$. The
term $(b)$ can now be expressed as the following telescoping sum. Note that
since $\lambda$ is fixed at $\lambda^{k-1}$, we temporarily define
$L(x,z)\coloneqq\mathcal{L}(x,z,\lambda^{k-1})$ to simplify exposition.

$\displaystyle(b)$
$\displaystyle=\big{[}L(x^{k},z^{k})-L(x^{k},z^{k-1})\big{]}$
$\displaystyle+\sum_{t=1}^{T}\big{[}L(x_{<t}^{k},x_{t}^{k},x_{>t}^{k-1},z^{k-1})-L(x_{<t}^{k},x_{t}^{k-1},x_{>t}^{k-1},z^{k-1})\big{]}$
$\displaystyle=\alpha+\sum_{t=1}^{T}(\beta_{t}+\gamma_{t})$ (26)

where the expressions for $\beta_{t}$ and $\gamma_{t}$ can be verified by
expanding out both sides of each summand:

$\displaystyle\alpha$ $\displaystyle=L(x^{k},z^{k})-L(x^{k},z^{k-1})$
$\displaystyle\leq-\frac{\tau_{z}}{2}\left\lVert\Delta z^{k}\right\rVert^{2}$
$\displaystyle\beta_{t}$
$\displaystyle=L(x_{<t}^{k-1},x_{t}^{k},x_{>t}^{k-1},z^{k-1})-L(x_{<t}^{k-1},x_{t}^{k-1},x_{>t}^{k-1},z^{k-1})$
$\displaystyle\leq-\frac{\tau_{x}}{2}\left\lVert\Delta
x^{k}_{t}\right\rVert^{2}_{A_{t}^{\top}A_{t}}$ $\displaystyle\gamma_{t}$
$\displaystyle=\left[L(x_{<t}^{k},x_{t}^{k},x_{>t}^{k-1},z^{k-1})-L(x_{<t}^{k-1},x_{t}^{k},x_{>t}^{k-1},z^{k-1})\right]$
$\displaystyle\phantom{=}-\left[L(x_{<t}^{k},x_{t}^{k-1},x_{>t}^{k-1},z^{k-1})-L(x_{<t}^{k-1},x_{t}^{k-1},x_{>t}^{k-1},z^{k-1})\right]$
$\displaystyle=\rho\sum_{s=1}^{t-1}(A_{t}\Delta x_{t}^{k})^{\top}(A_{s}\Delta
x_{s}^{k})$ $\displaystyle\leq\frac{\rho}{2}\sum_{s=1}^{t-1}\left(\left\lVert
A_{t}\Delta x_{t}^{k}\right\rVert^{2}+\left\lVert A_{s}\Delta
x_{s}^{k}\right\rVert^{2}\right)$ $\displaystyle\sum_{t=1}^{T}\gamma_{t}$
$\displaystyle\leq\frac{\rho(T-1)}{2}\sum_{t=1}^{T}\left\lVert\Delta
x_{t}^{k}\right\rVert^{2}_{A_{t}^{\top}A_{t}}$
$\displaystyle\leq\frac{\rho(T-1)}{2}\sum_{t=1}^{T}\left(\left\lVert\Delta
x_{t}^{k}\right\rVert^{2}_{A_{t}^{\top}A_{t}}+\left\lVert\Delta
x_{t}^{k-1}\right\rVert^{2}_{A_{t}^{\top}A_{t}}\right).$

The $\alpha$-inequality is derived as follows: substituting (20) in the
$\lambda$-update formula (12) yields

$\lambda^{k-1}+\rho[Ax^{k}+z^{k}-b]+\theta z^{k}+\tau_{z}\Delta z^{k}=0,$

which is equivalent to $\nabla F(z^{k})=0$, where
$F:z\mapsto\mathcal{L}(x^{k},z,\lambda^{k-1})+\frac{\tau_{z}}{2}\left\lVert
z-z^{k-1}\right\rVert^{2}$ is a convex quadratic function; that is, $z^{k}$
minimizes $F$ and hence, $F(z^{k})\leq F(z^{k-1})$ and this is seen to be
equivalent to the $\alpha$-inequality. The $\beta_{t}$-inequality is obtained
by exploiting the fact $x_{t}^{k}$ is locally optimal with a better objective
value than $x_{t}^{k-1}$ in problem (10). Finally, the $\gamma_{t}$-inequality
follows from $a_{1}a_{2}\leq\frac{1}{2}(a_{1}^{2}+a_{2}^{2})$ for arbitrary
real $a_{1},a_{2}$. The inequality (21) is then obtained by combining (22),
(23), (25) and (26). ∎

######  Lemma 3.

For all $k\geq 1$, we have ($\hat{\Phi}$ defined in (15)):

$\Phi^{k-1}\geq\Phi^{k}\geq\hat{\Phi}.$

###### Proof.

Since $\eta_{x},\eta_{z}>0$, Lemma 2 implies that $\Phi^{k}-\Phi^{k-1}\leq 0$
for all $k\geq 1$. Now suppose, for the sake of contradiction, that there
exists $l\geq 1$ such that $\Phi^{l}<\hat{\Phi}$. Therefore, we must have

$\Phi^{k}\leq\Phi^{l}<\hat{\Phi},\;\forall k\geq l.$ (27)

Consider the augmented Lagrangian (8) and Lyapunov functions (13). Here, all
terms except $\sum_{t=1}^{T}f_{t}(x_{t})$ and $\lambda^{\top}(Ax+z-b)$ are
non-negative. Therefore, using (15) and the $\lambda$-update formula (12), we
have for all $j\geq 1$:

$\displaystyle\Phi^{j}$
$\displaystyle\geq\hat{\Phi}+(\lambda^{j})^{\top}\left(Ax^{j}+z^{j}-b\right)$
$\displaystyle\geq\hat{\Phi}+\frac{1}{\rho}(\lambda^{j})^{\top}(\lambda^{j}-\lambda^{j-1})$
$\displaystyle\geq\hat{\Phi}+\frac{1}{2\rho}\left(\left\lVert\lambda^{j}\right\rVert^{2}-\left\lVert\lambda^{j-1}\right\rVert^{2}\right),$

where the inequality is obtained by noting that

$2a_{1}^{\top}a_{2}=\left\lVert a_{1}\right\rVert^{2}+\left\lVert
a_{2}\right\rVert^{2}-\left\lVert a_{1}-a_{2}\right\rVert^{2}\geq\left\lVert
a_{1}\right\rVert^{2}-\left\lVert a_{1}-a_{2}\right\rVert^{2},$

for arbitrary vectors $a_{1},a_{2}$ and by setting $a_{1}=\lambda^{j}$ and
$a_{2}=\lambda^{j}-\lambda^{j-1}$. Re-arranging the above inequality and
summing over $j\in[k]$, we obtain:

$\sum_{j=1}^{k}\big{(}\Phi^{j}-\hat{\Phi}\big{)}\geq-\frac{1}{2\rho}\left\lVert\lambda^{0}\right\rVert^{2},\;\forall
k\geq 1.$

However, (27) implies:

$\sum_{j=1}^{k}\big{(}\Phi^{j}-\hat{\Phi}\big{)}\leq\sum_{j=1}^{l}\big{(}\Phi^{j}-\hat{\Phi}\big{)}+(k-l)\big{(}\Phi^{l}-\hat{\Phi}\big{)},\;\forall
k\geq l.$

Since $\Phi^{l}-\hat{\Phi}<0$ (by hypothesis), the right-hand side of the
above inequality can be made arbitrarily smaller than 0 for sufficiently large
$k$, contradicting the previous inequality that
$\sum_{j=1}^{k}\big{(}\Phi^{j}-\hat{\Phi}\big{)}\geq-(2\rho)^{-1}\left\lVert\lambda^{0}\right\rVert^{2}$
for all $k\geq 1$. ∎

######  Lemma 4.

For all $K\geq 1$, there exists $j\in[K]$ such that

$\displaystyle\left\lVert\Delta z^{j}\right\rVert^{2}+\left\lVert\Delta
z^{j-1}\right\rVert^{2}\leq\left(\Phi^{1}-\Phi^{K}\right)\left(K\eta_{z}\right)^{-1},$
(28) $\displaystyle\sum_{t=1}^{T}\left\lVert\Delta
x_{t}^{j}\right\rVert_{A_{t}^{\top}A_{t}}^{2}\leq\left(\Phi^{1}-\Phi^{K}\right)\left(K\eta_{x}\right)^{-1}.$
(29)

###### Proof.

Let $j$ be the index that maximizes the right-hand side of (21) in Lemma 2
over $k\in[K]$. The first inequality follows by summing (21) over $k\in[K]$
and noting that $\eta_{x},\eta_{z}>0$. The second one follows similarly if we
ignore the contribution of $\|\Delta x_{t}^{j-1}\|_{A_{t}^{\top}A_{t}}\geq 0$.
∎

######  Proposition 1.

For all $K\geq 1$, there exists $j\in[K]$ such that the primal residual
$\pi^{j}\coloneqq\pi(x^{j})$, defined in (7a), satisfies

$\pi^{j}\leq\sqrt{\frac{2\big{(}\Phi^{1}-\hat{\Phi}\big{)}}{\theta}\left(1+\frac{2(\theta+\tau_{z})^{2}}{K\eta_{z}\rho}\right)}$

###### Proof.

First, we note that

$\Phi^{1}\geq\Phi^{k}\geq\mathcal{L}(x^{k},z^{k},\lambda^{k})\geq\hat{\Phi}+(c),$
(30)

where $(c)$ is equal to

$\displaystyle\frac{\theta}{2}\left\lVert
z^{k}\right\rVert^{2}+(\lambda^{k})^{\top}\left[Ax^{k}+z^{k}-b\right]+\frac{\rho}{2}\left\lVert
Ax^{k}+z^{k}-b\right\rVert^{2}.$

Note that the first inequality in (30) follows from Lemma 3; the second
follows from the definition (14) of $\Phi^{k}$; and the third follows from the
definitions (8) and (15) of the augmented Lagrangian function and
$\hat{\Phi}$, respectively. Substituting $\lambda^{k}$ from (20) in the
expression for $(c)$ above, we obtain

$\displaystyle(c)$ $\displaystyle=\frac{\theta}{2}\left\lVert
z^{k}\right\rVert^{2}-(\theta
z^{k})^{\top}\left[Ax^{k}+z^{k}-b\right]+\frac{\rho}{2}\left\lVert
Ax^{k}+z^{k}-b\right\rVert^{2}$ $\displaystyle\phantom{=}\;-(\tau_{z}\Delta
z^{k})^{\top}\left[Ax^{k}+z^{k}-b\right]$
$\displaystyle=\frac{\theta}{2}\left\lVert
Ax^{k}-b\right\rVert^{2}+\frac{\rho-\theta}{2}\left\lVert
Ax^{k}+z^{k}-b\right\rVert^{2}$ $\displaystyle\phantom{=}\;-(\tau_{z}\Delta
z^{k})^{\top}\left[Ax^{k}+z^{k}-b\right]$
$\displaystyle\geq\frac{\theta}{2}\big{(}\pi^{k}\big{)}^{2}-(\tau_{z}\Delta
z^{k})^{\top}\left[Ax^{k}+z^{k}-b\right]$
$\displaystyle=\frac{\theta}{2}\big{(}\pi^{k}\big{)}^{2}-(\tau_{z}\Delta
z^{k})^{\top}\big{(}\frac{1}{\rho}\Delta\lambda^{k}\big{)},$ (31)

where the first inequality follows from $\rho>\theta$ which can be inferred by
noting that $\eta_{z}>0$ implies (after completing the square in $\tau_{z}$)
that $(16\tau_{z}+16\theta-\rho)^{2}+(16\theta)^{2}-(16\theta-\rho)^{2}<0,$
implying $(16\theta)^{2}-(16\theta-\rho)^{2}<0$, and hence,
$\rho(\rho-32\theta)>0$, that is, $\rho>32\theta>\theta$. Equation (31) then
follows directly from the $\lambda$-update formula (12). The second term in
the right-hand side of (31) can be bounded using: _(i)_ the Cauchy-Schwarz
inequality, _(ii)_ $a_{1}a_{2}\leq\frac{1}{2}(a_{1}^{2}+a_{2}^{2})$ for real
$a_{1},a_{2}$, and _(iii)_ relation (25) from Lemma 2, as follows:

$\displaystyle(\tau_{z}\Delta
z^{k})^{\top}\big{(}\frac{1}{\rho}\Delta\lambda^{k}\big{)}$
$\displaystyle\leq\frac{1}{2\rho}\Big{(}\tau_{z}^{2}\left\lVert\Delta
z^{k}\right\rVert^{2}+\left\lVert\Delta\lambda^{k}\right\rVert^{2}\Big{)}$
$\displaystyle\leq\frac{(\theta+\tau_{z})^{2}}{\rho}\Big{(}\left\lVert\Delta
z^{k}\right\rVert^{2}+\left(\left\lVert\Delta
z^{k}\right\rVert^{2}+\left\lVert\Delta
z^{k-1}\right\rVert^{2}\right)\Big{)}.$ (32)

The claim now follows from inequalities (30), (31), (32), and from the $\Delta
z^{k}$-bound (28) established in Lemma 4. ∎

######  Proposition 2.

For all $K\geq 1$, there exists $j\in[K]$ such that the dual residual
$\delta_{t}^{j}\coloneqq\delta_{t}(x_{t}^{j},\lambda^{j})$, defined in (7b),
satisfies

$\delta_{t}^{j}\leq(\rho+\tau_{x})\left\lVert
A_{t}\right\rVert\sqrt{\frac{(T+1)\big{(}\Phi^{1}-\Phi^{K}\big{)}}{K\min\\{\eta_{x},\eta_{z}\\}}}$
(33)

for all $t\in[T]$.

###### Proof.

Fix $t\in[T]$. Since $x_{t}^{k}$ is locally optimal in problem (10), the
first-order optimality conditions imply, for $k\geq 1$:

$\left(\begin{aligned} \nabla
f_{t}(x_{t}^{k})+A_{t}^{\top}\lambda^{k-1}+\tau_{x}A_{t}^{\top}A_{t}\Delta
x_{t}^{k}\\\ +\rho A_{t}^{\top}\left[A_{t}x_{t}^{k}+A_{\neq t}x_{\neq
t}^{k-1}+z^{k-1}-b\right]\end{aligned}\right)\in-N_{X_{t}}(x_{t}^{k}).$

Substituting $\lambda^{k-1}=\lambda^{k}-\rho[Ax^{k}+z^{k}-b]$ from (12):

$\left(\begin{aligned} \nabla
f_{t}(x_{t}^{k})+A_{t}^{\top}\lambda^{k}+\tau_{x}A_{t}^{\top}A_{t}\Delta
x_{t}^{k}\\\ -\rho A_{t}^{\top}A_{\neq t}\Delta x_{\neq t}^{k}-\rho
A_{t}^{\top}\Delta z^{k}\end{aligned}\right)\in-N_{X_{t}}(x_{t}^{k}).$

Definition (7b) of $\delta_{t}^{k}=\mathop{\mathrm{dist}}\left(\nabla
f_{t}(x_{t}^{k})+A_{t}^{\top}\lambda^{k},-N_{X_{t}}(x_{t}^{k})\right)$ implies
that the latter quantity is bounded from above by:

$\displaystyle(\delta_{t}^{k})^{2}\leq\left\lVert-\rho A_{t}^{\top}A_{\neq
t}\Delta x_{\neq t}^{k}-\rho A_{t}^{\top}\Delta
z^{k}+\tau_{x}A_{t}^{\top}A_{t}\Delta x_{t}^{k}\right\rVert^{2}$
$\displaystyle\leq\left\lVert A_{t}\right\rVert^{2}\left(\rho\left\lVert
A_{\neq t}\Delta x_{\neq t}^{k}\right\rVert+\rho\left\lVert\Delta
z^{k}\right\rVert+\tau_{x}\left\lVert A_{t}\Delta
x_{t}^{k}\right\rVert\right)^{2}$
$\displaystyle\leq(T+1)(\rho+\tau_{x})^{2}\left\lVert
A_{t}\right\rVert^{2}\Big{(}\sum_{s=1}^{T}\left\lVert\Delta
x_{s}^{k}\right\rVert_{A_{s}^{\top}A_{s}}^{2}+\left\lVert\Delta
z^{k}\right\rVert^{2}\Big{)}.$

where the inequality on the second line follows from the Cauchy-Schwarz and
triangle inequalities, and the last inequality follows from
$\max\\{\rho,\tau_{x}\\}\leq(\rho+\tau_{x})$ and $(a_{1}+\ldots+a_{N})^{2}\leq
N(a_{1}^{2}+\ldots a_{N}^{2})$ for arbitrary real $a_{1},\ldots,a_{N}$. The
claim now follows from the bounds (28), (29) in Lemma 4. ∎

##  Appendix B Local convergence proof

Throughout this section, we suppose that the conditions outlined in Theorem 2
hold. Also, for $k\geq 1$, we define
$\mu^{k}\coloneqq(\mu_{1}^{k},\ldots,\mu_{T}^{k})$, where $\mu_{t}^{k}$ is the
optimal Lagrange multiplier vector of the $x_{t}$-subproblem (10) at iteration
$k$.

The key steps of the proof are as follows. We first introduce the following
perturbed variant of problem (2) that is parameterized by $p\in\mathbb{R}^{m}$
and $d=(d_{1},\ldots,d_{T},d_{z})\in\mathbb{R}^{n+m}$.

$\displaystyle\mathop{\text{minimize}}_{x_{1},\ldots,x_{T},z}$
$\displaystyle\sum_{t=1}^{T}\left[f_{t}(x_{t})-d_{t}^{\top}x_{t}\right]+\frac{\theta}{2}\left\lVert
z\right\rVert^{2}-d_{z}^{\top}z$ ($S(p,d)$) subject to $\displaystyle
c_{t}(x_{t})=0,\;\;t\in[T],$ $\displaystyle Ax+z=b+p.$

Observe that $S(0,0)$ coincides precisely with problem (2). By defining
$p^{k}$ and $d^{k}$ to be the primal and dual residuals of problem (2) at
$(x^{k},z^{k})$ as per eq. (19), Lemma 5 and 6 show that
$(x^{k},z^{k},\mu^{k},\lambda^{k})$ must coincide precisely with the optimal
solution of $S(p^{k},d^{k})$, whenever the former is close to
$(x^{*},z^{*},\mu^{*},\lambda^{*})$ and $p^{k}$, $d^{k}$ are sufficiently
close to $0$. This allows us to bound the difference in the augmented
Lagrangian function values of problem (2) evaluated at its optimal solution
and at an arbitrary iterate (Lemma 7 and 8). In Lemma 9, we establish that it
is sufficient to have $(Dx^{k},z^{k},\lambda^{k})$ to be close to
$(Dx^{*},z^{*},\lambda^{*})$. Finally, Proposition 3 proves that whenever the
latter is true, then the sequence $\\{(x^{k},z^{k})\\}$ must converge to
$(x^{*},z^{*})$. This requires some intermediate results that are proved in
Lemma 10 and 11. The proof of Theorem 2 follows directly from that of
Proposition 3. Indeed, its first part is already true due to Theorem 1,
whereas its second part follows from Assumption (A7) if we set $w_{z}=0$
therein.

######  Lemma 5.

There exist $C^{1}$ functions $\hat{x}$, $\hat{z}$, $\hat{\mu}$,
$\hat{\lambda}$ with domain $\mathbb{R}^{m}\times\mathbb{R}^{n+m}$ and
constants $\epsilon_{1},\epsilon_{2}>0$ such that

  1. 1.

$(\hat{x}(0,0),\hat{z}(0,0),\hat{\mu}(0,0),\hat{\lambda}(0,0))=(x^{*},z^{*},\mu^{*},\lambda^{*})$.

  2. 2.

For all $(p,d)\in B_{\epsilon_{1}}(0)$, $(\hat{x}(p,d),\hat{z}(p,d))$ is a
local minimizer of $S(p,d)$ where the LICQ and SOSC assumptions are satisfied
with optimal Lagrange multipliers $(\hat{\mu}(p,d),\hat{\lambda}(p,d))$.

  3. 3.

The functions $\hat{x}$, $\hat{z}$, $\hat{\mu}$, $\hat{\lambda}$ map to
locally unique points in the sense that, for all $(p,d)\in
B_{\epsilon_{1}}(0)$, if $(\bar{x},\bar{z})$ is a local minimizer (or a first-
order stationary point) of $S(p,d)$ with optimal Lagrange multipliers
$(\bar{\mu},\bar{\lambda})$, and if
$(\bar{x},\bar{z},\bar{\lambda},\bar{\mu})\in
B_{\epsilon_{2}}((x^{*},z^{*},\mu^{*},\lambda^{*}))$, then
$(\bar{x},\bar{z},\bar{\lambda},\bar{\mu})=(\hat{x}(p,d),\hat{z}(p,d),\hat{\mu}(p,d),\hat{\lambda}(p,d))$.

###### Proof.

The LICQ (A6) and SOSC assumptions (A7) are satisfied at $(x^{*},z^{*})$, The
claims then follow from a direct application of [42, Theorem 2.1] to $S(p,d)$.
∎

######  Lemma 6.

For all $k\geq 1$ such that $(p^{k},d^{k})\in B_{\epsilon_{1}}(0)$ and
$(x^{k},z^{k},\mu^{k},\lambda^{k})\in
B_{\epsilon_{2}}((x^{*},z^{*},\mu^{*},\lambda^{*}))$, where
$\epsilon_{1},\epsilon_{2}$ are defined in Lemma 5 and $p^{k}$ and $d^{k}$
from (19), we have:

$(\hat{x}(p^{k},d^{k}),\hat{z}(p^{k},d^{k}),\hat{\mu}(p^{k},d^{k}),\hat{\lambda}(p^{k},d^{k}))=(x^{k},z^{k},\mu^{k},\lambda^{k}).$

###### Proof.

We first show that $(x^{k},z^{k},\mu^{k},z^{k})$ satisfy the KKT conditions of
$S(p^{k},d^{k})$. To see this, first note that $x_{t}^{k}$ is locally optimal
in problem (10). Under assumption (A8), this implies that it satisfies the
first-order optimality condition:

$\displaystyle\left(\begin{aligned} \nabla f_{t}(x_{t}^{k})+\nabla
c_{t}(x_{t}^{k})\mu_{t}^{k}+A_{t}^{\top}\lambda^{k-1}+\tau_{x}A_{t}^{\top}A_{t}\Delta
x_{t}^{k}\\\ +\rho A_{t}^{\top}\left[A_{t}x_{t}^{k}+A_{\neq t}x_{\neq
t}^{k-1}+z^{k-1}-b\right]\end{aligned}\right)=0.$

Using the $\lambda$-update (12) and eq. (19b), this is equivalent to:

$\displaystyle\nabla f_{t}(x_{t}^{k})+\nabla
c_{t}(x_{t}^{k})\mu_{t}^{k}+A_{t}^{\top}\lambda^{k}-d_{t}^{k}=0.$

Second, note that Lemma 1 along with equation (19c) implies:

$\displaystyle\theta z^{k}+\lambda^{k}=d_{z}^{k}.$

Finally, equation (19a) can be equivalently written as follows:

$\displaystyle Ax^{k}+z^{k}=b+p^{k}.$

Along with $c_{t}(x_{t}^{k})=0$, the last three equations are precisely the
KKT conditions of $S(p^{k},d^{k})$. The result now follows from the third part
of Lemma 5. ∎

######  Lemma 7.

There exist constants $\rho_{1},\epsilon_{3}>0$ such that for all
$\rho>\rho_{1}$ and $k\geq 1$ where $(x^{k},z^{k})\in
B_{\epsilon_{3}}((x^{*},z^{*}))$, we have:

$\displaystyle\sum_{t\in[T]}f_{t}(x_{t}^{k})+\frac{\theta}{2}\left\lVert
z^{k}\right\rVert^{2}+(\lambda^{*})^{\top}p^{k}+\frac{\rho}{4}\left\lVert
p^{k}\right\rVert^{2}$ (34)
$\displaystyle\;\;\geq\sum_{t\in[T]}f_{t}(x_{t}^{*})+\frac{\theta}{2}\left\lVert
z^{*}\right\rVert^{2}.$

###### Proof.

Define the (fully) augmented Lagrangian function of $S(0,0)$ with respect to
the optimal multipliers $(\mu^{*},\lambda^{*})$, as follows
$\mathcal{S}:\mathbb{R}^{n}\times\mathbb{R}^{m}\mapsto\mathbb{R}$, where
$\mathcal{S}(x,z)$ is given by:

$\Lambda_{\theta}(x,z,\mu^{*},\lambda^{*})+\frac{\rho}{4}\sum_{t=1}^{T}\left\lVert
c_{t}(x_{t})\right\rVert^{2}+\frac{\rho}{4}\left\lVert Ax+z-b\right\rVert^{2}$

The definition of $(x^{*},z^{*})$ as a local minimizer of problem $S(0,0)$
satisfying the LICQ assumption, means that $\nabla\mathcal{S}(x^{*},z^{*})=0$
and $c_{t}(x^{*}_{t})=0$. The latter also implies that the Hessian
$\nabla^{2}\mathcal{S}(x^{*},z^{*})$ is given by:

$\begin{bmatrix}\nabla^{2}_{xx}\Lambda(x^{*},z^{*},\mu^{*},\lambda^{*})&0\\\
0&\theta\mathrm{I}\end{bmatrix}+\frac{\rho}{2}\begin{bmatrix}J(x^{*})^{\top}J(x^{*})&A^{\top}\\\
A&\mathrm{I}\end{bmatrix}.$

Assumption (A7) and [43, Lemma 3.2.1] ensure the existence of $\rho_{1}>0$
such that $\nabla^{2}\mathcal{S}(x^{*},z^{*})$ is positive definite and that
its minimum eigenvalue is larger than some fixed $\epsilon_{3}^{\prime}>0$,
for all $\rho>\rho_{1}$. Continuity of $\nabla^{2}\mathcal{S}$ along with the
second-order sufficient conditions for unconstrained minimization then imply
the existence of $\epsilon_{3}>0$ such that $(x^{*},z^{*})$ is a local
minimizer of $\mathcal{S}$ in some neighborhood of radius $\epsilon_{3}$
around $(x^{*},z^{*})$. Note that the radius $\epsilon_{3}$ is independent of
$\rho$ as long as the minimum eigenvalue of
$\nabla^{2}\mathcal{S}(x^{*},z^{*})$ remains larger than
$\epsilon_{3}^{\prime}>0$. The statement of the lemma then follows by noting
that $\mathcal{S}(x^{k},z^{k})\geq\mathcal{S}(x^{*},z^{*})$ and after
substituting the expressions for $\mathcal{S}(x,z)$ and for $p^{k}$ from
(19a). ∎

######  Lemma 8.

There exist constants $\rho_{2},\epsilon_{4}>0$ such that for all
$\rho>\rho_{2}$ and $k\geq 1$ where $(p^{k},d^{k})\in B_{\epsilon_{4}}(0)$ and
$(x^{k},z^{k},\mu^{k},\lambda^{k})\in
B_{\epsilon_{2}}((x^{*},z^{*},\mu^{*},\lambda^{*}))$, with $\epsilon_{2}$
defined in Lemma 5, we have:

$\displaystyle\sum_{t\in[T]}\big{(}f_{t}(x_{t}^{k})-d_{t}^{k}x_{t}^{k}\big{)}+\frac{\theta}{2}\left\lVert
z^{k}\right\rVert^{2}-(d^{k}_{z})^{\top}z^{k}+(\lambda^{k})^{\top}p^{k}-\frac{\rho}{4}\left\lVert
p^{k}\right\rVert^{2}$ (35)
$\displaystyle\;\;\leq\sum_{t\in[T]}\big{(}f_{t}(x_{t}^{*})-d_{t}^{k}x_{t}^{*}\big{)}+\frac{\theta}{2}\left\lVert
z^{*}\right\rVert^{2}-(d^{k}_{z})^{\top}z^{*}.$

###### Proof.

Define the primal functional corresponding to $S(p,d)$, as follows
$\mathcal{Q}:\mathbb{R}^{m}\times\mathbb{R}^{n+m}\mapsto\mathbb{R}$, where
$\mathcal{Q}(p,d)$ is given by:

$\displaystyle\sum_{t=1}^{T}\left[f_{t}(\hat{x}_{t}(p,d))-d_{t}^{\top}\hat{x}_{t}(p,d)\right]+\frac{\theta}{2}\left\lVert\hat{z}(p,d)\right\rVert^{2}-d_{z}^{\top}\hat{z}(p,d)\;.$

Now [5, Proposition 1.28] implies that $\nabla_{p}Q(p,d)=-\hat{\lambda}(p,d)$,
which is continuously differentiable from Lemma 5. Therefore,
$\nabla^{2}_{pp}Q$ is continuous and [43, Lemma 3.2.1] ensures the existence
of $\rho_{2}>0$ such that $\nabla^{2}_{pp}Q(0,0)+\frac{\rho}{2}\mathrm{I}$ is
positive definite and that its minimum eigenvalue is bounded strictly away
from $0$ for all $\rho>\rho_{2}$. By a similar argument as in the proof of
Lemma 7, continuity of $\nabla^{2}_{pp}Q$ implies that there exists
$\epsilon_{4}^{\prime}>0$ (independent of $\rho$) such that
$\nabla^{2}_{pp}Q(p,d)+\frac{\rho}{2}\mathrm{I}$ is positive definite, and
hence the function $F_{d}:p\mapsto Q(p,d)+\frac{\rho}{4}\left\lVert
p\right\rVert^{2}$ is convex, whenever $(p,d)\in
B_{\epsilon_{4}^{\prime}}(p,d)$. Now fix
$\epsilon_{4}=\min\\{\epsilon_{4}^{\prime},\epsilon_{1}\\}>0$ and $d=d^{k}$.
Convexity of $F_{d_{k}}$ implies:

$F_{d_{k}}(p^{k})+(0-p^{k})^{\top}\nabla F_{d_{k}}(p^{k})\leq F_{d_{k}}(0).$

The statement of the lemma now follows by _(i)_ substituting
$F_{d_{k}}(p^{k})=Q(p^{k},d^{k})+\frac{\rho}{4}\left\lVert
p^{k}\right\rVert^{2}$, replacing $Q(p^{k},d^{k})$ using its definition, and
noting $(\hat{x}(p^{k},d^{k}),\hat{z}(p^{k},d^{k}))=(x^{k},z^{k})$ from Lemma
6; _(ii)_ substituting $\nabla
F_{d_{k}}(p^{k})=-\hat{\lambda}(p^{k},d^{k})+\frac{\rho}{2}p^{k}=-\lambda^{k}+\frac{\rho}{2}p^{k}$
(from Lemma 6); and, _(iii)_ $F_{d_{k}}(0)=Q(0,d^{k})$ is less than or equal
to the right-hand side of (35), since $(x^{*},z^{*})$ is feasible (but
possibly suboptimal) in $S(0,d^{k})$. ∎

######  Lemma 9.

There exist $\rho_{3},\epsilon_{5}>0$ such that for all $\rho>\rho_{3}$ and
$k\geq 1$ where $(Dx^{k-1},z^{k-1},\lambda^{k-1})\in
B_{\epsilon_{5}}((Dx^{*},z^{*},\lambda^{*}))$, we have
$(x^{k},z^{k},\mu^{k},\lambda^{k})\in
B_{\min\\{\epsilon_{2},\epsilon_{3}\\}}((x^{*},z^{*},\mu^{*},\lambda^{*})),$
where $\epsilon_{2},\epsilon_{3}$ are defined in Lemma 5 and 7.

###### Proof.

Observe that for fixed choices of $\rho,\theta,\tau_{x},\tau_{z}$ in Algorithm
1, the $x$-update step (10) at iteration $k$ depends only on the values of
$Dx^{k-1}$, $z^{k-1}$ and $\lambda^{k-1}$ at the previous iteration;
therefore, we can define
$G:(Dx^{k-1},z^{k-1},\lambda^{k-1})\mapsto(x^{k},\mu^{k})$ to be the
corresponding mapping. Observe now that the claim follows trivially if we can
show that $G$ is continuous and satisfies
$G(Dx^{*},z^{*},\lambda^{*})=(x^{*},\mu^{*})$. Indeed, if this is true, then
observe that the $z$-update (11) and $\lambda$-update (12) formulas also have
the same properties: they define continuous maps and satisfy $z^{k}=z^{*}$ and
$\lambda^{k}=\lambda^{*}$ whenever $x^{k}=x^{*}$, $z^{k-1}=z^{*}$ and
$\lambda^{k-1}=\lambda^{*}$, since $Ax^{*}+z^{*}=b$ and $\lambda^{*}=-\theta
z^{*}$ (from the KKT conditions of problem (2)).

Now fix $t\in[T]$ and
$(Dx^{k-1},z^{k-1},\lambda^{k-1})=(Dx^{*},z^{*},\lambda^{*})$. Observe that
problem (10) satisfies the LICQ assumption at $x_{t}^{*}$, because of (A6).
Moreover, since $Ax^{*}+z^{*}=b$, observe that $(x_{t}^{*},\mu_{t}^{*})$ is
also a first-order stationary point of problem (10):

$\displaystyle\left(\begin{aligned} \nabla f_{t}(x_{t}^{*})+\nabla
c_{t}(x_{t}^{*})\mu_{t}^{*}+A_{t}^{\top}\lambda^{*}+\tau_{x}A_{t}^{\top}A_{t}(x_{t}^{*}-x_{t}^{*})\\\
+\rho A_{t}^{\top}\left[A_{t}x_{t}^{*}+A_{\neq t}x_{\neq
t}^{*}+z^{*}-b\right]\end{aligned}\right)=0.$

To show that it satisfies SOSC, note that Assumption (A7) and [43, Lemma
3.2.1] ensure the existence of $\rho_{3}>0$ such that

$\begin{bmatrix}\nabla^{2}_{xx}\Lambda(x^{*},z^{*},\mu^{*},\lambda^{*})&0\\\
0&\theta\mathrm{I}\end{bmatrix}+\rho\begin{bmatrix}A^{\top}A&A^{\top}\\\
A&\mathrm{I}\end{bmatrix}.$

is positive definite on the domain
$\\{(w_{x},w_{z})\in\mathbb{R}^{n}\times\mathbb{R}^{m}\setminus\\{0\\}:\nabla
c(x^{*})^{\top}w_{x}=0\\}$ for all $\rho>\rho_{3}$. This means that for
$w_{x_{t}}\neq 0$, $\nabla c_{t}(x_{t}^{*})^{\top}w_{x_{t}}=0$, $w_{z}=0$, and
$w_{x_{s}}=0$ for $s\neq t$, we have:

$w_{x_{t}}^{\top}\left[\nabla^{2}_{x_{t}x_{t}}\Lambda(x^{*},z^{*},\mu^{*},\lambda^{*})+(\rho+\tau_{x})A_{t}^{\top}A_{t}\right]w_{x_{t}}>0,$

which is precisely the SOSC condition for problem (10) at
$(x_{t}^{*},\mu_{t}^{*})$, and therefore, $(x^{*},\mu^{*})$ is a strict local
solution.

Continuity of the mapping $G$ now follows directly from Assumption (A9) and
classical NLP sensitivity [42, Theorem 2.1] applied to problem (10) with
$Dx_{t}^{k-1}$, $z^{k-1}$ and $\lambda^{k-1}$ viewed as perturbation
parameters. ∎

######  Lemma 10.

The sequences $\\{p^{k}\\}$ and $\\{d^{k}\\}$ converge to $0$.

###### Proof.

After summing inequality (22) over $k\in[K]$ and noting
$\Phi^{K}\geq\hat{\Phi}$ from Lemma 3, we obtain that the sequences $\\{\Delta
z^{k}\\}$ and $\\{A_{t}\Delta x_{t}^{k}\\}$ for $t\in[T]$ must all converge to
$0$. Therefore, the definitions (19b) and (19c) imply that $\\{d^{k}\\}$ also
converges to $0$. Equations (12) and (24) imply
$p^{k}=\frac{1}{\rho}\Delta\lambda^{k}=-\frac{\theta+\tau_{z}}{\rho}\Delta
z^{k}+\frac{\tau_{z}}{\rho}\Delta z^{k-1}$, which means that the sequence
$\\{p^{k}\\}$ must also converge to $0$. ∎

######  Lemma 11.

  1. 1.

Any two vectors $a_{1},a_{2}$ of equal dimension satisfy
$2a_{1}^{\top}a_{2}=\left\lVert a_{1}+a_{2}\right\rVert^{2}-\left\lVert
a_{1}\right\rVert^{2}-\left\lVert a_{2}\right\rVert^{2}$.

  2. 2.

For any
$w=(w_{1},\ldots,w_{T})\in\mathbb{R}^{n_{1}}\times\ldots\times\mathbb{R}^{n_{T}}$,
we have:

$\displaystyle(\rho+\tau_{x})\sum_{t=1}^{T}\left\lVert
w_{t}\right\rVert^{2}_{A_{t}^{\top}A_{t}}-\rho\left\lVert
w\right\rVert^{2}_{A^{\top}A}=\left\lVert Dw\right\rVert^{2}_{R},$
$\displaystyle R\coloneqq(\rho+\tau_{x})\mathrm{I}-\rho EE^{\top}\succ 0,$
$\displaystyle E^{\top}\coloneqq\begin{bmatrix}\mathrm{I}_{m\times
m}&\ldots&\mathrm{I}_{m\times m}\end{bmatrix}\in\mathbb{R}^{m\times Tm}.$

###### Proof.

The first part follows by expanding its right-hand side. To prove the second
part, note that $D=\mathop{\mathrm{diag}}(A_{1},\ldots,A_{T})$ can be used to
verify that $\sum_{t=1}^{T}\left\lVert
w_{t}\right\rVert^{2}_{A_{t}^{\top}A_{t}}=w^{\top}D^{\top}Dw.$ Similarly, the
definition of $E$ can be used to verify that $A=E^{\top}D$ and hence,
$\left\lVert w\right\rVert^{2}_{A^{\top}A}=w^{\top}D^{\top}EE^{\top}Dw.$ This
proves the claimed equation. To show that $R\succ 0$, observe that
$E\;E^{\top}$ can be equivalently expressed as the Kronecker product
$\mathrm{e}\;\mathrm{e}^{\top}\otimes\mathrm{I}_{m\times m}$, where
$\mathrm{e}$ is the vector of ones in $\mathbb{R}^{T}$. Therefore, its
eigenvalues are given by pairwise products of the eigenvalues of
$\mathrm{e}\;\mathrm{e}^{\top}$ (which are $0$ and $T$) and
$\mathrm{I}_{m\times m}$ (which is $1$). Therefore, the eigenvalues of
$E\;E^{\top}$ are $0$ and $T$, and hence, those of
$(\rho+\tau_{x})\mathrm{I}-\rho E\;E^{\top}$ are $\rho+\tau_{x}$ and
$\rho+\tau_{x}-\rho T$, both of which are positive since $\eta_{x}>0$, see
(16). ∎

######  Proposition 3.

There exist constants $\bar{\rho},\epsilon^{\prime}>0$ such that if
$\rho>\bar{\rho}$ and
$\left\lVert(Dx^{0},z^{0},\lambda^{0})-(Dx^{*},z^{*},\lambda^{*})\right\rVert<\epsilon^{\prime}$,
then the sequence $\\{(x^{k},z^{k},\lambda^{k},\Delta z^{k})\\}$ converges to
$(x^{*},z^{*},\lambda^{*},0)$ at a sublinear rate.

###### Proof.

Set $\bar{\rho}=\max\\{\rho_{1},\rho_{2},\rho_{3}\\}$, where the latter are
defined in Lemma 7, 8, 9. Also, assume for the moment that the other
conditions of Lemma 7 and 8 are also satisfied for all $k\geq 1$; we shall
shortly show how $\epsilon^{\prime}$ can be chosen to ensure this. Adding the
corresponding equations (34) and (35), we obtain:

$\begin{aligned}
&\underset{(d)}{\underbrace{(\lambda^{*}-\lambda^{k})^{\top}p^{k}+\frac{\rho}{2}\left\lVert
p^{k}\right\rVert^{2}}}+\underset{(e)}{\underbrace{(d_{z}^{k})^{\top}(z^{k}-z^{*})}}\\\
&\quad+\underset{(f)}{\underbrace{\sum_{t=1}^{T}(d_{t}^{k})^{\top}(x_{t}^{k}-x_{t}^{*})}}\end{aligned}\geq
0.$

We now we have the following

$\displaystyle(d)$
$\displaystyle=\frac{1}{\rho}(\lambda^{*}-\lambda^{k})^{\top}\Delta\lambda^{k}+\frac{1}{2\rho}\left\lVert\Delta\lambda^{k}\right\rVert^{2}$
$\displaystyle=\frac{1}{2\rho}\big{[}\left\lVert\lambda^{*}-\lambda^{k-1}\right\rVert^{2}-\left\lVert\lambda^{*}-\lambda^{k}\right\rVert^{2}\big{]},$

where the first equality follows from definition (19a) and the
$\lambda$-update formula (12), and the second equality follows from the first
part of Lemma 11. Subsequently,

$\displaystyle(e)$ $\displaystyle=\tau_{z}(\Delta z^{k})^{\top}(z^{*}-z^{k})$
$\displaystyle=\frac{\tau_{z}}{2}\big{[}\left\lVert
z^{*}-z^{k-1}\right\rVert^{2}-\left\lVert
z^{*}-z^{k}\right\rVert^{2}-\left\lVert\Delta z^{k}\right\rVert^{2}\big{]},$

where again the first equality follows from definition (19c) and the second
follows from the first part of Lemma 11. Finally,

$\displaystyle(f)$ $\displaystyle=\sum_{t=1}^{T}\big{[}\rho A\Delta
x^{k}+\rho\Delta z^{k}-(\rho+\tau_{x})A_{t}\Delta
x_{t}^{k}\big{]}^{\top}A_{t}(x_{t}^{k}-x_{t}^{*})$
$\displaystyle=(f_{1})+(f_{2}),$

where the first equality follows from definition (19b) and after adding and
subtracting $\rho A_{t}\Delta x_{t}^{k}$ from the latter. For the first
component of $(f)$ we have

$\displaystyle(f_{1})$ $\displaystyle=-\rho(A\Delta
x^{k})^{\top}A(x^{*}-x^{k})$
$\displaystyle\phantom{=}+(\rho+\tau_{x})\sum_{t=1}^{T}(A_{t}\Delta
x_{t}^{k})^{\top}A_{t}(x_{t}^{*}-x_{t}^{k})$ $\displaystyle=$
$\displaystyle-\frac{\rho}{2}\left[\left\lVert
x_{t}^{*}-x_{t}^{k-1}\right\rVert_{A^{\top}A}^{2}-\left\lVert
x_{t}^{*}-x_{t}^{k}\right\rVert_{A^{\top}A}^{2}-\left\lVert\Delta
x_{t}^{k}\right\rVert_{A^{\top}A}^{2}\right]$
$\displaystyle+\frac{\rho+\tau_{x}}{2}\sum_{t=1}^{T}\left[\begin{aligned}
&\left\lVert
x_{t}^{*}-x_{t}^{k-1}\right\rVert^{2}_{A_{t}^{\top}A_{t}}-\left\lVert
x_{t}^{*}-x_{t}^{k}\right\rVert^{2}_{A_{t}^{\top}A_{t}}\\\ &-\left\lVert\Delta
x_{t}^{k}\right\rVert^{2}_{A_{t}^{\top}A_{t}}\end{aligned}\right]$
$\displaystyle=$ $\displaystyle\frac{1}{2}\big{[}\left\lVert
D(x^{*}-x^{k-1})\right\rVert_{R}^{2}-\left\lVert
D(x^{*}-x^{k})\right\rVert_{R}^{2}-\left\lVert D\Delta
x^{k}\right\rVert_{R}^{2}\big{]}$

where the second and third equalities follow from the first and second parts
of Lemma 11, respectively. For the second component of $f$ we have

$\displaystyle(f_{2})$ $\displaystyle=\rho(\Delta z^{k})^{\top}A(x^{k}-x^{*})$
$\displaystyle=(\Delta z^{k})^{\top}\Delta\lambda^{k}+\rho(\Delta
z^{k})^{\top}(z^{*}-z^{k})$ $\displaystyle=-(\theta+\tau_{z})\left\lVert\Delta
z^{k}\right\rVert^{2}+\tau_{z}(\Delta z^{k})^{\top}\Delta z^{k-1}$
$\displaystyle\phantom{=}+\frac{\rho}{2}\big{[}\left\lVert
z^{*}-z^{k-1}\right\rVert^{2}-\left\lVert
z^{*}-z^{k}\right\rVert^{2}-\left\lVert\Delta z^{k}\right\rVert^{2}\big{]}$
$\displaystyle\leq\frac{\rho}{2}\big{[}\left\lVert
z^{*}-z^{k-1}\right\rVert^{2}-\left\lVert
z^{*}-z^{k}\right\rVert^{2}\big{]}+\frac{\tau_{z}}{2}\left\lVert\Delta
z^{k}\right\rVert^{2},$

where the second equality follows by substituting $Ax^{*}=b-z^{*}$ and
$Ax^{k}=b-z^{k}+\frac{1}{\rho}\Delta\lambda^{k}$, the third equality follows
from (24) and the first part of Lemma 11, and the inequality follows by noting
first that $2(\Delta z^{k})^{\top}\Delta z^{k-1}\leq\|{\Delta
z^{k}}\|^{2}+\|{\Delta z^{k-1}}\|^{2}$ and then that
$\theta+(\tau_{z}/2)+\rho>0$.

Combining the equations for $(d),(e),(f),(f_{1})$, the inequality for
$(f_{2})$, and $(d)+(e)+(f)\geq 0$, we obtain:

$\displaystyle\left\lVert(Dx^{k},z^{k},\lambda^{k})-(Dx^{*},z^{*},\lambda^{*})\right\rVert^{2}_{*}+\tau_{z}\left\lVert\Delta
z^{k}\right\rVert^{2}+\left\lVert D\Delta x^{k}\right\rVert^{2}_{R}$ (36)
$\displaystyle\leq\left\lVert(Dx^{k-1},z^{k-1},\lambda^{k-1})-(Dx^{*},z^{*},\lambda^{*})\right\rVert^{2}_{*}+\tau_{z}\left\lVert\Delta
z^{k-1}\right\rVert^{2},$

where we define the norm:

$\left\lVert(Dx,z,\lambda)\right\rVert_{*}\coloneqq\sqrt{\left\lVert
Dx\right\rVert^{2}_{R}+(\rho+\tau_{z})\left\lVert
z\right\rVert^{2}+(1/\rho)\left\lVert\lambda\right\rVert^{2}}.$

Summing inequality (36) over $k\in[j]$ and using $\Delta
z^{0}=(-\lambda^{0}-\theta z^{0})/\tau_{z}=0$ (by hypothesis), we obtain for
all $j\geq 1$:

$\displaystyle\left\lVert(Dx^{j},z^{j},\lambda^{j})-(Dx^{*},z^{*},\lambda^{*})\right\rVert^{2}_{*}+\tau_{z}\left\lVert\Delta
z^{j}\right\rVert^{2}+\left\lVert D\Delta x^{j}\right\rVert^{2}_{R}$
$\displaystyle\leq\left\lVert(Dx^{0},z^{0},\lambda^{0})-(Dx^{*},z^{*},\lambda^{*})\right\rVert^{2}_{*}$

The equivalence of norms implies that whenever the Euclidean norm
$\left\lVert(Dx^{0},z^{0},\lambda^{0})-(Dx^{*},z^{*},\lambda^{*})\right\rVert$
is less than $\epsilon^{\prime}$, the right-hand side of the above inequality
is also sufficiently small, and the terms on the left-hand side are even
smaller. In particular, by choosing a small $\epsilon^{\prime}$, we can
ensure: _(i)_ $(Dx^{j},z^{j},\lambda^{j})$ remains close to
$(Dx^{*},z^{*},\lambda^{*})$; and _(ii)_ $D\Delta x^{j}$ and $\Delta z^{j}$
are close to $0$, which implies $p^{j},d^{j}$ are also close to $0$ (see
argument in proof of Lemma 10). Finally, note that (i) and (ii) satisfy the
conditions of Lemma 9 (which in turn allows us to satisfy those of Lemma 7)
and Lemma 8, respectively.

To be precise, the equivalence of norms implies there exist $0<c_{1}\leq
c_{2}$ such that
$c_{1}\left\lVert\cdot\right\rVert_{*}\leq\left\lVert\cdot\right\rVert\leq
c_{2}\left\lVert\cdot\right\rVert_{*}$. Now set
$\epsilon^{\prime}=\min\\{\frac{c_{1}}{c_{0}}\epsilon_{1},\frac{c_{1}}{c_{2}}\epsilon_{5}\\}<\epsilon_{5}$,
where $\epsilon_{1}$ and $\epsilon_{5}$ are defined in Lemma 5 and 9,
respectively, and $c_{0}\coloneqq
4(\rho+\theta+\tau_{x}+\tau_{z}+1)(\frac{1}{\rho}+1)(\frac{1}{\sqrt{\tau_{z}}}+1)(T+1)(\left\lVert
D\right\rVert+1)(c_{2}+1)$ is sufficiently large. Then, it can be verified
that for all $k\geq 0$, we have: $(Dx^{k},z^{k},\lambda^{k})\in
B_{\epsilon_{5}}(Dx^{*},z^{*},\lambda^{*})$, which verifies the conditions of
Lemma 9 and hence of Lemma 7, and $(p^{k},d^{k})\in B_{\epsilon_{1}}(0)$,
which verifies the condition of Lemma 8. Now, since the conditions of Lemma 6
are also satisfied, and since $\\{p^{k}\\}$ and $\\{d^{k}\\}$ converge to $0$
(from Lemma 10), this proves that
$\\{(x^{k},z^{k},\lambda^{k})\\}=\\{(\hat{x}(p^{k},d^{k}),\hat{z}(p^{k},d^{k}),\hat{\lambda}(p^{k},d^{k}))\\}$
converge to
$(\hat{x}(0,0),\hat{z}(0,0),\hat{\lambda}(0,0))=(x^{*},z^{*},\lambda^{*})$
from Lemma 5. Finally, the convergence rate follows from (36) as follows:

$\frac{\left\lVert(Dx^{k},z^{k},\lambda^{k},\Delta
z^{k})-(Dx^{*},z^{*},\lambda^{*},0)\right\rVert^{2}_{\dagger}}{\left\lVert(Dx^{k-1},z^{k-1},\lambda^{k-1},\Delta
z^{k-1})-(Dx^{*},z^{*},\lambda^{*},0)\right\rVert^{2}_{\dagger}}\leq 1,$

where we have defined the norm:

$\left\lVert(Dx,z,\lambda,\Delta
z)\right\rVert_{\dagger}\coloneqq\sqrt{\left\lVert(Dx,z,\lambda)\right\rVert_{*}^{2}+\tau_{z}\left\lVert\Delta
z\right\rVert^{2}}.$

The above is equivalent to the definition of Q-sublinear convergence [5,
Proposition 1.1d] which proves the claim. ∎

## References

  * [1] L. Grüne and J. Pannek, _Nonlinear Model Predictive Control_. Springer, 2017.
  * [2] J. Gondzio and A. Grothey, “Exploiting structure in parallel implementation of interior point methods for optimization,” _Computational Management Science_ , vol. 6, no. 2, pp. 135–160, 2009.
  * [3] N. Chiang, C. G. Petra, and V. M. Zavala, “Structured nonconvex optimization of large-scale energy systems using PIPS-NLP,” in _2014 Power Systems Computation Conference_. IEEE, 2014, pp. 1–7.
  * [4] D. Kourounis, A. Fuchs, and O. Schenk, “Toward the next generation of multiperiod optimal power flow solvers,” _IEEE Transactions on Power Systems_ , vol. 33, no. 4, pp. 4005–4014, 2018.
  * [5] D. P. Bertsekas, _Constrained optimization and Lagrange multiplier methods_. Academic press, 2014.
  * [6] B. Houska, J. Frasch, and M. Diehl, “An augmented lagrangian based algorithm for distributed nonconvex optimization,” _SIAM Journal on Optimization_ , vol. 26, no. 2, pp. 1101–1127, 2016.
  * [7] J.-H. Hours and C. N. Jones, “A parametric nonconvex decomposition algorithm for real-time and distributed NMPC,” _IEEE Transactions on Automatic Control_ , vol. 61, no. 2, pp. 287–302, 2016.
  * [8] S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein, “Distributed optimization and statistical learning via the alternating direction method of multipliers,” _Foundations and Trends in Machine Learning_ , vol. 3, no. 1, pp. 1–122, 2011.
  * [9] G. Li and T. K. Pong, “Global convergence of splitting methods for nonconvex composite optimization,” _SIAM Journal on Optimization_ , vol. 25, no. 4, pp. 2434–2460, 2015.
  * [10] T. Lin, S. Ma, and S. Zhang, “On the global linear convergence of the ADMM with multiblock variables,” _SIAM Journal on Optimization_ , vol. 25, no. 3, pp. 1478–1497, 2015.
  * [11] S. Magnússon, P. C. Weeraddana, M. G. Rabbat, and C. Fischione, “On the Convergence of Alternating Direction Lagrangian Methods for Nonconvex Structured Optimization Problems,” _IEEE Transactions on Control of Network Systems_ , vol. 3, no. 3, pp. 296–309, 2016.
  * [12] M. Hong, Z.-Q. Luo, and M. Razaviyayn, “Convergence analysis of alternating direction method of multipliers for a family of nonconvex problems,” _SIAM Journal on Optimization_ , vol. 26, no. 1, pp. 337–364, 2016.
  * [13] M. Hong and Z.-Q. Luo, “On the linear convergence of the alternating direction method of multipliers,” _Mathematical Programming_ , vol. 162, no. 1-2, pp. 165–199, 2017.
  * [14] J. G. Melo and R. D. Monteiro, “Iteration-complexity of a Jacobi-type non-Euclidean ADMM for multi-block linearly constrained nonconvex programs,” _arXiv preprint arXiv:1705.07229_ , 2017.
  * [15] J. Bolte, S. Sabach, and M. Teboulle, “Proximal alternating linearized minimization for nonconvex and nonsmooth problems,” _Mathematical Programming_ , vol. 146, no. 1, pp. 459–494, 2014.
  * [16] Q. Liu, X. Shen, and Y. Gu, “Linearized ADMM for nonconvex nonsmooth optimization with convergence analysis,” _IEEE Access_ , vol. 7, pp. 76 131–76 144, 2019.
  * [17] Y. Wang, W. Yin, and J. Zeng, “Global convergence of ADMM in nonconvex nonsmooth optimization,” _Journal of Scientific Computing_ , vol. 78, no. 1, pp. 29–63, 2019.
  * [18] B. Jiang, T. Lin, S. Ma, and S. Zhang, “Structured nonconvex and nonsmooth optimization: algorithms and iteration complexity analysis,” _Computational Optimization and Applications_ , vol. 72, no. 1, pp. 115–157, 2019.
  * [19] K. Sun and X. A. Sun, “A two-level distributed algorithm for nonconvex constrained optimization,” _arXiv preprint arXiv:1902.07654_ , 2019.
  * [20] W. Tang and P. Daoutidis, “Fast and stable nonconvex constrained distributed optimization: the ELLADA algorithm,” _Optimization and Engineering_ , pp. 1–43, 2021.
  * [21] Y. Yang, G. Hu, and C. J. Spanos, “A proximal linearization-based decentralized method for nonconvex problems with nonlinear constraints,” _arXiv preprint arXiv:2001.00767_ , 2020.
  * [22] Q. Shi and M. Hong, “Penalty Dual Decomposition Method for Nonsmooth Nonconvex Optimization—Part I: Algorithms and Convergence Analysis,” _IEEE Transactions on Signal Processing_ , vol. 68, pp. 4108–4122, 2020.
  * [23] D. Zhu, L. Zhao, and S. Zhang, “A first-order primal-dual method for nonconvex constrained optimization based on the augmented lagrangian,” _arXiv preprint arXiv:2007.12219_ , 2020.
  * [24] S. M. Harwood, “Analysis of the alternating direction method of multipliers for nonconvex problems,” in _SN Operations Research Forum_ , vol. 2, no. 1. Springer, 2021, pp. 1–29.
  * [25] G. Scutari, F. Facchinei, and L. Lampariello, “Parallel and distributed methods for constrained nonconvex optimization—part i: Theory,” _IEEE Transactions on Signal Processing_ , vol. 65, no. 8, pp. 1929–1944, 2016.
  * [26] X. Yi, S. Zhang, T. Yang, T. Chai, and K. H. Johansson, “Linear convergence of first-and zeroth-order primal-dual algorithms for distributed nonconvex optimization,” _arXiv preprint arXiv:1912.12110_ , 2019.
  * [27] W. Deng, M.-J. Lai, Z. Peng, and W. Yin, “Parallel multi-block ADMM with $o(1/k)$ convergence,” _Journal of Scientific Computing_ , vol. 71, no. 2, pp. 712–736, 2017.
  * [28] G. Chen and M. Teboulle, “A proximal-based decomposition method for convex minimization problems,” _Mathematical Programming_ , vol. 64, no. 1, pp. 81–101, 1994.
  * [29] G. Banjac, K. Margellos, and P. J. Goulart, “On the Convergence of a Regularized Jacobi Algorithm for Convex Optimization,” _IEEE Transactions on Automatic Control_ , vol. 63, no. 4, pp. 1113–1119, 2018.
  * [30] N. Chatzipanagiotis and M. M. Zavlanos, “On the convergence of a distributed Augmented Lagrangian method for nonconvex optimization,” _IEEE Transactions on Automatic Control_ , vol. 62, no. 9, pp. 4405–4420, 2017.
  * [31] R. T. Rockafellar and R. J.-B. Wets, _Variational Analysis_. Springer, Berlin, Heidelberg, 1998.
  * [32] J. Nocedal and S. Wright, _Numerical optimization_. Springer New York, 2006.
  * [33] S. Frank and S. Rebennack, “An introduction to optimal power flow: Theory, formulation, and examples,” _IIE transactions_ , vol. 48, no. 12, pp. 1172–1197, 2016.
  * [34] A. Maffei, D. Meola, G. Marafioti, G. Palmieri, L. Iannelli, G. Mathisen, E. Bjerkan, and L. Glielmo, “Optimal power flow model with energy storage, an extension towards large integration of renewable energy sources.” _IFAC Proceedings Volumes_ , vol. 47, no. 3, pp. 9456–9461, 2014, 19th IFAC World Congress.
  * [35] A. Castillo, C. Laird, C. A. Silva-Monroy, J.-P. Watson, and R. P. O’Neill, “The unit commitment problem with ac optimal power flow constraints,” _IEEE Transactions on Power Systems_ , vol. 31, no. 6, pp. 4853–4866, 2016\.
  * [36] D. K. Molzahn, I. A. Hiskens _et al._ , “A survey of relaxations and approximations of the power flow equations,” _Foundations and Trends® in Electric Energy Systems_ , vol. 4, no. 1-2, pp. 1–221, 2019.
  * [37] R. D. Zimmerman, C. E. Murillo-Sánchez, and R. J. Thomas, “Matpower: Steady-state operations, planning, and analysis tools for power systems research and education,” _IEEE Transactions on power systems_ , vol. 26, no. 1, pp. 12–19, 2010.
  * [38] ISO New England, “Hourly Real-Time System Demand,” https://www.iso-ne.com/isoexpress/web/reports/load-and-demand/-/tree/dmnd-rt-hourly-sys, Online.
  * [39] I. Dunning, J. Huchette, and M. Lubin, “JuMP: A modeling language for mathematical optimization,” _SIAM Review_ , vol. 59, no. 2, pp. 295–320, 2017.
  * [40] A. Wächter and L. T. Biegler, “On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming,” _Mathematical programming_ , vol. 106, no. 1, pp. 25–57, 2006.
  * [41] S. S. Vazhkudai _et al._ , “The Design, Deployment, and Evaluation of the CORAL Pre-Exascale Systems,” in _SC18: International Conference for High Performance Computing, Networking, Storage and Analysis_ , 2018, pp. 661–672.
  * [42] A. V. Fiacco, “Sensitivity analysis for nonlinear programming using penalty methods,” _Mathematical programming_ , vol. 10, no. 1, pp. 287–311, 1976\.
  * [43] D. Bertsekas, _Nonlinear programming_. Athena Scientific, 1999.

Anirudh Subramanyam is a postdoctoral researcher in the Mathematics and
Computer Science Division at Argonne National Laboratory. He obtained his
bachelor’s degree from the Indian Institute of Technology, Bombay and his
Ph.D. from Carnegie Mellon University, both in chemical engineering. His
research interests are in computational methods for nonlinear and discrete
optimization under uncertainty with applications in energy, transportation and
process systems.
---
Youngdae Kim is a postdoctoral researcher in the Mathematics and Computer
Science Division at Argonne National Laboratory. He received the B.Sc. and
M.Sc. degrees in computer science and engineering from Pohang University of
Science and Technology, Pohang, South Korea, and the Ph.D. degree in computer
sciences from the University of Wisconsin-Madison. His research interests
include distributed optimization methods using hardware accelerators with
applications in power systems.
---
Michel Schanen is an assistant computer scientist in the Mathematics and
Computer Science Division at Argonne National Laboratory. He received his
M.Sc. in 2008 from Rheinisch-Westfälische Technische Hochschule (RWTH) Aachen,
Germany and in 2014 his Ph.D. in computer science from RWTH Aachen, Germany.
His interests lie in automatic differentiation and applications in large-scale
computing.
---
François Pacaud is a postdoctoral researcher in the Mathematics and Computer
Science Division at Argonne National Laboratory. He obtained his M.Sc. in 2015
from Mines ParisTech, Paris, France, and in 2018 his Ph.D. in applied
mathematics from the École des Ponts ParisTech, Paris, France. His research
interests encompasses stochastic and nonlinear optimization, with application
in energy systems.
---
Mihai Anitescu is a senior computational mathematician in the Mathematics and
Computer Science Division at Argonne National Laboratory and a professor in
the Department of Statistics at the University of Chicago. He obtained his
engineer diploma (electrical engineering) from the Polytechnic University of
Bucharest in 1992 and his Ph.D. in applied mathematical and computational
sciences from the University of Iowa in 1997. He specializes in the areas of
numerical optimization, computational science, numerical analysis and
uncertainty quantification in which he has published more than 100 papers in
scholarly journals and book chapters. He has been recognized for his work in
applied mathematics by his selection as a SIAM Fellow in 2019.
---

The submitted manuscript has been created by UChicago Argonne, LLC, Operator
of Argonne National Laboratory (“Argonne”). Argonne, a U.S. Department of
Energy Office of Science laboratory, is operated under Contract No. DE-
AC02-06CH11357. The U.S. Government retains for itself, and others acting on
its behalf, a paid-up nonexclusive, irrevocable worldwide license in said
article to reproduce, prepare derivative works, distribute copies to the
public, and perform publicly and display publicly, by or on behalf of the
Government. The Department of Energy will provide public access to these
results of federally sponsored research in accordance with the DOE Public
Access Plan (http://energy.gov/downloads/doe-public-access-plan).