# Distributed Dual Gradient Tracking for Priority-Considered Load Shedding

Ismi Rosyiana Fitri, Jung-Su Kim Ismi Rosyiana Fitri and Jung-Su Kim are with
the Department of Electrical and Information Engineering, Seoul National
University of Science and Technology, Seoul 01811, Korea e-mail:
<EMAIL_ADDRESS>This work has been submitted to the IEEE for possible
publication. Copyright may be transferred without notice, after which this
version may no longer be accessible.

###### Abstract

This paper studies two fundamental problems in power systems: the economic
dispatch problem (EDP) and load shedding. For the EDP, an extension of the
problem considering the transmission losses is presented. Because the
optimization problem is non-convex owing to quadratic equality constraints
modeling the transmission losses, a convex relaxation method is presented.
Under a particular assumption, it is shown that a distributed algorithm can be
designed to solve the considered EDP. Furthermore, this work aims to handle an
overloading problem by employing an optimal load shedding method. Emphasis is
placed on scheduling the load shedding when there exist some priorities on the
loads. First, a novel optimization problem is proposed to obtain the load
shedding according to the predefined priority order. Then, a distributed
algorithm solving the optimization problem is presented. Finally, this paper
presents how to integrate the proposed algorithms, i.e., the EDP and the load
shedding, in a distributed manner. Simulation results are presented to
demonstrate the effectiveness of the proposed method.

##  I Introduction

###  I-A Background and motivation

There has been a growing interest in developing distributed energy resource
(DER) systems to reduce gas emissions and mitigate global warming. Over the
last few years, the economic dispatch problem (EDP) has been studied to
utilize and manage DER efficiently. The problem is formulated to allocate
power from the generators such that the aggregate generation cost is minimized
while achieving supply-demand balance [1].

Traditionally, the EDP is solved in a centralized manner using global system-
wide information. However, such an approach suffers from single-point failure;
thus, it is inherently not suitable for distributed energy systems [2]. To
deal with this problem, distributed control approaches are widely studied to
make up for the weak point of the centralized controller [3].

Owing to the uncertainty and varying characteristics of DERs, the available
power from the DERs and energy storage may not be sufficient to meet the load
requirements. Thus, load shedding may be necessary to avoid a total collapse
of the electrical supply [4]. The improved communication technologies and
sensing capabilities of a distributed energy system give the operator a chance
to handle this problem through demand-side management. For example, under a
demand response (DR) program, the utilities may request that energy users
curtail a portion or all of their load power demand in exchange for some
incentives, discounts, or other concessions. In light of this, it is
appropriate to consider varying load shedding scenarios due to various
customers’ choices on the offered load shedding program in real-time
operation.

This paper aims to develop a distributed energy management that is applicable
for ensuring the distributed energy systems’ reliability and efficiency. We
are particularly interested in solving two fundamental problems in power
systems: the EDP and load shedding. This work attempts to allocate the optimal
dispatched power through the EDP when transmission loss is considered. In
addition, it aims to design an algorithm for scheduling the load shedding when
there are loads that agree to be shed first due to some prior contract. A
distributed energy management is then formulated so that the proposed
algorithm for the load shedding can be helpful for ensuring the sufficiency of
the available power to supply the remaining demand, meaning that the
feasibility of the EDP is guaranteed.

###  I-B Literature review and problem descriptions

####  I-B1 Distributed EDP

Recent developments in distributed algorithms for the EDP have heightened the
need for distributed agents to agree on the incremental costs [5]. In [6, 7,
8], consensus-based algorithms are developed particularly for a quadratic
function representing the generation cost. In order to consider general convex
cost functions, the algorithms are designed using the dual problem of the EDP,
e.g., in [9, 10, 11]. Furthermore, it is shown in [12] that the convergence of
the existing algorithms can be improved when the generation costs are also
smooth functions.

The works mentioned above do not consider the effects of transmission losses
that exist in practical implementation. The studies in [13] and [14] have
given an account of transmission losses in the EDP. The approach is able to
yield the supply-demand balance by calculating the power loss at every
iteration of the EDP’s algorithm and then adding it to the constant demand. In
[15], the loss model in a quadratic function is fully accounted for obtaining
the optimal dispatched power. Because of the complicated quadratic function
modeling the power losses, the study considers an algorithm which includes an
inner loop of consensus steps per iteration. It is well known that an inner
loop places a more considerable communication burden and needs extra
coordination among agents to terminate the inner loop. In [16, 17], the loss
model is simplified. Hence, the EDP can be solved in a distributed manner, and
the algorithm does not require any inner loop steps. Given the existing works
above, this paper is interested in studying an extension of the traditional
EDP that considers the quadratic function loss model but still allows us to
solve the problem without any need for a consensus inner loop.

####  I-B2 Priority-considered load shedding

To date, several studies have examined a suitable scheme for scheduling the
load shedding, e.g., [18, 19]. Heuristically, the priority-based load shedding
can be obtained by gradually allocating the load shedding from the highest to
the lowest-priority customers [20, 21]. In [22, 23], the scheduled load
shedding is determined by minimizing some convex functions representing the
cost or the adversity of the load shedding. Larger weights are assigned to
high-priority customers so that the resulting load shedding is in accordance
with the pre-selected priority order. In response to the time-varying energy
market and new forecast data, the desired priority list or the buses’ maximum
allowed load shed might change. As a result, the pre-selected weights may no
longer work; thus, they must be re-tuned, which is non-trivial in real-time
operation. This work aims to design a distributed algorithm for a priority-
considered load shedding that minimizes the need for tuning parameters in
order to handle the changes in the desired load shedding scenario.

###  I-C Problem description and statement of contribution

The main contributions of this paper are as follows:

  * •

In this paper, the EDP is formulated such that transmission losses are taken
into account. The proposed EDP is a non-convex optimization problem because of
a quadratic equality constraint denoting the supply-demand balance requirement
under consideration of transmission losses. This paper shows that a convex
relaxation method can be used to obtain the solution of the formulated non-
convex optimization.

  * •

Under some assumptions, it is shown that the EDP can be solved in a
distributed manner such that the method does not require any inner loop
algorithms.

  * •

This paper proposes a novel optimization for the priority-considered load
shedding such that it does not require a re-tuning procedure to deal with
possibly time-varying priority orders. A distributed algorithm to solve the
proposed optimization problem is presented.

Moreover, several supporting consensus algorithms are presented for the
purpose of integrating the two proposed algorithms, i.e., the EDP and the load
shedding, in a distributed manner. Note that when an overloading condition
occurs, an EDP is not feasible due to the lack of available power. Thus, the
load shedding can be applied to restore the feasibility of the EDP.

The paper is constructed as follows. Section II describes the problem
definitions of the EDP and the priority-considered load shedding. Section III
and IV present the algorithm development for the EDP and the priority-
considered load shedding, respectively. Section V discussed the necessary
additional algorithm for integrating both proposed methods. The proposed
algorithms are validated through numerical simulations, which are given in
Section VI.

####  I-C1 Notation

We use boldface to distinguish vector $\bm{x}$ in $\mathbb{R}^{n}$ from the
scalar $x$ in $\mathbb{R}$. Given a vector $\bm{x}$, we denote $\|\bm{x}\|$ as
its Euclidean norm, that is $\|\bm{x}\|=\sqrt{\bm{x}^{T}\bm{x}}$. Given matrix
$A\in\mathbb{R}^{n\times n}$, $A\geq 0\leavevmode\nobreak\ (A>0)$ is used to
denote a positive semi-definite (definite) matrix $A$. The cardinality of a
set $\mathcal{M}$ is denoted by $|\mathcal{M}|$. Consider a set $\mathcal{X}$
defined as $\mathcal{X}=\\{x\in\mathbb{R}|x^{min}\leq x\leq x^{max}\\}$. We
use $\mbox{relint}(\mathcal{X})$ to denote the relative interior of set
$\mathcal{C}$ given by
$\mbox{relint}(\mathcal{X})=\\{x\in\mathbb{R}|x^{min}<x<x^{max}\\}$ for any
$x^{min}$ and $x^{max}$ satisfying $x^{min}<x^{max}$. In this paper, an
overbar is used to denote the optimal solution of an optimization problem. For
example, $\accentset{\rule{2.79996pt}{0.8pt}}{x}$ represents the optimal
solution of $\min_{x}F(x)$. In order to avoid a heavy notation, the
superscript is used in this paper. For instance, $\eta^{(1)},\eta^{(2)}$, and
$\eta^{(3)}$ represent three different variables. For scalars
$x_{1},x_{2},\ldots,x_{n}$, $\bm{x}=[x_{i}]$, $i\in\\{1,\ldots,n\\}$, is the
aggregate vector, i.e., $\bm{x}=[x_{1},x_{2},\ldots,x_{n}]^{T}$. Hence,
$\bm{x}=[x^{(i)}]$, $i\in\\{1,\ldots,n\\}$, is equivalent to
$\bm{x}=[x^{(1)},x^{(2)},\ldots,x^{(n)}]^{T}$.

##  II Problem Definition and Assumptions

Consider the power network with $n$ buses consisting of a generator or (and) a
load where each bus is incorporated with a controller, that is an agent, such
that each agent is able to communicate with at least one other agent. The
communication topology is modeled as an undirected graph
$\mathcal{G}=(\mathcal{V},\mathcal{E})$, where
$\mathcal{V}:=\\{1,\ldots,n\\}$. Denoting $\mathcal{N}_{i}$ as the set of
neighbors of the $i$th agent, $j\in\mathcal{N}_{i}$ implies
$i\in\mathcal{N}_{j}$. This paper aims to develop a distributed algorithm for
scheduling the generation or (and) the load shedding for each bus. For this
purpose, the agents are allowed to exchange their local information with their
neighboring agents through communication network $\mathcal{G}$ satisfying the
following assumption.

######  Assumption 1.

$\mathcal{G}$ is an undirected and connected graph [24].

###  II-A EDP considering transmission losses

Consider a set of buses $\mathcal{V}:=\\{1,\ldots,n\\}$, where $x_{i}$ and
$d_{i}$, $i\in\mathcal{V}$, are the $i$th bus’s dispatched power and forecast
load demand, respectively. In this paper, the EDP is formulated such that
power losses are taken into account.

####  II-A1 Transmission loss model

Denoting $\bm{x}=[x_{i}]$, $i\in\mathcal{V}$, in this paper, the power losses
are modeled as follows

$\displaystyle\Upsilon(\bm{x})=\sum_{i\in\mathcal{V}}\sum_{j\in\mathcal{V}}x_{i}b_{ij}x_{j}+\sum_{i\in\mathcal{V}}b_{0i}x_{i}+b_{00},$
(1)

where $B=[b_{ij}]\in\mathbb{R}^{n\times n}$ and $B=B^{T}$. Without loss of
generality and for simplicity, $b_{0i}=b_{00}=0,i=\\{1,\ldots,n\\}$, is used
hereafter.

######  Assumption 2.

$B\geq 0$.

Note that under this assumption, there exists a real positive semi-definite
matrix $R=[r_{ij}]\in\mathbb{R}^{n\times n}$ such that $R^{T}R=B$.

######  Assumption 3.

The $i$th agent knows the $i$th column of $R$.

######  Remark 1.

The $b$ terms (i.e., $b_{ij},b_{0i}$, and $b_{00}$) in (1) are called loss-
coefficient which mostly can only be computed by the central controller as
extensive amount of global information is required. For example, in [25, Ch.
13.3], the power losses are computed according to the parameters of
transmission lines. For the purpose of solving an EDP problem considering the
transmission loss in a distributed manner, the operator inevitably has to
broadcast the necessary information regarding the loss model to the agents.
For instance, the studies in [15, 14] require the operator to broadcast the
value of $b_{i1},b_{i2},\ldots,b_{in}$ to the $i$th agent since $b_{ij}$ can
not be computed locally by any agent (including the $i$th agent). In
Assumption 3, we assume the $i$th agent receives the information of the $i$th
column of $R$. In the light of distributed systems, such assumptions make the
algorithm depend on the information from the operator. Nevertheless, since the
loss coefficient do not change frequently, it is sufficient for the grid
operator to send the information to the agents from time to time [25, Ch.
13.3].

####  II-A2 EDP formulation

This work aims to solve the following problem

$\displaystyle\min_{\begin{subarray}{c}\bm{x}=[x_{i}],x_{i}\in\mathcal{X}_{i}\\\
i=1,\ldots,n\end{subarray}}$
$\displaystyle\sum_{i\in\mathcal{V}}C_{i}(x_{i}),$ (2a) subject to
$\displaystyle\Upsilon(\bm{x})+\sum\limits_{i\in\mathcal{V}}d_{i}-\sum\limits_{i\in\mathcal{V}}x_{i}=0,$
(2b)

where $\Upsilon(\mathbf{x})$ is given in (1), $C_{i}(\cdot)$ denotes the
generation cost of the $i$th generator, and

$\displaystyle\mathcal{X}_{i}:=\\{x_{i}|x_{i}^{min}\leq x_{i}\leq
x_{i}^{max}\\}.$

For the buses with no generator, $x_{i}^{min}=x_{i}^{max}=0$ and $C_{i}=0$,
otherwise $0\leq x_{i}^{min}<x_{i}^{max}$. Moreover, if the $i$th bus has no
load, then $d_{i}=0$. Note that (2b) represents the supply-demand balance
constraint considering the transmission losses.

####  II-A3 Assumptions for the EDP

The following describes the assumptions considered in this paper.

######  Assumption 4.

For all $i\in\mathcal{V}$, the function
$C_{i}:\mathbb{R}\rightarrow\mathbb{R}$ is proper, closed, and
$c_{i}$-strongly convex, i.e., for all $x,y\in\mathbb{R}$ we have

$\displaystyle C_{i}(y)\geq C_{i}(x)+\nabla
C_{i}(x)\cdot(y\\!-\\!x)+\frac{c_{i}}{2}(y-x)^{2}.$

Moreover, $\nabla C_{i}(x_{i})\geq 0$ for all $x_{i}\in\mathcal{X}_{i}$ and
$\nabla C_{i}(x_{i})>0$ for all $x_{i}\in\mbox{relint}(\mathcal{X}_{i})$.

######  Assumption 5.

Denoting $\bm{x}^{min}\\!=\\![x_{i}^{min}]$,
$\bm{x}^{max}\\!=\\![x_{i}^{max}]$, for all $i\in\mathcal{V}$, it follows that

  1. A5.1

$\sum_{i\in\mathcal{V}}x_{i}^{min}<\sum_{i\in\mathcal{V}}d_{i}+\Upsilon(\bm{x}^{min})$,

  2. A5.2

$\Upsilon(\bm{x}^{min})\\!\leq\\!\sum_{i\in\mathcal{V}}(x_{i}^{min})$,
$\Upsilon(\bm{x}^{max})\\!\leq\\!\sum_{i\in\mathcal{V}}(x_{i}^{max})$,

where $0\leq x_{i}^{min}\leq x_{i}^{max}$ and $d_{i}\geq 0$.

Note that it is easy to satisfy Assumption 4. For instance, it is widely known
that the generation cost can be modeled as
$C(x)=\mathsf{a}x^{2}+\mathsf{b}x+\mathsf{c}$ with $\mathsf{a}_{i}\geq
0,\mathsf{b}_{i}\geq 0$, $\mathsf{c}\in\mathbb{R}$, and $x>0$. Moreover,
Assumption A5.1 indicates that the required supply-demand balance (2b) is
feasible under the minimum-capacity constraint, and A5.2 means that the total
transmission losses are smaller than the total available power. Under the
assumption that (2) is feasible, we aim to solve the problem using a fully
distributed algorithm. The next section presents the problem formulation for
the load shedding. Note that the load shedding can be applied for handling the
case when an overloading condition occurs.

###  II-B Priority-considered load shedding

Let $y_{tot}$ be the total amount of load shedding to be scheduled in the load
buses and $y_{i}$ be the shed load at the $i$th bus such that
$\sum_{i=1}^{n}y_{i}=y_{tot}$. Denoting $y_{i}^{max}$ as the maximum allowed
load shedding at the $i$th bus, it follows that $0\leq y_{i}\leq y_{i}^{max}$.
In this paper, it is assumed that the load buses are divided into the
following categories.

  * •

Prioritized shedding buses ($i\in\mathcal{M}$): The $i$th bus in this category
is shed according to its predefined priority order $p_{i}$, where $p_{i}$
denotes the priority order for the $i$th bus and the set $\mathcal{M}$ is the
set of buses belonging to this category.

  * •

Regular load buses ($i\\!\in\\!\mathcal{V}\backslash\mathcal{M}$): Upon
shedding the buses belonging to this set, the customers suffer from not using
their appliances and $D_{i}(y_{i})$ represents the value of damages or loss
due to the shed load $y_{i}$. Moreover, to alleviate the pain of load
shedding, each customer receives an incentive payment $\mathsf{r}_{i}y_{i}$,
$\mathsf{r}_{i}>0$.

####  II-B1 LS (Load Shedding) Problem

In this work, we aim to design a distributed algorithm for scheduling the load
shedding such that the following conditions are met:

  1. 1.

It is desirable that the load shedding is determined in the buses belonging to
the set $\mathcal{M}$ first, before cutting the loads in
$\mathcal{V}\backslash\mathcal{M}$.

  2. 2.

Among the buses belonging to the set $\mathcal{M}$, $y_{i}$ is selected
according to predefined priority $p_{i}$ without knowing the priority order of
other customers.

  3. 3.

Due to the DR, the customer in the set $\mathcal{V}\backslash\mathcal{M}$ aims
to maximize its utility represented by $\mathsf{r}_{i}y_{i}-D_{i}(y_{i})$
[26]. In other words, $y_{i}$ is determined such that
$\sum_{i\in\mathcal{V\backslash M}}D_{i}(y_{i})-\mathsf{r}_{i}y_{i}$ is
minimized.

####  II-B2 Notation for the priority-considered load shedding

Given the desired priority order $p_{i}$ for each $i$th bus,
$i\in\mathcal{M}$, suppose that the lowest order of priority is $m$, that is,
$\max_{i}p_{i}=:m$. Introduce the set $\mathcal{M}_{\ell}$ given by
$\mathcal{M}_{\ell}:=\\{i|p_{i}=\ell\\},\ell=1,2,\ldots,m$. Then,
$\mathcal{M}$ is their union, i.e.,
$\mathcal{M}=\mathcal{M}_{1}\cup\ldots\cup\mathcal{M}_{m}$.

####  II-B3 Assumption for the load shedding

######  Assumption 6.

For all $i\in\mathcal{V}$, the function
$D_{i}:\mathbb{R}\rightarrow\mathbb{R}$ is proper, closed, and strongly
convex.

In [26], a quadratic function is used to model the load shedding cost
$D_{i}(y_{i})$.

##  III Proposed EDP Considering Power Losses

This section is entirely dedicated to present the solution of the considered
EDP described in the previous section. Let us assume that the EDP (2) is
feasible and that the following holds.

######  Assumption 7.

[Slater’s condition] There exists a point
$(\tilde{x}_{1},\ldots,\tilde{x}_{n})$ that lies in the interior of
$\mathcal{X}_{1}\times\mathcal{X}_{2}\times\ldots\times\mathcal{X}_{n}$ such
that
$\Upsilon(\tilde{\bm{x}})+\sum_{i\in\mathcal{V}}d_{i}-\sum_{i\in\mathcal{V}}\tilde{x}_{i}=0$
where $\tilde{\bm{x}}=[\tilde{x}_{i}]$.

###  III-A Development of a distributed algorithm for the EDP

This section describes several equivalent optimization problems that are
helpful for solving (2) in a distributed manner.

####  III-A1 An equivalent optimization problem for the EDP

Let us introduce a variable $u_{i}\in\mathbb{R}$, $i\in\mathcal{V}$. Note that
the following optimization problem has the same solution as (2).

$\displaystyle\min_{\begin{subarray}{c}\bm{x}=[x_{i}],\bm{u}=[u_{i}],\\\
x_{i}\in\mathcal{X}_{i},\forall i\in\mathcal{V}\end{subarray}}$
$\displaystyle\sum_{i\in\mathcal{V}}C_{i}(x_{i}),$ (3a) subject to
$\displaystyle\sum\limits_{i\in\mathcal{V}}u_{i}^{2}+\sum\limits_{i\in\mathcal{V}}d_{i}-\sum\limits_{i\in\mathcal{V}}x_{i}=0,$
(3b) $\displaystyle R\bm{x}=\bm{u},$ (3c)

where $R^{T}R=B$. The optimization problem is defined such that a new slack
variable $u_{i}$ is added and the loss model $\bm{x}^{T}B\bm{x}$ is
substituted by $\bm{u}^{T}\bm{u}$ where $\bm{u}=R\bm{x}$. Since the feasible
set and the objective function of (3) are the same as those of (2), the
optimal solution of (2) can be found by solving (3).

####  III-A2 Relaxation method for the non-convex EDP

In (3), the problem is non-convex due to the quadratic equality constraint
(3b); thus, obtaining the optimal solution can be non-trivial. To deal with
this problem, consider the following optimization problem

$\displaystyle\min_{\begin{subarray}{c}\bm{x}=[x_{i}],\bm{u}=[u_{i}],\\\
x_{i}\in\mathcal{X}_{i},\forall i\in\mathcal{V}\end{subarray}}$
$\displaystyle\sum_{i\in\mathcal{V}}C_{i}(x_{i}),$ (4a) subject to
$\displaystyle\sum\limits_{i\in\mathcal{V}}u_{i}^{2}+\sum\limits_{i\in\mathcal{V}}d_{i}-\sum\limits_{i\in\mathcal{V}}x_{i}\leq
0,$ (4b) $\displaystyle R\bm{x}=\bm{u}.$ (4c)

Note that this problem is convex since the quadratic equality in (3b) is
relaxed into an inequality. Nevertheless, the following theorem shows that the
optimal solution of the relaxed problem (4) is also that of the EDP (3),
equivalently (2).

######  Theorem 1.

[Relaxation] Suppose that Assumptions 2, 4, 5, and 7 hold. If
$(\accentset{\rule{2.79996pt}{0.8pt}}{x}_{i},\accentset{\rule{2.79996pt}{0.8pt}}{u}_{i}),i\in\mathcal{V},$
is the optimal solution of (4), then,
$(\accentset{\rule{2.79996pt}{0.8pt}}{x}_{i},\accentset{\rule{2.79996pt}{0.8pt}}{u}_{i}),i\in\mathcal{V},$
is also the optimal solution of (3).

###### Proof.

Let $r_{ij}$ be the element of $R$ on the $i$th row and the $j$th column, and
the Lagrangian function of (4) given as

$\displaystyle\mathcal{L}(\bm{x},\bm{u},\lambda,\bm{\xi})\\!\\!\\!$
$\displaystyle=$
$\displaystyle\\!\\!\\!\sum_{i\in\mathcal{V}}C_{i}(x_{i})+\lambda\sum\limits_{i\in\mathcal{V}}\left(u_{i}^{2}+d_{i}-x_{i}\right)$
(5)
$\displaystyle+\sum_{j\in\mathcal{V}}\xi^{(j)}\left(-u_{j}+\sum\limits_{i\in\mathcal{V}}r_{ji}x_{i}\right),$

where $x_{i}\in\mathcal{X}_{i},\bm{x}=[x_{i}],\bm{u}=[u_{i}]$,
$i\in\mathcal{V},$ $\bm{\xi}=[\xi^{(j)}]$, $j\in\mathcal{V},$ and $\lambda\geq
0$. Denote
$(\accentset{\rule{2.79996pt}{0.8pt}}{x}_{i},\accentset{\rule{2.79996pt}{0.8pt}}{u}_{i})$
and
$(\accentset{\rule{2.79996pt}{0.8pt}}{\lambda},\accentset{\rule{2.79996pt}{0.8pt}}{\bm{\xi}})$
as the primal and dual optimal solution of (4), respectively. Since (4) is
convex and the Slater’s condition holds, the KKT conditions give the necessary
and sufficient condition for the optimal solution of (4). In this paper, the
theorem is proven by showing that the following holds.

$\displaystyle\accentset{\rule{2.79996pt}{0.8pt}}{x}_{i}\in\mathcal{X}_{i},\quad\\!(\accentset{\rule{2.79996pt}{0.8pt}}{x}_{i},\accentset{\rule{2.79996pt}{0.8pt}}{u}_{i})=\operatorname*{argmin}\limits_{x_{i}\in\mathcal{X}_{i},u_{i}}\mathcal{L}(\bm{x},\bm{u},\lambda,\bm{\xi}),i\in\mathcal{V},$
(6a)
$\displaystyle\accentset{\rule{2.79996pt}{0.8pt}}{\lambda}>0,\sum_{i\in\mathcal{V}}\accentset{\rule{2.79996pt}{0.8pt}}{u}_{i}^{2}+d_{i}-\accentset{\rule{2.79996pt}{0.8pt}}{x}_{i}=0,$
(6b)
$\displaystyle\sum_{i\in\mathcal{V}}r_{ji}\accentset{\rule{2.79996pt}{0.8pt}}{x}_{i}=\accentset{\rule{2.79996pt}{0.8pt}}{u}_{j},\leavevmode\nobreak\
j=1,\ldots,n.$ (6c)

Observe that the equality in (6b) shows that the optimal solution of (4) is
also the optimal solution of (3). The optimality conditions given in (6) is
proven in the Appendix. The proof shows that the monotonicity of the objective
functions $C_{i}(x_{i})$ on $\mbox{relint}(X_{i})$ plays a key role in
obtaining the result. ∎

In light of Theorem 1, the optimal solution (3) can be obtained through the
relaxed problem (4).

####  III-A3 Utilizing a redundant constraint for the distributed EDP

Consider the following optimization problem

$\displaystyle\min_{\begin{subarray}{c}\bm{x}=[x_{i}],\bm{u}=[u_{i}],\\\
x_{i}\in\mathcal{X}_{i},\forall i\in\mathcal{V}\end{subarray}}$
$\displaystyle\sum_{i\in\mathcal{V}}C_{i}(x_{i}),$ (7) subject to
$\displaystyle\eqref{eqbalrel},\leavevmode\nobreak\
\eqref{pxu},u_{i}\in\mathcal{U},\leavevmode\nobreak\ i\in\mathcal{V},$

where $\mathcal{U}:=\\{u|-u^{max}\leq u\leq u^{max}\\}$ and

$\displaystyle u^{max}\\!=\\!\sum_{i\in\mathcal{V}}(r_{i}^{max}x_{i}^{max}),$
(8)

and $r_{i}^{max}:=\max_{j=1,\ldots,n}|r_{ji}|$. Observe that (7) is formulated
by adding a new constraint $u_{i}\in\mathcal{U}$ to (4), and that the
following holds

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\sum_{j\in\mathcal{V}}r_{ij}x_{j}^{max}\leq
u^{max},\leavevmode\nobreak\ \forall i,j\in\mathcal{V}.$

Moreover, since $x_{i}^{min}\leq x_{i}^{max}$,
$-u^{max}\leq\sum_{j\in\mathcal{V}}r_{ij}x_{j}^{min}$, the constraint
$u_{i}\in\mathcal{U}$ is always satisfied as long as $x_{i}$ and $u_{i}$
satisfy the constraints in (4): $x_{i}\in\mathcal{X}_{i}$ and (4c). In other
words, the new constraint $u_{i}\in\mathcal{U}$ does not alter the feasible
region (7). As a result, the optimal solution of (7) is equivalent to that of
(4). In view this, the optimal solution of (4) can be obtained by solving (7).
Clearly, it is redundant to specify $u_{i}\in\mathcal{U}$ given that
$x_{i}\in\mathcal{X}_{i}$ and (4c) hold. However, it is shown later (c.f.
Theorem 2) that the constraint $u_{i}\in\mathcal{U}$ is useful for designing a
distributed algorithm to obtain the optimal solution of (7).

####  III-A4 Dual problem

Let $\accentset{\rule{2.79996pt}{0.8pt}}{\mathcal{S}}$ be the set of the
optimal solutions of (7). In view of $\mathcal{X}_{i}$, the feasible region of
(7) is closed and bounded. From Assumption 4, the objective function is convex
and continuous. Thus, according to the Weierstrass theorem [27], there exists
an optimal point
$(\accentset{\rule{2.79996pt}{0.8pt}}{\bm{x}},\accentset{\rule{2.79996pt}{0.8pt}}{\bm{u}})\in\accentset{\rule{2.79996pt}{0.8pt}}{\mathcal{S}}$
that achieves a global optimum. Moreover, under the Slater’s condition given
in Assumption 7, the strong duality holds and the set of the dual optimal
solutions is nonempty. Hence, an algorithm based on a primal-dual method can
be used for obtaining the optimal solution of (7). To this end, consider
$\lambda_{i},\xi_{i}^{(1)},\ldots,\xi_{i}^{(n)}$ as the $i$th agent’s
estimates of the Lagrangian multipliers $\lambda,\xi^{(1)},\ldots,\xi^{(n)}$,
and the following local Lagrangian function

$\displaystyle\mathcal{L}_{i}(x_{i},u_{i},\lambda_{i},\bm{\xi}_{i})\\!\\!\\!$
$\displaystyle=$
$\displaystyle\\!\\!\\!C_{i}(x_{i})+\lambda_{i}(u_{i}^{2}+d_{i}-x_{i})+\xi_{i}^{(i)}(r_{ii}x_{i}-u_{i})+\\!\\!\\!\\!\sum_{j\in\mathcal{V},j\not=i}\xi_{i}^{(j)}r_{ji}x_{i},$

where $\bm{\xi}_{i}=[\xi_{i}^{(1)},\ldots,\xi_{i}^{(n)}]^{T}$,
$x_{i}\in\mathcal{X}_{i}$, and $u_{i}\in\mathcal{U}$. Observe that
$\sum_{i\in\mathcal{V}}\mathcal{L}_{i}=\mathcal{L}$ if and only if
$\bm{\xi}_{i}=\bm{\xi}_{j}$ and $\lambda_{i}=\lambda_{j}$ hold for all
$i,j\in\mathcal{V}$. Denoting
$h_{i}(\lambda_{i},\bm{\xi}_{i})=\inf_{x_{i}\in\mathcal{X}_{i},u_{i}\in\mathcal{U}}\mathcal{L}_{i}(x_{i},u_{i},\lambda_{i},\bm{\xi}_{i})$,
the dual problem of (7) can be formulated as follows:

$\displaystyle\min\limits_{\begin{subarray}{c}\lambda_{i}\geq
0,\bm{\xi}_{i},\\\ i=1,\ldots,n\end{subarray}}$
$\displaystyle\sum_{i\in\mathcal{V}}q_{i}(\lambda_{i},\bm{\xi}_{i})$ (9a)
subject to $\displaystyle\lambda_{i}=\lambda_{j},\leavevmode\nobreak\
\leavevmode\nobreak\ \forall i,j\in\mathcal{V},$ (9c)
$\displaystyle\bm{\xi}_{i}=\bm{\xi}_{j},\leavevmode\nobreak\
\leavevmode\nobreak\ \forall i,j\in\mathcal{V},$

where $q_{i}(\lambda_{i},\bm{\xi}_{i}):=-h_{i}(\lambda_{i},\bm{\xi}_{i})$.
Note that the objective function of (9) is a summation of several convex
functions. Utilizing this structure, the focus of the problem is placed on
designing a distributed algorithm for (9).

###  III-B Proposed distributed algorithm for the EDP

In this section, the algorithm for solving the dual problem (9) is designed
based on the distributed subgradient method [28]. In addition to
$x_{i},u_{i},\lambda_{i},$ and $\bm{\xi}_{i}$, let us introduce two variables
$v_{i}\in\mathbb{R}$ and $\bm{w}_{i}\in\mathbb{R}^{n}$, $i\in\mathcal{V}$. The
sequences of $x_{i}(k),u_{i}(k)$ for $k\geq 0$ are defined as follows

$\displaystyle\\!\\!v_{i}(k)\\!\\!\\!\\!\\!$ $\displaystyle=$
$\displaystyle\\!\\!\\!\\!\\!\sum\limits_{j\in\mathcal{V}}a_{ij}\lambda_{j}(k),$
(10a) $\displaystyle\\!\\!\bm{w}_{i}(k)\\!\\!\\!\\!\\!$ $\displaystyle=$
$\displaystyle\\!\\!\\!\\!\\!\sum\limits_{j\in\mathcal{V}}a_{ij}\bm{\xi}_{j}(k),$
(10b) [resume]
$\displaystyle[x_{i}(k\\!+\\!1),u_{i}(k\\!+\\!1)]\\!=\\!\inf_{\begin{subarray}{c}x_{i}\in\mathcal{X}_{i},u_{i}\in\mathcal{U}\end{subarray}}\\!\\!\\!\mathcal{L}_{i}(x_{i},u_{i},v_{i}(k),\bm{w}_{i}(k)),$
(11a)
$\displaystyle\lambda_{i}(k\\!+\\!1)\\!=\\!\mathcal{P}_{\mathcal{R}}\\!\left[v_{i}(k)\\!+\\!\alpha(k)\left(u_{i}^{2}(k\\!+\\!1)\\!+\\!d_{i}\\!-\\!x_{i}(k\\!+\\!1)\right)\right]\\!\\!,\leavevmode\nobreak\
\leavevmode\nobreak\ $ (11b)
$\displaystyle\bm{\xi}_{i}(k\\!+\\!1)\\!=\\!\bm{w}_{i}(k)+\alpha(k)\bm{g}_{i}(x_{i}(k+1),u_{i}(k+1)),$
(11c)

where $\alpha(k)$ is a stepsize sequence, $\mathcal{P}_{\mathcal{R}}$ is a
projection onto the set $\mathcal{R}=\\{\lambda|\lambda\geq 0\\}$,
$\bm{g}_{i}(x_{i},u_{i}):=[g_{i}^{(1)}(x_{i},u_{i}),\ldots,g_{i}^{(n)}(x_{i},u_{i})]^{T}$,
and

$\displaystyle
g_{i}^{(j)}(x_{i},u_{i})=\left\\{\\!\\!\begin{array}[]{cc}r_{ii}x_{i}-u_{i}&\mbox{if
}j=i,\\\ r_{ji}x_{i}&\mbox{otherwise}.\end{array}\right.$ (14)

Moreover, $A:=[a_{ij}]$ is a doubly stochastic matrix, i.e.,
$\sum_{i\in\mathcal{V}}a_{ij}=\sum_{j\in\mathcal{V}}a_{ji}=1$ with $a_{ii}>0$,
such that $a_{ij}>0$ if and only if $j\in\mathcal{N}_{i}$ otherwise $a_{ij}=0$

In view of (10a)-(10b), the $i$th agent broadcasts $\lambda_{i}(k)$ and
$\bm{\xi}_{i}(k)$ to its neighboring agents which are used to update
$v_{j}(k)$ and $\bm{w}_{j}(k)$, $j\in\mathcal{N}_{i}$. In step (11b)-(11c),
$\lambda_{i}$ and $\bm{\xi}_{i}$ are updated by moving the current values
along the subgradient of $q_{i}$ at $v_{i}(k)$ and $\bm{w}_{i}(k)$ with step
size $\alpha(k)$. Projection $\mathcal{P}_{\mathcal{R}}[\cdot]$ is used to
make $\lambda_{i}$ nonnegative. With Assumption 3 in mind, the proposed
algorithm (10) can be implemented in a distributed way if $u_{max}$ is known
to all agents. Note that, from the definition of $u_{max}$ given in (8), an
average consensus method can be applied to generate $u_{max}$ in a distributed
manner before running (10). Furthermore, the proposed method does not need any
inner loop algorithm at each iteration $k$.

###  III-C Convergence analysis

Let
$(\accentset{\rule{2.79996pt}{0.8pt}}{x}_{i},\accentset{\rule{2.79996pt}{0.8pt}}{u}_{i})$
and
$(\accentset{\rule{2.79996pt}{0.8pt}}{\lambda},\accentset{\rule{2.79996pt}{0.8pt}}{\bm{\xi}})$
be the optimal solution of (7) and (9) respectively, $i\in\mathcal{V}$. The
following shows the convergence of (10).

######  Theorem 2.

Suppose that the EDP (2) is feasible, and that Assumptions 1, 2, 4, 5, and 7
hold true. Consider the sequences (10) with a step size $\alpha(k)$ satisfying
$\sum_{k=0}^{\infty}\alpha(k)=\infty$ and
$\sum_{k=0}^{\infty}\alpha^{2}(k)<\infty$. Then, for all $i\in\mathcal{V}$,
$(\lambda_{i}(k),\bm{\xi}_{i}(k))\\!\\!\rightarrow\\!\\!(\accentset{\rule{2.79996pt}{0.8pt}}{\lambda},\accentset{\rule{2.79996pt}{0.8pt}}{\bm{\xi}})$
and
$(x_{i}(k),u_{i}(k))\rightarrow(\accentset{\rule{2.79996pt}{0.8pt}}{x}_{i},\accentset{\rule{2.79996pt}{0.8pt}}{u}_{i})$
as $k\rightarrow\infty$.

###### Proof.

The feasibility of (2) implies the feasibility of (7). Note that the dual
problem (9) has the similar form to the optimization problem considered in
[28]. According to [28, Proposition 4], the convergence of
$(\lambda_{i}(k),\bm{\xi}_{i}(k))$ to dual optimal point
$(\accentset{\rule{2.79996pt}{0.8pt}}{\lambda},\accentset{\rule{2.79996pt}{0.8pt}}{\bm{\xi}})$
can be shown when the subgradient of $q_{i}$ is uniformly bounded. Observe
that $|x_{i}(k)|\leq x_{i}^{max}$ and $|u_{i}^{(j)}(k)|\leq u^{max}$ for all
$k\geq 1$. Having this, the subgradient of $q_{i}$ is uniformly bounded, i.e.,
there exists scalar $\Gamma>0$ such that
$|u_{i}(k)^{2}+d_{i}-x_{i}(k)|\leq\Gamma$ and
$|g_{i}^{(j)}(x_{i}(k),u_{i}(k))|\leq\Gamma$ for all $k\geq 1$,
$i,j\in\\{1,\ldots,n\\}$. As a result,
$(\lambda_{i}(k),\bm{\xi}_{i}(k))\\!\\!\rightarrow\\!\\!(\accentset{\rule{2.79996pt}{0.8pt}}{\lambda},\accentset{\rule{2.79996pt}{0.8pt}}{\bm{\xi}})$
as $k\rightarrow\infty$.

Since (6) holds for the primal-dual optimal solution of (4), it also holds for
the primal-dual optimal solution of (7). Given that $C_{i}$ is strongly convex
and $\accentset{\rule{2.79996pt}{0.8pt}}{\lambda}>0$,
$\mathcal{L}_{i}(x_{i},u_{i},\lambda_{i},\bm{\xi}_{i})$ is a strongly convex
function for fixed $\lambda_{i}=\accentset{\rule{2.79996pt}{0.8pt}}{\lambda}$
and ${\bm{\xi}_{i}}=\accentset{\rule{2.79996pt}{0.8pt}}{\bm{\xi}}$. In other
words, there exists $\beta_{i}>0$ satisfying

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\frac{\beta_{i}}{2}\left\|\left[\begin{array}[]{c}\accentset{\rule{2.79996pt}{0.8pt}}{x}_{i}\\\
\accentset{\rule{2.79996pt}{0.8pt}}{u}_{i}\end{array}\right]-\left[\begin{array}[]{c}{x}_{i}(k+1)\\\
u_{i}(k+1)\end{array}\right]\right\|^{2}$ (19)
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\leq\mathcal{L}_{i}(\accentset{\rule{2.79996pt}{0.8pt}}{x}_{i},\accentset{\rule{2.79996pt}{0.8pt}}{u}_{i},\accentset{\rule{2.79996pt}{0.8pt}}{\lambda},\accentset{\rule{2.79996pt}{0.8pt}}{\bm{\xi}})\\!-\mathcal{L}_{i}(x_{i}(k\\!+\\!1),u_{i}(k\\!+\\!1),\accentset{\rule{2.79996pt}{0.8pt}}{\lambda},\accentset{\rule{2.79996pt}{0.8pt}}{\bm{\xi}})$
$\displaystyle\\!\\!\\!\\!\\!\\!-\nabla\mathcal{L}_{i}(x_{i}(k\\!+\\!1),u_{i}(k\\!+\\!1),\accentset{\rule{2.79996pt}{0.8pt}}{\lambda},\accentset{\rule{2.79996pt}{0.8pt}}{\bm{\xi}})^{T}\left(\left[\begin{array}[]{c}\\!\\!\\!\accentset{\rule{2.79996pt}{0.8pt}}{x}_{i}\\!\\!\\!\\\
\\!\\!\\!\accentset{\rule{2.79996pt}{0.8pt}}{u}_{i}\\!\\!\\!\end{array}\right]\\!-\\!\left[\begin{array}[]{c}\\!\\!\\!{x}_{i}(k\\!+\\!1)\\!\\!\\!\\\
\\!\\!\\!u_{i}(k\\!+\\!1)\\!\\!\\!\end{array}\right]\right)\\!.\leavevmode\nobreak\
$ (24)

Given that $x_{i}(k+1)$ and $u_{i}(k+1)$ satisfy (11a), observe that
$\lim_{k\rightarrow\infty}(v_{i}(k),\bm{w}_{i}(k))=(\accentset{\rule{2.79996pt}{0.8pt}}{\lambda},\accentset{\rule{2.79996pt}{0.8pt}}{\bm{\xi}})$
implies
$\nabla\mathcal{L}(x_{i}(k\\!+\\!1),u_{i}(k\\!+\\!1),\accentset{\rule{2.79996pt}{0.8pt}}{\lambda},\accentset{\rule{2.79996pt}{0.8pt}}{\bm{\xi}})\rightarrow\bm{0}$
as $k\rightarrow\infty$. With this and (19) in mind, we have

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\lim\limits_{k\rightarrow\infty}\sum_{i\in\mathcal{V}}\frac{\beta_{i}}{2}\left\|\left[\begin{array}[]{c}\accentset{\rule{2.79996pt}{0.8pt}}{x}_{i}\\\
\accentset{\rule{2.79996pt}{0.8pt}}{u}_{i}\end{array}\right]-\left[\begin{array}[]{c}{x}_{i}(k+1)\\\
u_{i}(k+1)\end{array}\right]\right\|^{2}$ (29)
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\leq\lim\limits_{k\rightarrow\infty}\sum_{i\in\mathcal{V}}\left(\mathcal{L}_{i}(\accentset{\rule{2.79996pt}{0.8pt}}{x}_{i},\accentset{\rule{2.79996pt}{0.8pt}}{u}_{i},\accentset{\rule{2.79996pt}{0.8pt}}{\lambda},\accentset{\rule{2.79996pt}{0.8pt}}{\bm{\xi}})\\!-\mathcal{L}_{i}(x_{i}(k\\!+\\!1),u_{i}(k\\!+\\!1),\accentset{\rule{2.79996pt}{0.8pt}}{\lambda},\accentset{\rule{2.79996pt}{0.8pt}}{\bm{\xi}})\right).$

Claim 1: The strong-duality of (7) yields

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\lim\limits_{k\rightarrow\infty}\\!\sum_{i\in\mathcal{V}}\\!\\!\mathcal{L}_{i}(\accentset{\rule{2.79996pt}{0.8pt}}{x}_{i},\accentset{\rule{2.79996pt}{0.8pt}}{u}_{i},\accentset{\rule{2.79996pt}{0.8pt}}{\lambda},\accentset{\rule{2.79996pt}{0.8pt}}{\bm{\xi}})\\!-\\!\mathcal{L}_{i}(x_{i}(k\\!+\\!1),u_{i}(k\\!+\\!1),\accentset{\rule{2.79996pt}{0.8pt}}{\lambda},\accentset{\rule{2.79996pt}{0.8pt}}{\bm{\xi}})\\!=\\!0.$
(30)

The proof of the claim is given in Appendix. In view of (29) and (30),
$(x_{i}(k),u_{i}(k))\rightarrow(\accentset{\rule{2.79996pt}{0.8pt}}{x}_{i},\accentset{\rule{2.79996pt}{0.8pt}}{u}_{i})$
as $k\rightarrow\infty$. ∎

##  IV Priority-considered Load Shedding

In this section, at first, a convex optimization problem is proposed for
addressing the desired load shedding result described in LS Problem. Then, a
distributed algorithm for solving the proposed optimization problem is
presented.

###  IV-A Proposed optimization problem

Let us assume that the total amount of load shedding $y_{tot}$ is known at
least by the central operator, and that each $i$th agent has a knowledge on
its average value $s_{i}$, i.e.,

$\displaystyle s_{i}:=y_{tot}/n.$

Consider a new variable $z_{i},i\in\mathcal{M},$ denoting a nonnegative slack
variable which is used to maximize $y_{i}$ belonging to the set $\mathcal{M}$.
Moreover, define

$\displaystyle\leavevmode\nobreak\ \mathcal{Z}_{i}:=\\{z_{i}|0\leq z_{i}\leq
ns_{i}\\},\quad\mbox{and}\quad\mathcal{Y}_{i}:=\\{y_{i}|0\leq y_{i}\leq
y_{i}^{max}\\}.$

For all $i\in\mathcal{M}_{\ell},\ell=1,\ldots,m,$ consider the function
$f_{i}(z_{i},y_{i})$ defined by

$\displaystyle f_{i}(z_{i},y_{i}):=\kappa
z_{i}^{2}+\left(y_{i}-\frac{ns_{i}}{p_{i}}\right)^{2}.$

Recall that $p_{i}$ is the desired priority order of the $i$th bus, e.g.,
$p_{i}=1\Leftrightarrow i\in\mathcal{M}_{1}$. Having this, the proposed
optimization problem for the priority-considered load shedding is given as
follows:

$\displaystyle\min_{\begin{subarray}{c}y_{i},\forall i\in\mathcal{V}\\\
z_{i},\forall i\in\mathcal{M}\end{subarray}}\\!\\!\\!\\!$
$\displaystyle\\!\\!\\!\\!\\!\sum\limits_{i\in\mathcal{V}\backslash\mathcal{M}}\\!\\!\\!\\!\\!\left(D_{i}(y_{i})-\mathsf{r}_{i}y_{i}\right)+\\!\\!\sum\limits_{i\in\mathcal{M}}\\!\\!f_{i}(z_{i},y_{i})\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ $ subject to
$\displaystyle\\!\\!\\!\\!\\!\\!\leavevmode\nobreak\
y_{i}\in\mathcal{Y}_{i},\leavevmode\nobreak\ \forall
i\in\mathcal{V},\leavevmode\nobreak\
z_{i}\in\mathcal{Z}_{i},\leavevmode\nobreak\ \forall i\in\mathcal{M}$ (31d)
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\mbox{if }m=0:$
$\displaystyle\\!\\!\sum_{i\in\mathcal{V}}s_{i}-\sum_{i\in\mathcal{V}}y_{i}=0,$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\mbox{if }m\not=0:$
$\displaystyle\\!\\!\sum\limits_{i\in\mathcal{V}}s_{i}-\\!\\!\\!\sum\limits_{i\in\mathcal{M}_{1}}\\!\\!\left(z_{i}\\!+\\!y_{i}\right)=0,$
$\displaystyle\\!\\!\\!\\!\\!\\!\sum\limits_{i\in\mathcal{M}_{\ell-1}}\\!\\!\\!\\!\\!z_{i}\\!-\\!\\!\\!\sum\limits_{i\in\mathcal{M}_{\ell}}\\!\\!\\!\left(z_{i}\\!+\\!y_{i}\right)\\!=\\!0,\leavevmode\nobreak\
\\!\forall\ell\\!=\\!2,\ldots,m,$
$\displaystyle\\!\\!\\!\\!\sum\limits_{i\in\mathcal{M}_{m}}\\!\\!\\!\\!z_{i}\\!-\\!\\!\\!\sum_{i\in\mathcal{V}\backslash\mathcal{M}}y_{i}=0,$

where $\kappa\geq 1$ is a tuning parameter. Note that the proposed
optimization problem is convex and consists of $m+1$ equality constraints.

####  IV-A1 Priority-considered load shedding through (31)

The objective function consists of two terms. The first term
$\sum_{i\in\mathcal{V}\backslash\mathcal{M}}D_{i}(y_{i})-\mathsf{r}_{i}y_{i}$
aims to maximize the utility of the customers belonging to the set
$\mathcal{V}\backslash\mathcal{M}$. Moreover, the following shows that
$f_{i}(z_{i},y_{i})$ is used to make the optimal solution of (31) adhere to
the desired priority list.

  1. (i)

Substituting of (31d) and (31d) into (31d) yields (31d), i.e.,
$\sum_{i\in\mathcal{M}}{y_{i}}=y_{tot}$.

  2. (ii)

With (31d) in mind, minimizing $\kappa\sum_{i\in\mathcal{M}_{1}}\\!z_{i}^{2}$
is equivalent to minimizing
$\sum_{i\in\mathcal{V}}s_{i}-\sum_{i\in\mathcal{M}_{1}}\\!y_{i}$. Hence,
giving $\kappa>\nabla D_{i}(y_{i})$ for all
$i\in\mathcal{V}\backslash\mathcal{M}$ assigns more cost to
$\kappa\sum_{i\in\mathcal{M}_{1}}\\!z_{i}^{2}=\kappa\left(\sum_{i\in\mathcal{V}}s_{i}-\sum_{i\in\mathcal{M}_{1}}\\!y_{i}\right)^{2}$
compared to $\sum_{i\in\mathcal{V}\backslash\mathcal{M}}D_{i}(y_{i})$. This
means sufficiently large $\kappa$ results in more load shedding on the buses
belonging to the set $\mathcal{M}_{1}$ compared to other buses.

  3. (iii)

Equality (31d) for $\ell=2$ can be written as
$\sum_{i\in\mathcal{M}_{2}}z_{i}=\sum_{i\in\mathcal{M}_{1}}z_{i}-\sum_{i\in\mathcal{M}_{2}}\\!y_{i}$.
Substituting $\sum_{i\in\mathcal{M}_{1}}z_{i}$ with its definition given in
(31d) yields
$\sum_{i\in\mathcal{M}_{2}}z_{i}=\sum_{i\in\mathcal{V}}s_{i}-\sum_{i\in\mathcal{M}_{1}\cup\mathcal{M}_{2}}\\!y_{i}$.
Hence, minimizing $\kappa\sum_{i\in\mathcal{M}_{2}}\\!z_{i}^{2}$ with
$\kappa>\nabla D_{i}(y_{i})$ for all $i\in\mathcal{V}\backslash\mathcal{M}$
can be seen as an attempt to allocate more load shedding on the buses
belonging to the set $\mathcal{M}_{1}\cup\mathcal{M}_{2}$.

  4. (iv)

Note that from the points (ii) and (iii), $y_{i},i\in\mathcal{M}_{1},$ is
prioritized through the minimization of
$\kappa\sum_{i\in\mathcal{M}_{1}}\\!z_{i}^{2}$ and
$\kappa\sum_{i\in\mathcal{M}_{2}}\\!z_{i}^{2}$. Meanwhile,
$y_{i},i\in\mathcal{M}_{2},$ is prioritized through the minimization of
$\kappa\sum_{i\in\mathcal{M}_{2}}\\!z_{i}^{2}$ only. In light of this,
$y_{i},i\in\mathcal{M}_{1},$ is prioritized more than
$y_{i},i\in\mathcal{M}_{2}$.

  5. (v)

This reasoning can be extended to show that the constraint (31d) helps to
generate $y_{i}$ according to its desired priority order $p_{i}$.

In addition, the objective function $\left(y_{i}-{ns_{i}}/{p_{i}}\right)^{2}$
is used to make the term strongly convex. Hence, the optimal solution
$\accentset{\rule{2.79996pt}{0.8pt}}{y}_{i},i\in\mathcal{M},$ can be obtained
using a primal-dual approach (c.f. Theorem 3). By definition, the function
minimizes the difference between $y_{i}$ and ${y_{tot}}/{\ell}$,
$i\in\mathcal{M}_{\ell},\ell=1,\ldots,m$.

####  IV-A2 Feasibility of (31)

The following describes the condition when (31) is feasible.

######  Assumption 8.

$\sum_{i\in\mathcal{V}}y_{i}^{max}>\sum_{i\in\mathcal{V}}s_{i}=y_{tot}$.

Note that since the constraint set in (31) is compact and the objective
function is continuous, the set of the optimal solution of (31) is nonempty.
Moreover, under Assumption 8, there exists
$\tilde{y}_{1},\ldots,\tilde{y}_{n}$ that lie in the interior of
$\mathcal{Y}_{1}\times\ldots\times\mathcal{Y}_{n}$ satisfying
$\sum_{i\in\mathcal{V}}s_{i}-\tilde{y}_{i}=0$, meaning that the strong duality
holds for (31).

###  IV-B Distributed priority-considered load shedding

Here, we present an algorithm to solve (31) in a distributed manner

####  IV-B1 Dual Decomposition

Let $\bm{y}=[y_{i}]$, ${i\in\mathcal{V}}$, and $\bm{z}=[z_{i}]$,
${i\in\mathcal{M}}$, consider the Lagrangian function of (31) given as

$\displaystyle\mathcal{L}^{ls}(\bm{y},\bm{z},\eta^{(1)},\ldots,\eta^{(m+1)})=\sum\limits_{i\in\mathcal{V}\backslash\mathcal{M}}\\!\\!\\!\\!\\!\left(\\!D_{i}(y_{i})\\!-\\!\mathsf{r}_{i}y_{i}\\!\right)\\!+\\!\\!\\!\sum\limits_{i\in\mathcal{M}}\\!\\!f_{i}(z_{i},y_{i})+\sum_{i\in\mathcal{V}}\sum_{\nu=1}^{m+1}\eta^{(\nu)}\mathsf{g}_{i}^{(\nu)}(z_{i},y_{i}),$
(32)

where $y_{i}\in\mathcal{Y}_{i}$ for all ${i\in\mathcal{V}}$ and
$z_{i}\in\mathcal{Z}_{i}$ for all ${i\in\mathcal{M}}$. Moreover, with a slight
abuse of notation,

$\displaystyle\mathsf{g}_{i}^{(m+1)}(z_{i},y_{i})$ $\displaystyle=$
$\displaystyle\left\\{\begin{array}[]{cl}s_{i}-y_{i}&m=0,\\\
z_{i}&m\not=0,i\in\mathcal{M}_{m},\\\
-y_{i}&m\not=0,i\in\mathcal{V}\backslash\mathcal{M},\\\
0&\mbox{otherwise},\end{array}\right.$ (33e)
$\displaystyle\mathsf{g}_{i}^{(1)}(z_{i},y_{i})$ $\displaystyle=$
$\displaystyle\left\\{\begin{array}[]{cl}s_{i}-y_{i}-z_{i}&i\in\mathcal{M}_{1}\\\
s_{i}&\mbox{otherwise}.\end{array}\right.$ (33h)

For $\nu=2,\ldots,m$, $\mathsf{g}_{i}^{(\nu)}(y_{i},z_{i})$ is defined as
follows

[resume] $\displaystyle\mathsf{g}_{i}^{(\nu)}(z_{i},y_{i})$ $\displaystyle=$
$\displaystyle\left\\{\begin{array}[]{cl}-y_{i}-z_{i}&i\in\mathcal{M}_{\ell},\nu=\ell\\\
z_{i}&i\in\mathcal{M}_{\ell},\nu=\ell+1\\\
0&\mbox{otherwise},\end{array}\right.$ (34d)

where $\ell\in\\{1,\ldots,m\\}$. Denoting $\eta_{i}^{(\nu)}$ as the $i$th
agent’s estimate of $\eta^{(\nu)}$, and $\bm{\eta}_{i}=[\eta_{i}^{(\nu)}]$,
$\nu=1,\ldots,m+1$, consider the following local Lagrangian function

$\displaystyle\mathcal{L}^{ls}_{i}(y_{i},z_{i},\bm{\eta}_{i}):=F_{i}(y_{i},z_{i})+\sum_{\nu=1}^{m+1}\eta_{i}^{(\nu)}\mathsf{g}_{i}^{(\nu)}(z_{i},y_{i})$
(35)

where, with a slight abuse of notation,
$F_{i}(z_{i},y_{i})=D_{i}(y_{i})-\mathsf{r}_{i}y_{i}$ for
$i\in\mathcal{V}\backslash\mathcal{M}$, and
$F_{i}(z_{i},y_{i})=f_{i}(z_{i},y_{i})$ for $i\in\mathcal{M}$. Note that
$\sum_{i\in\mathcal{V}}\mathcal{L}_{i}^{ls}=\mathcal{L}^{ls}$ if and only if
$\bm{\eta}_{i}=\bm{\eta}_{j}$ holds for all $i,j\in\mathcal{V}$. Thus, the
dual problem of (31) can be formulated as follows

$\displaystyle\max\limits_{\begin{subarray}{c}\bm{\eta}_{i}\\\
i=1,\ldots,n\end{subarray}}$
$\displaystyle\sum_{i\in\mathcal{V}}h_{i}^{ls}(\bm{\eta}_{i})$ (36) subject to
$\displaystyle\bm{\eta}_{i}=\bm{\eta}_{j},\leavevmode\nobreak\
\leavevmode\nobreak\ \forall i,j\in\mathcal{V}$

where
$h_{i}^{ls}(\bm{\eta}_{i}):=\inf_{y_{i}\in\mathcal{Y}_{i},z_{i}\in\mathcal{Z}_{i}}\mathcal{L}_{i}^{ls}(y_{i},z_{i},\bm{\eta}_{i})$.
As in the previous section, a method based on the distributed sub-gradient
approach is used to solve (36) in a distributed manner.

####  IV-B2 Distributed algorithm

Let us introduce $\phi_{i}^{(\nu)},\nu=1,\ldots,m+1$, $i\in\mathcal{V}$.
Denoting $\bm{\eta}_{i}=[\eta_{i}^{(\nu)}]$ and
$\bm{\phi}_{i}=[\phi_{i}^{(\nu)}],\nu=1,\ldots,m+1$, consider the following
sequences for solving (31)

$\displaystyle[z_{i}(k+1),y_{i}(k+1)]=\inf_{y_{i}\in\mathcal{Y}_{i},z_{i}\in\mathcal{Z}_{i}}\mathcal{L}_{i}^{ls}(y_{i},z_{i},\bm{\phi}_{i}(k)),$
(37a)

where

[resume] $\displaystyle\bm{\phi}_{i}(k+1)\\!\\!$ $\displaystyle=$
$\displaystyle\\!\\!\sum_{j\in\mathcal{V}}a_{ij}\bm{\eta}_{j}(k+1),$ (38a)
$\displaystyle\bm{\eta}_{i}(k+1)\\!\\!$ $\displaystyle=$
$\displaystyle\\!\\!\bm{\phi}_{i}(k)+\alpha(k)\bm{\mathsf{g}}_{i}(z_{i}(k\\!+\\!1),y_{i}(k\\!+\\!1)),$
(38b)

and
$\bm{\mathsf{g}}_{i}(z_{i},y_{i})=[\mathsf{g}_{i}^{(1)}(z_{i},y_{i}),\ldots,\mathsf{g}_{i}^{(m+1)}(z_{i},y_{i})]^{T}$.
The proposed algorithm (37) can be implemented in a distributed manner under
the assumption that $m$ is known to all agents.

####  IV-B3 Convergence Analysis

The following shows the convergence of (37).

######  Theorem 3.

Suppose that Assumptions 1, 6 and 8 hold, and that $\alpha(k)$ satisfies
$\sum_{k=0}^{\infty}\alpha(k)=\infty$ and
$\sum_{k=0}^{\infty}\alpha^{2}(k)<\infty$. Then, the sequences (37) solves
optimization problem (31), i.e.,
$\bm{\eta}_{i}(k)\rightarrow\accentset{\rule{2.79996pt}{0.8pt}}{\bm{\eta}}$
and
$(z_{i}(k),y_{i}(k))\rightarrow(\accentset{\rule{2.79996pt}{0.8pt}}{z}_{i},\accentset{\rule{2.79996pt}{0.8pt}}{y}_{i})$,
as $k\rightarrow\infty$ for all $i\in\mathcal{V}$.

See Appendix for the proof. When the desired priority list is changed, note
that (37) requires only one tuning parameter $\kappa>\nabla D_{i}(y_{i})$ for
all $y_{i}\in\mathcal{Y}_{i},i\in\mathcal{V}\backslash\mathcal{M}$, under
assumption that the value of $m$ and $s_{i}$ are informed to all agents at
every initial step of (37). Given this, a question arises on developing a
method such that $m$ and $s_{i}$ can be obtained by exchanging information
between the neighbor agents; thus, the communication between the central
controller and the agents can be reduced. The algorithm is discussed in the
next section.

##  V DISTRIBUTED ENERGY MANAGEMENT : EDP + LOAD SHEDDING

This section presents how to integrate the proposed algorithms given in the
previous sections. To this end, algorithms used to generate the value of
$y_{tot}$ and $m$ in a distributed manner are proposed.

###  V-A Solving an overloading condition in the EDP

In Section III, it is assumed that the EDP is feasible. However, there can be
a case when there is no $\tilde{x}_{i}\in\mathcal{X}_{i}$ such that
$\Upsilon(\tilde{\bm{x}})+\sum_{i\in\mathcal{V}}d_{i}-\sum_{i\in\mathcal{V}}\tilde{x}_{i}=0$
holds, meaning that an overloading condition occurs. To handle this, this
section aims to discuss how to obtain $\sum_{i\in\mathcal{V}}s_{i}=y_{tot}$ in
a distributed manner so that the feasibility of the EDP can be recovered.

####  V-A1 A feasibility problem

To find the necessary load shedding $y_{tot}=\sum_{i\in\mathcal{V}}s_{i}$,
consider the following feasibility problem

$\displaystyle\begin{array}[]{cc}\mbox{{P1:}}\nobreak\leavevmode\nobreak\leavevmode\nobreak\leavevmode\nobreak\leavevmode\nobreak\leavevmode\nobreak\leavevmode&\begin{array}[]{c}\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\min\limits_{\begin{subarray}{c}s_{i}\geq
0,\bm{x}=[x_{i}],\\\
x_{i}\in\mathcal{X}_{i},i\in\mathcal{V}\end{subarray}}\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \sum_{i\in\mathcal{V}}s_{i}^{2}\\\ \mbox{subject
to\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
}\Upsilon(\bm{x})+\sum\limits_{i\in\mathcal{V}}\left(d_{i}-x_{i}-s_{i}\right)=0.\end{array}\end{array}$

Let us denote $s_{i}^{*}$ and $x_{i}^{*}$ as the optimal solution of this
problem. Note that ${s}_{i}^{*}={s}_{j}^{*}$, $i,j\in\mathcal{V}$. Moreover,
the problem yields the optimal sum of load shedding
$y_{tot}=\sum_{i\in\mathcal{V}}{s}_{i}^{*}$ because $s_{i}^{*}$ is the
smallest value such that
$\sum_{i\in\mathcal{V}}{s}_{i}^{*}=\Upsilon(\bm{x}^{*})+\sum_{i\in\mathcal{V}}\left(d_{i}-{x}_{i}^{*}\right)$
holds, $\bm{x}^{*}=[x_{i}^{*}],i\in\mathcal{V}$.

####  V-A2 Relaxation of feasibility problem P1

Consider the relaxed problem of P1 defined as

$\displaystyle\begin{array}[]{cc}\mbox{{P2:}}\nobreak\leavevmode\nobreak\leavevmode&\begin{array}[]{c}\\!\\!\\!\\!\\!\\!\min\limits_{\begin{subarray}{c}s_{i}\geq
0,\bm{u}=[u_{i}],\bm{x}=[x_{i}],\\\
x_{i}\in\mathcal{X}_{i},i\in\mathcal{V}\end{subarray}}\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\
\sum_{i\in\mathcal{V}}s_{i}^{2}+\tau\sum_{i\in\mathcal{V}}x_{i}^{2},\\\
\mbox{subject to\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
}\\!\\!\\!\sum\limits_{i\in\mathcal{V}}\left(u_{i}^{2}+d_{i}-x_{i}-s_{i}\right)\leq\\!0,\\\
R\bm{x}=\bm{u},\end{array}\end{array}$

where $\tau>0$ is a tuning parameter. Denoting
$\accentset{\rule{2.79996pt}{0.8pt}}{x}_{i},\accentset{\rule{2.79996pt}{0.8pt}}{u}_{i}$
and $\accentset{\rule{2.79996pt}{0.8pt}}{s}_{i}$ as the optimal solution of
this problem, it follows that
$\accentset{\rule{2.79996pt}{0.8pt}}{s}_{i}=\accentset{\rule{2.79996pt}{0.8pt}}{s}_{j}$,
$i,j\in\mathcal{V}$.

Note that P2 is convex as the equality constraint in P1 is relaxed into an
inequality. Moreover, the term $\tau\sum_{i\in\mathcal{V}}x_{i}^{2}$ is
included as an additional objective function with $\tau>0$. Having this, the
objective function of P2 is monotonically increasing on the set
$\\{x_{i},s_{i}|x_{i}\in\mbox{relint}(\mathcal{X}_{i}),s_{i}>0\\}$. Hence,
Theorem 1 can be re-established to show that the optimal solution of this
problem meets equality
$\sum_{i\in\mathcal{V}}\left(\accentset{\rule{2.79996pt}{0.8pt}}{u}_{i}^{2}+d_{i}-\accentset{\rule{2.79996pt}{0.8pt}}{x}_{i}-\accentset{\rule{2.79996pt}{0.8pt}}{s}_{i}\right)=0$,
i.e.,
$\sum_{i\in\mathcal{V}}\accentset{\rule{2.79996pt}{0.8pt}}{s}_{i}=\Upsilon(\accentset{\rule{2.79996pt}{0.8pt}}{\bm{x}}_{i})+\sum_{i\in\mathcal{V}}\left(d_{i}-\accentset{\rule{2.79996pt}{0.8pt}}{x}_{i}\right)$
where
$\accentset{\rule{2.79996pt}{0.8pt}}{\bm{x}}_{i}=[\accentset{\rule{2.79996pt}{0.8pt}}{x}_{i}],i\in\mathcal{V}$.
However, due to the additional term $\tau\sum_{i\in\mathcal{V}}x_{i}^{2}$, we
have
$\sum_{i\in\mathcal{V}}\accentset{\rule{2.79996pt}{0.8pt}}{s}_{i}\geq\sum_{i\in\mathcal{V}}s_{i}^{*}$,
meaning that the optimal solution of P2 is not the optimal solution of P1.
Nevertheless, small value of
$\|s_{i}^{*}-\accentset{\rule{2.79996pt}{0.8pt}}{s}_{i}\|$ can be obtained by
setting $\tau>0$ sufficiently small. With some minor modifications, the
distributed algorithm given in (10) can be adopted for obtaining
$\accentset{\rule{2.79996pt}{0.8pt}}{s}_{i}$ using only local information. The
details are omitted due to page limit.

###  V-B Supporting algorithm for obtaining $m$

####  V-B1 Review of a consensus algorithm

At first, a consensus algorithm studied in [29] is reviewed. Suppose that each
$i$th agent has two local variables $\theta_{i},\omega_{i}\in\mathbb{R}$ and
applies a consensus algorithm given by

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\theta_{i}(k+1)\\!\\!\\!$
$\displaystyle=$
$\displaystyle\\!\\!\\!\omega_{i}(k)+\frac{1}{2}\sum_{j\in\mathcal{N}_{i}}\frac{\omega_{j}(k)-\omega_{i}(k)}{\max\\{|\mathcal{N}|_{i},|\mathcal{N}_{j}|\\}},$
(43a) $\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\omega_{i}(k+1)\\!\\!\\!$
$\displaystyle=$
$\displaystyle\\!\\!\\!\theta_{i}(k)\\!+\\!\left(\\!1\\!-\\!\frac{2}{9n+1}\\!\right)\\!\\!(\theta_{i}(k+1)-\theta_{i}(k)).$
(43b)

Then, the following lemma holds.

######  Lemma 1.

[29, Theorem 2.1] Suppose that Assumption 1 holds and that each node in an
undirected graph $\mathcal{G}$ implements (43) with
$\omega_{i}(0)=\theta_{i}(0)$. Then,
$\lim_{k\rightarrow\infty}\theta_{i}(k)=\frac{1}{n}\sum_{i=1}^{n}\omega_{i}(0)=:\theta_{ss}$
and

$\displaystyle\|\bm{\theta}(k)-\theta_{ss}\bm{1}\|_{2}^{2}\leq
2\left(1-\frac{1}{9n}\right)^{k-1}\|\bm{\theta}(0)-\theta_{ss}\bm{1}\|_{2}^{2}$

where $\bm{\theta}=[\theta_{i}],i=1,\ldots,n.$

Lemma 1 implies that (43) converges to the average of the initial values as
$k\rightarrow\infty$. Moreover,
$\|\bm{\theta}(k)-\theta_{ss}\bm{1}\|_{\infty}\leq\epsilon$ after
$O(n\ln(\|\bm{\theta}(0)-\theta_{ss}\bm{1}\|/\epsilon))$ iterations. The
following section describes the approach based on (43) for obtaining the value
of $m$ in a distributed manner.

####  V-B2 Applying consensus algorithm (43) to generate $m$

Let us consider the following value for the initial condition of (43)

$\displaystyle\omega_{i}(0)$ $\displaystyle=$
$\displaystyle\left\\{\begin{array}[]{cl}{np_{i}}/{|M_{\ell}|}&i\in\mathcal{M}_{\ell},\ell\in[1,m],\\\
{n}/{|\mathcal{V}\backslash\mathcal{M}|}&i\in\mathcal{V}\backslash\mathcal{M},\end{array}\right.$
(46)

with $\theta_{i}(0)=\omega_{i}(0)$. Then,
$\lim_{k\rightarrow\infty}\theta_{i}(k)=\theta_{ss}=1+\sum_{\ell=1}^{m}\ell=(2+{m(m+1)})/{2}$,
$i\in\mathcal{V}$. It is easy to see that $\theta_{ss}$ is always unique for
any $m\in[0,n]$. For example, $\theta_{ss}=1$ when $m=0$, $\theta_{ss}=2$ when
$m=1$, and so on. Since $m$ is upper-bounded by $n$, the consensus value
$\theta_{ss}$ can be decoded to obtain the value of $m$.

In practice, it may be hard to obtain $\theta_{i}=\theta_{ss}$ since it holds
asymptotically. Thus, we may need to use $\theta_{i}(k)$ instead of
$\theta_{ss}$. From Lemma 1, $\theta_{i}(k)=\theta_{ss}\pm\epsilon$ for all
$k\geq K_{n}(\epsilon)$ where
$K_{n}(\epsilon):=O(n\ln(\|\bm{\theta}(0)-\theta_{ss}\bm{1}\|/\epsilon))$. The
following describes the value of
$\theta_{i}(k)\\!\\!=\\!\theta_{ss}\\!\pm\\!\epsilon$ for given different
values of $m$.

Table I: $\theta_{i}(k)\\!\\!=\\!\theta_{ss}\\!\pm\\!\epsilon$ for given different value of $m$. $\theta_{i}(k)=\theta_{ss}\pm\epsilon$ |  $1\pm\epsilon$ |  $2\pm\epsilon$ |  $\ldots$ |  $(2+n^{2}+n\pm 2\epsilon)/{2}$
---|---|---|---|---
$m$ |  0 |  1 |  $\ldots$ |  $\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ n\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ $

In view of Table I, choosing $\epsilon$ satisfying $0<\epsilon<0.5$ yields an
unique set $\theta_{i}(k)\in[\theta_{ss}-\epsilon,\theta_{ss}+\epsilon]$ for
any $m$, $m\in[0,n]$. Having this, the value of $m$ can be still attained by
running the algorithm (43) for at least $K_{n}(\epsilon)$ iterations where
$\epsilon\in(0,0.5)$.

######  Remark 2.

The proposed algorithm (43) can be implemented under the assumption that the
$i$th agent has knowledge of the cardinality of the set $\mathcal{M}_{\ell}$
(if $i\in\mathcal{M}_{\ell}$) or $\mathcal{V}\backslash\mathcal{M}$ (if
$i\in\mathcal{V}\backslash\mathcal{M}$), $\ell\in\\{1,\ldots,m\\}$. In
practice, this information can be obtained from the central controller. Note
that when the desired priority list is changed, then, the central controller
broadcast the new value of $p_{i}$ and $|\mathcal{M}_{\ell}|$ only to the
agents that are affected by the modification. In light of this, the method
requires less communication between the agents and the central controller,
compared to the case when $m$ is broadcast to all agents. Moreover, from the
initial condition (46), the information regarding the choice of the load
shedding program of the $i$th customer is not disclosed to other agents,
meaning that the privacy is preserved.

###  V-C Proposed distributed energy management

The integration of the proposed EDP and the load shedding method is formally
written in Algorithm 1. Bear in mind that the required demand for the EDP is
updated in Algorithm 1 by considering the load shedding, i.e.,
$d_{i}\leftarrow d_{i}-s_{i}$ in Stage 2.

Algorithm 1 Distributed EDP + load shedding

Input from the operator: $p_{i}$, $|\mathcal{M}_{\ell}|$ if
$i\in\mathcal{M}_{\ell}$, $|\mathcal{V}\backslash\mathcal{M}|$ otherwise
Stage 1a : Finding $m$

Initialize $\theta_{i}(0)=\omega_{i}(0)$ using (46)

Run (43) for $k\geq K_{n}(\epsilon)$

Output : All agent know $m$ by using Table I

Stage 1b : Estimating $y_{tot}$

Solve P2 in a distributed manner

Output : Each agent holds $s_{i}$

Stage 2: Determining $x_{i},y_{i}$

Input : $m,s_{i}$

Set $d_{i}\leftarrow d_{i}-s_{i}$

Initialize $\lambda_{i}(0)\\!\geq\\!0,\bm{\xi}_{i}(0)\in\mathbb{R}^{n}$,
$\bm{\eta}_{i}(0)\in\mathbb{R}^{m+1}$

repeat$\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
k=0,1,\ldots$

Exchange $\lambda_{i}(k),\bm{\xi}_{i}(k),\bm{\eta}_{i}(k)$ to neighboring
agents.

Run (10a)-(10b) and (38a)

Calculates $x_{i}(k+1)$ and $y_{i}(k+1)$ using (11a) and (37a)

Update $\lambda_{i}(k+1),\bm{\xi}_{i}(k+1),\bm{\eta}_{i}(k)$ using (11b)-(11c)
and (38b)

until converge

In the next section, numerical simulations are used to demonstrate that the
proposed algorithms yield the desired output of the EDP and the priority-
considered load shedding.

##  VI Numerical Simulations

Note that Stage 2 of Algorithm 1 plays an essential role in determining the
dispatched power and load shedding at each bus. Moreover, Stages 1a and 1b can
be seen as two supporting algorithms. Having this, due to page limit, this
section only presents the simulation results for the algorithms under Stage 2
of Algorithm 1, i.e., the methods given in (10) and (37) are validated. For
this purpose, consider the IEEE-30 bus [30]. The system consists of 6
generators and 24 load buses. Furthermore, it is assumed that the
communication between the buses is modeled as an undirected graph
$\mathcal{G}=(\mathcal{V},\mathcal{E})$, $\mathcal{V}=\\{1,\ldots,30\\}$ where
the communication paths do not necessarily coincide with the power flow.

###  VI-A The proposed EDP considering transmission losses

Table II: Generators parameters ($MU$=money unit). Generator |  $\mathsf{a}_{i}$ $(MU/MW^{2})$ |  $\mathsf{b}_{i}$ $(MU/MW)$ | $x_{i}^{max}(MW)$
---|---|---|---
1 | 0.08 | 2 | 20
2 | 0.06 | 3 | 10
3 | 0.07 | 4 | 30
4 | 0.06 | 4 | 15
5 | 0.08 | 2.5 | 10
6 | 0.08 | 2.5 | 8

The simulation results of the proposed distributed EDP given in (10) is
presented in this section. Table II describes $x_{i}^{max}$ and the parameters
of the cost function of the $i$th generator with
$C_{i}(x_{i})=\mathsf{a}_{i}x_{i}^{2}+\mathsf{b}_{i}x_{i}$. The minimum
generation $x_{i}^{min}$ for all generators are 5MW. The $B$-matrix is
obtained from [30] and is given as follows:

$\displaystyle
B\\!=\\!10^{-2}\\!\\!\left[\begin{array}[]{cccccc}\\!\\!\\!13.82\\!&\\!-2.99\\!&\\!0.44\\!&\\!-0.22\\!&\\!-0.10\\!&\\!-0.08\\!\\!\\!\\\
\\!\\!\\!-2.99\\!&\\!4.87\\!&\\!-0.25\\!&\\!0.04\\!&\\!0.16\\!&\\!0.41\\!\\\
\\!\\!\\!0.44\\!&\\!-0.25\\!&\\!1.82\\!&\\!-0.70\\!&\\!-0.66\\!&\\!-0.66\\!\\\
\\!\\!\\!-0.22\\!&\\!0.04\\!&\\!-0.70\\!&\\!1.37\\!&\\!0.50\\!&\\!0.33\\!\\\
\\!\\!\\!-0.10\\!&\\!0.16\\!&\\!-0.66\\!&\\!0.50\\!&\\!1.09\\!&\\!0.05\\!\\\
\\!\\!\\!-0.08\\!&\\!0.41\\!&\\!-0.66\\!&\\!0.33\\!&\\!0.05\\!&\\!2.44\\!\end{array}\right]$

At first, the proposed relaxed problem (7) is solved in a centralized way to
show that the optimal solution $\accentset{\rule{2.79996pt}{0.8pt}}{x}_{i}$
yields
$\Upsilon(\accentset{\rule{2.79996pt}{0.8pt}}{\bm{x}})+\sum_{i\in\mathcal{V}}d_{i}-\sum_{i\in\mathcal{V}}\accentset{\rule{2.79996pt}{0.8pt}}{x}_{i}=0$,
$\accentset{\rule{2.79996pt}{0.8pt}}{\bm{x}}=[\accentset{\rule{2.79996pt}{0.8pt}}{x}_{i}]$.
The results are depicted in Table III. This implied that the supply-demand
balance is attained using the relaxed-problem (7).

Table III: Optimal solution of (7) obtained using a centralized method.  $\sum_{i\in\mathcal{V}}d_{i}$ | $\accentset{\rule{2.79996pt}{0.8pt}}{\bm{x}}=[\accentset{\rule{2.79996pt}{0.8pt}}{x}_{i}]$ | $\Upsilon(\accentset{\rule{2.79996pt}{0.8pt}}{\bm{x}})-\bm{1}^{T}\accentset{\rule{2.79996pt}{0.8pt}}{\bm{x}}$
---|---|---
36 | [5, 6.0836, 8.8734, 7.31, 8.2366, 6.57] | -36
48 | [5, 7.406, 14.844, 11.544, 10, 8] | -48
55.2 | [5, 9.4079, 19.5281, 15, 10, 8] | -55.2

For the purpose of solving (7) in a distributed setting, the algorithm (10) is
employed with $\alpha_{i}(k)={100}/{k^{0.6}}$ and the load demand is 2MW at
all load buses, meaning that $\sum_{i\in\mathcal{V}}d_{i}=48$ since there are
24 load buses. Figure 1 shows that the total power mismatch converges to zero
as iteration increases. Moreover, Figure 2 describes the dispatched power at
each generator. Note that the optimal solution from the distributed algorithm
converges to the solution by a centralized method shown in Table III.
Moreover, Table IV shows the dispatched power obtained using the proposed
distributed method (10) for given different load demands. Observe that the
results are the same as the solution obtained using a centralized approach
described in Table III.

Figure 1: The dispatched power converges to the value such that $\sum_{i\in\mathcal{V}}d_{i}=\sum_{i\in\mathcal{V}}x_{i}-\Upsilon(\bm{x})$, $\bm{x}=[x_{i}]$. Figure 2: For the given load demand $\sum_{i\in\mathcal{V}}d_{i}=48$, the dispatched power $x_{i}(k)$ converges to $[5,7.38,14.78,11.64,10,8]$.  Table IV: Optimal solution of (7) obtained using the proposed distributed algorithm (10).  $\sum_{i\in\mathcal{V}}d_{i}$ | $\accentset{\rule{2.79996pt}{0.8pt}}{\bm{x}}=[\accentset{\rule{2.79996pt}{0.8pt}}{x}_{i}]$ | $\Upsilon(\accentset{\rule{2.79996pt}{0.8pt}}{\bm{x}})-\bm{1}^{T}\accentset{\rule{2.79996pt}{0.8pt}}{\bm{x}}$
---|---|---
36 | [5, 6.05, 8.82, 7.34, 8.2323, 6.6437] | -35.9999
48 | [5, 7.38, 14.78, 11.64, 10, 8] | -48.0035
55.2 | [5, 9.407, 19.5321, 15, 10, 8] | -55.2027

###  VI-B The proposed priority considered load shedding

Numerical simulations are carried out to investigate the effectiveness of
distributed algorithm (37) to solve LS Problem defined in Section II. The
algorithm is applied with $\alpha(k)=1000/(k+500)$ to allocate the load
shedding for 24 buses, i.e., $i\in\\{7,\ldots,30\\}$. In this study, $m=3$ is
considered where $\mathcal{M}_{1}=\\{7\\}$, $\mathcal{M}_{2}=\\{8,9\\}$, and
$\mathcal{M}_{3}=\\{10\\}$. This means that we would like to shed the load in
the $7$th, $8$th, $9$th, and $10$th bus first before shedding the remaining 20
buses. Moreover, the bus belonging to the set $\mathcal{M}_{1}$ has the
highest, and $\mathcal{M}_{3}$ has the lowest priority to be shed. Consider
$\mathsf{r}_{i}=3$ and $D_{i}(y_{i})=\frac{1}{2}\mathsf{q}_{i}y_{i}^{2}$ where
$\mathsf{q}_{i}=1$ for all $i=\\{11,12,21,22\\}$, $\mathsf{q}_{i}=2$ for all
$i=\\{13,14,23,24\\}$, $\mathsf{q}_{i}=2.5$ for all $i=\\{15,16,25,26\\}$,
$\mathsf{q}_{i}=3$ for all $i=\\{17,18,27,28\\}$, and $\mathsf{q}_{i}=4$ for
all $i=\\{19,20,29,30\\}$.

The algorithm is employed by setting $\kappa=40$ and $y_{i}^{max}=1.2$,
$i\in\mathcal{V}$. At first, let us consider the case when
$\sum_{i\in\mathcal{V}}s_{i}=y_{tot}=1.8$. Figure 3 shows that
$y_{7}=y_{7}^{max}=1.2$ and $y_{8}=y_{9}=0.3$. Since $y_{i}\approx 0$ for all
$i\in\mathcal{V}\backslash\\{7,8,9\\}$, we have
$\sum_{i\in\mathcal{V}}y_{i}=y_{tot}$. Note that $y_{8}=y_{9}=0.3$ because the
load shedding at the bus belonging to the set $\mathcal{M}_{1}$ is not enough
to meet the required load shedding $y_{tot}$. Owing to this, the priority-
considered load shedding scenario is achieved.

Next, simulations are carried out by considering three different values of
$y_{tot}$. The results are shown in Figure 4. When, $y_{tot}=1$, only the bus
belonging to the set $\mathcal{M}_{1}$ is shed. In the case $y_{tot}=4$, we
have $y_{7}=y_{8}=y_{9}=1.2$ and $y_{10}=0.4$. Moreover, when $y_{tot}=6$, we
have $y_{7}=y_{8}=y_{9}=y_{10}=\sum_{i\in
V\backslash\\{7,8,9,10\\}}y_{i}=1.2$. Owing to this, the proposed method is
able to yield the load shedding according to the desired priority list. Figure
5 shows the convergence of the algorithm in the case when $y_{tot}=6$. Note
that $y_{i}$ converges to the value such that its summation over all
$i\in\mathcal{V}$ coincides with the desired value
$\sum_{i\in\mathcal{V}}s_{i}=y_{tot}=6$.

Figure 3: Let $\mathcal{M}_{1}=\\{7\\},\mathcal{M}_{2}=\\{8,9\\}$,
$y_{7}^{max}=1.2$, and $y_{tot}=1.8$, simulation results show that the desired
priority-considered load shedding is obtained.  Figure 4: Simulation results
for different $y_{tot}:1,4,$ and $6$. Given that
$\mathcal{M}_{1}=\\{7\\},\mathcal{M}_{2}=\\{8,9\\}$,
$\mathcal{M}_{3}=\\{10\\}$, and $y_{i}^{max}=1.2$ for all $i\in\mathcal{V},$
the scheduled load shedding follows the desired priority list. Figure 5: For
the case when $\sum_{i\in\mathcal{V}}s_{i}=y_{tot}=6$ the value of
$\sum_{i\in\mathcal{V}}y_{i}-\sum_{i\in\mathcal{V}}s_{i}$ converges to zero as
the iterations increase.

##  VII Conclusion

This paper formulates an optimization problem for the EDP to deal with the
transmission loss power modeled in a quadratic function. In addition, an
optimization problem is proposed to consider some priorities in the load
shedding program. Under some assumptions, distributed algorithms for both
problems are proposed. Moreover, a feasibility problem is adopted to find the
necessary load shedding for handling an overloading condition in a distributed
manner.

## References

  * [1] Ji-Yuan Fan and Lan Zhang, “Real-time economic dispatch with line flow and emission constraints using quadratic programming,” _IEEE Transactions on Power Systems_ , vol. 13, no. 2, pp. 320–325, May 1998.
  * [2] Qiang Wan, W. Zhang, Y. Xu, and I. Khan, “Distributed control for energy management in a microgrid,” in _2016 IEEE/PES Transmission and Distribution Conference and Exposition (T &D)_, Dallas, TX, 2016, pp. 1–5.
  * [3] T. Yang, D. Wu, H. Fang, W. Ren, H. Wang, Y. Hong, and K. H. Johansson, “Distributed energy resource coordination over time-varying directed communication networks,” _IEEE Transactions on Control of Network Systems_ , vol. 6, no. 3, pp. 1124–1134, 2019.
  * [4] F. Lin, “Worst-case load shedding in electric power networks,” _IEEE Transactions on Control of Network Systems_ , vol. 6, no. 3, pp. 1269–1277, 2019\.
  * [5] Z. Zhang and M. Chow, “Convergence analysis of the incremental cost consensus algorithm under different communication network topologies in a smart grid,” _IEEE Transactions on Power Systems_ , vol. 27, no. 4, pp. 1761–1768, 2012.
  * [6] S. Kar and G. Hug, “Distributed robust economic dispatch in power systems: A consensus + innovations approach,” in _2012 IEEE Power and Energy Society General Meeting_ , San Diego, CA, 2012, pp. 1–8.
  * [7] Z. Yang, J. Xiang, and Y. Li, “Distributed consensus based supply–demand balance algorithm for economic dispatch problem in a smart grid with switching graph,” _IEEE Transactions on Industrial Electronics_ , vol. 64, no. 2, pp. 1600–1610, 2017.
  * [8] R. Wang, Q. Li, B. Zhang, and L. Wang, “Distributed consensus based algorithm for economic dispatch in a microgrid,” _IEEE Transactions on Smart Grid_ , vol. 10, no. 4, pp. 3630–3640, 2019.
  * [9] T. T. Doan and C. L. Beck, “Distributed lagrangian methods for network resource allocation,” in _2017 IEEE Conference on Control Technology and Applications (CCTA)_ , Mauna Lani, HI, 2017, pp. 650–655.
  * [10] T. Yang, J. Lu, D. Wu, J. Wu, G. Shi, Z. Meng, and K. H. Johansson, “A distributed algorithm for economic dispatch over time-varying directed networks with delays,” _IEEE Transactions on Industrial Electronics_ , vol. 64, no. 6, pp. 5095–5106, 2017.
  * [11] J. Xu, S. Zhu, Y. C. Soh, and L. Xie, “A dual splitting approach for distributed resource allocation with regularization,” _IEEE Transactions on Control of Network Systems_ , vol. 6, no. 1, pp. 403–414, 2019\.
  * [12] A. Nedić, A. Olshevsky, and W. Shi, “Improved convergence rates for distributed resource allocation,” in _2018 IEEE Conference on Decision and Control (CDC)_ , Miami Beach, FL, 2018, pp. 172–177.
  * [13] G. Binetti, A. Davoudi, F. L. Lewis, D. Naso, and B. Turchiano, “Distributed consensus-based economic dispatch with transmission losses,” _IEEE Transactions on Power Systems_ , vol. 29, no. 4, pp. 1711–1720, 2014\.
  * [14] K. Kim, “Distributed learning algorithms and lossless convex relaxation for economic dispatch with transmission losses and capacity limits,” _Mathematical Problems in Engineering_ , vol. 2019, pp. 1–11, 2019.
  * [15] H. Xing, Y. Mou, M. Fu, and Z. Lin, “Distributed algorithm for economic power dispatch including transmission losses,” in _2015 European Control Conference (ECC)_ , 2015, pp. 1076–1081.
  * [16] C. Zhao, J. He, P. Cheng, and J. Chen, “Consensus-based energy management in smart grid with transmission losses and directed communication,” _IEEE Transactions on Smart Grid_ , vol. 8, no. 5, pp. 2049–2061, Sep. 2017.
  * [17] S. Lee and H. Shim, “Distributed algorithm for economic dispatch problem with separable losses,” _IEEE Control Systems Letters_ , vol. 3, no. 3, pp. 685–690, July 2019.
  * [18] R. Faranda, A. Pievatolo, and E. Tironi, “Load shedding: A new proposal,” _IEEE Transactions on Power Systems_ , vol. 22, no. 4, pp. 2086–2093, 2007.
  * [19] A. Hussain, V. Bui, and H. Kim, “An effort-based reward approach for allocating load shedding amount in networked microgrids using multiagent system,” _IEEE Transactions on Industrial Informatics_ , vol. 16, no. 4, pp. 2268–2279, 2020.
  * [20] T. Kato, H. Takahashi, K. Sasai, G. Kitagata, H.-M. Kim, and T. Kinoshita, “Priority-based hierarchical operational management for multiagent-based microgrids,” _Energies_ , vol. 7, no. 4, p. 2051–2078, Mar 2014.
  * [21] A. Hussain, V. Bui, and H. Kim, “A resilient and privacy-preserving energy management strategy for networked microgrids,” _IEEE Transactions on Smart Grid_ , vol. 9, no. 3, pp. 2127–2139, 2018.
  * [22] G. S. Grewal, J. W. Konowalec, and M. Hakim, “Optimization of a load shedding scheme,” _IEEE Industry Applications Magazine_ , vol. 4, no. 4, pp. 25–30, 1998.
  * [23] D. F. Teshome, P. F. Correia, and K. L. Lian, “Stochastic Optimization for Network-Constrained Power System Scheduling Problem,” _Mathematical Problems in Engineering_ , vol. 2015, 2015.
  * [24] F. L. Lewis, H. Zhang, and K. Hengster-Movric, _Cooperative Control of Multi-Agent Systems: Optimal and Adaptive Design Approaches_. Springer-Verlag, 2014.
  * [25] J. J. Grainger and W. D. Stevenson, _Power System Analysis_. New York: McGraw-Hill, 1994.
  * [26] N. Li, L. Chen, and S. H. Low, “Optimal demand response based on utility maximization in power networks,” in _2011 IEEE Power and Energy Society General Meeting_ , Detroit, MI, USA, 2011.
  * [27] D. P. Bertsekas and J. N. Tsitsiklis, _Parallel and Distributed Computation: Numerical Methods_. Englewood Cliffs, New Jersey: Prentice-Hall, 1989.
  * [28] A. Nedic, A. Ozdaglar, and P. A. Parrilo, “Constrained consensus and optimization in multi-agent networks,” _IEEE Transactions on Automatic Control_ , vol. 55, no. 4, pp. 922–938, 2010.
  * [29] A. Olshevsky, “Linear time average consensus and distributed optimization on fixed graphs,” _SIAM Journal on Control and Optimization_ , vol. 55, no. 6, pp. 3990–4014, 2017.
  * [30] L. Wang and C. Singh, “Balancing risk and cost in fuzzy economic dispatch including wind power penetration based on particle swarm optimization,” _Electric Power Systems Research_ , vol. 78, no. 8, pp. 1361 – 1368, 2008\.
  * [31] A. Nedić and A. Ozdaglar, “Approximate primal solutions and rate analysis for dual subgradient methods,” _SIAM Journal on Optimization_ , vol. 19, no. 4, pp. 1757–1780, 2009.

###  -A Proof of Theorem 1

Note that $\mathcal{L}(\bm{x},\bm{u},\lambda,\bm{\xi})$ given in (5) is the
relaxed Lagrangian function of (4) since the constraint
$x_{i}\in\mathcal{X}_{i}$ is considered implicitly. For the purpose of showing
(6), let us consider the full Lagrangian function of (4) as follows:

$\displaystyle\hat{\mathcal{L}}(\bm{x},\bm{u},\lambda,\bm{\xi},\bm{\varphi},\bm{\sigma})\\!\\!\\!$
$\displaystyle:=$
$\displaystyle\\!\\!\\!\sum_{i\in\mathcal{V}}\left(\sigma_{i}(x_{i}\\!-\\!x_{i}^{max})+\varphi_{i}(x_{i}^{min}-x_{i})\right)+\mathcal{L}(\bm{x},\bm{u},\lambda,\bm{\xi})$

where $\bm{\varphi}=[\varphi_{i}],\bm{\sigma}=[\sigma_{i}],$$i\in\mathcal{V}$.
The following is complementary slackness condition corresponding to inequality
constraint (4b)

$\displaystyle\accentset{\rule{2.79996pt}{0.8pt}}{\lambda}\geq
0,\quad\accentset{\rule{2.79996pt}{0.8pt}}{\lambda}\left(\sum\limits_{i\in\mathcal{V}}(\accentset{\rule{2.79996pt}{0.8pt}}{u}_{i}^{2}+d_{i}-\accentset{\rule{2.79996pt}{0.8pt}}{x}_{i})\right)=0.$
(48)

Let us consider the case when $\accentset{\rule{2.79996pt}{0.8pt}}{\lambda}=0$
and check the existence of $\accentset{\rule{2.79996pt}{0.8pt}}{x}_{i}$,
$\accentset{\rule{2.79996pt}{0.8pt}}{u}_{i}$,
$\accentset{\rule{2.79996pt}{0.8pt}}{\sigma}_{i}$, and
$\accentset{\rule{2.79996pt}{0.8pt}}{\varphi}_{i},i\in\mathcal{V},$ satisfying
the KKT conditions. At first, $\partial\hat{\mathcal{L}}/\partial x_{i}=0$ and
$\partial\hat{\mathcal{L}}/\partial u_{i}=0$ can be written as follows:

$\displaystyle\\!\\!\\!\partial\hat{\mathcal{L}}/\partial x_{i}$
$\displaystyle=$ $\displaystyle\nabla
C_{i}(\accentset{\rule{2.79996pt}{0.8pt}}{x}_{i})\\!-\\!\accentset{\rule{2.79996pt}{0.8pt}}{\lambda}+\\!\sum_{j\in\mathcal{V}}\accentset{\rule{2.79996pt}{0.8pt}}{\xi}^{(j)}r_{ji}\\!+\\!\accentset{\rule{2.79996pt}{0.8pt}}{\sigma}_{i}\\!-\\!\accentset{\rule{2.79996pt}{0.8pt}}{\varphi}_{i}\\!=\\!0,$
(49) $\displaystyle\\!\\!\\!\partial\hat{\mathcal{L}}/\partial u_{i}$
$\displaystyle=$ $\displaystyle
2\accentset{\rule{2.79996pt}{0.8pt}}{\lambda}u_{i}-\accentset{\rule{2.79996pt}{0.8pt}}{\xi}^{(i)}=0.$
(50)

Since $\accentset{\rule{2.79996pt}{0.8pt}}{\lambda}=0$, from (50) we have
$\accentset{\rule{2.79996pt}{0.8pt}}{\xi}^{(j)}=0$ for all $j\in\mathcal{V}$.
Substituting it to (49) yields $\nabla
C_{i}(\accentset{\rule{2.79996pt}{0.8pt}}{x}_{i})+\accentset{\rule{2.79996pt}{0.8pt}}{\sigma}_{i}-\accentset{\rule{2.79996pt}{0.8pt}}{\varphi}_{i}=0$.

Having this, let us consider the case where
$\accentset{\rule{2.79996pt}{0.8pt}}{\varphi}_{i}=0,\accentset{\rule{2.79996pt}{0.8pt}}{\sigma}_{i}=0$,
$i\in\mathcal{V}$. This means
$\accentset{\rule{2.79996pt}{0.8pt}}{x}_{i}\in\mbox{relint}(\mathcal{X}_{i})$
and yields $\nabla C_{i}(\accentset{\rule{2.79996pt}{0.8pt}}{x_{i}})=0$. Note
that this contradicts the assumption that $\nabla C_{i}(x_{i})>0$ for all
$x_{i}\in\mbox{relint}(\mathcal{X}_{i})$ described in Assumption 4. Next,
consider the case where $\accentset{\rule{2.79996pt}{0.8pt}}{\varphi}_{i}=0$
and $\accentset{\rule{2.79996pt}{0.8pt}}{\sigma}_{i}\not=0,$ meaning that
$\accentset{\rule{2.79996pt}{0.8pt}}{x}_{i}=x_{i}^{max}$ for all
$i\in\mathcal{V}$. It follows that
$C_{i}(\accentset{\rule{2.79996pt}{0.8pt}}{x}_{i})+\accentset{\rule{2.79996pt}{0.8pt}}{\sigma}_{i}=0$.
Since $\accentset{\rule{2.79996pt}{0.8pt}}{\sigma}_{i}>0$, it yields $\nabla
C_{i}(\accentset{\rule{2.79996pt}{0.8pt}}{x}_{i}^{max})<0$ which can not be
true according to Assumption 4. Furthermore, let us consider the case where
$\accentset{\rule{2.79996pt}{0.8pt}}{\varphi}_{i}\not=0$ and
$\accentset{\rule{2.79996pt}{0.8pt}}{\sigma}_{i}=0,$ meaning that
$\accentset{\rule{2.79996pt}{0.8pt}}{x}_{i}=x_{i}^{min}$, $i\in\mathcal{V}$.
In this case, we have $\nabla
C_{i}(\accentset{\rule{2.79996pt}{0.8pt}}{x}_{i}^{min})-\accentset{\rule{2.79996pt}{0.8pt}}{\varphi}_{i}=0$,
which can be true. In conclusion, we have
$\accentset{\rule{2.79996pt}{0.8pt}}{x}_{i}=x_{i}^{min}$ when $\lambda=0$. By
using (6c) and $B=R^{T}R$, we have
$\sum_{i\in\mathcal{V}}\accentset{\rule{2.79996pt}{0.8pt}}{u}_{i}^{2}=\Upsilon(\bm{x}_{i}^{min})$.The
following holds according to Assumption 5

$\displaystyle\sum\limits_{i\in\mathcal{V}}\accentset{\rule{2.79996pt}{0.8pt}}{u}_{i}^{2}+d_{i}-\accentset{\rule{2.79996pt}{0.8pt}}{x}_{i}=\Upsilon(\bm{x}_{i}^{min})+\sum\limits_{i\in\mathcal{V}}\left(d_{i}-x_{i}^{min}\right)>0.$

Note that this contradicts the assumption that
$\accentset{\rule{2.79996pt}{0.8pt}}{\lambda}=0$ which equivalent to
$\sum_{i\in\mathcal{V}}\accentset{\rule{2.79996pt}{0.8pt}}{u}_{i}^{2}+d_{i}-\accentset{\rule{2.79996pt}{0.8pt}}{x}_{i}<0$.
Owing to this, $\accentset{\rule{2.79996pt}{0.8pt}}{\lambda}>0$, meaning that
equality constraint (6b) holds.

###  -B Proof of Claim 1

Because $x_{i}(k+1)$ and $u_{i}(k+1)$ satisfy (11a), i.e.,
$[x_{i}(k\\!+\\!1),u_{i}(k\\!+\\!1)]\\!=\\!\inf_{{x_{i}\in\mathcal{X}_{i},u_{i}\in\mathcal{U}}}\\!\mathcal{L}_{i}(x_{i},u_{i},v_{i}(k),\bm{w}_{i}(k))$
, for any $k\geq 0,i\in\mathcal{V}$, we have

$\displaystyle 0$ $\displaystyle\leq$
$\displaystyle-\mathcal{L}_{i}(x_{i}(k\\!+\\!1),u_{i}(k\\!+\\!1),v_{i}(k),\bm{w}_{i}(k))$
(51)
$\displaystyle+\mathcal{L}_{i}(\accentset{\rule{2.79996pt}{0.8pt}}{x}_{i},\accentset{\rule{2.79996pt}{0.8pt}}{u}_{i},v_{i}(k),\bm{w}_{i}(k)).$

Due to the same reason,
$q_{i}(v_{i}(k),\bm{w}_{i}(k))=-h_{i}(v_{i}(k),\bm{w}_{i}(k))=-\inf_{x_{i}\in\mathcal{X}_{i},u_{i}\in\mathcal{U}}\mathcal{L}_{i}(x_{i},u_{i},v_{i}(k),\bm{w}_{i}(k))$
can be rewritten as

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!q_{i}(v_{i}(k),\bm{w}_{i}(k))\\!\\!=\\!-\mathcal{L}_{i}(x_{i}(k\\!+\\!1),\\!{u}_{i}(k\\!+\\!1),v_{i}(k),\bm{w}_{i}(k)).$
(52)

The zero-duality gap implies that
$\sum_{i\in\mathcal{V}}C_{i}(\accentset{\rule{2.79996pt}{0.8pt}}{x}_{i})=\sum_{i\in\mathcal{V}}h_{i}(\accentset{\rule{2.79996pt}{0.8pt}}{\lambda},\accentset{\rule{2.79996pt}{0.8pt}}{\bm{\xi}})=-\sum_{i\in\mathcal{V}}q_{i}(\accentset{\rule{2.79996pt}{0.8pt}}{\lambda},\accentset{\rule{2.79996pt}{0.8pt}}{\bm{\xi}})$.
As a result, we have

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\sum_{i\in\mathcal{V}}\mathcal{L}_{i}(\accentset{\rule{2.79996pt}{0.8pt}}{x}_{i},\accentset{\rule{2.79996pt}{0.8pt}}{u}_{i},v_{i}(k),\bm{w}_{i}(k))\\!=\\!$
$\displaystyle-\sum_{i\in\mathcal{V}}\\!q_{i}(\accentset{\rule{2.79996pt}{0.8pt}}{\lambda},\accentset{\rule{2.79996pt}{0.8pt}}{\bm{\xi}})\\!+\\!\sum_{i\in\mathcal{V}}v_{i}(k)(\accentset{\rule{2.79996pt}{0.8pt}}{u}_{i}^{2}\\!+\\!d_{i}\\!-\\!\accentset{\rule{2.79996pt}{0.8pt}}{x}_{i})\\!$
$\displaystyle+\\!\sum_{i\in\mathcal{V}}w_{i}^{(i)}(k)(r_{ii}\accentset{\rule{2.79996pt}{0.8pt}}{x}_{i}\\!-\\!\accentset{\rule{2.79996pt}{0.8pt}}{u}_{i})+\\!\sum_{i\in\mathcal{V}}\sum_{j\in\mathcal{V},j\not=i}\\!\\!\\!\\!\\!\\!w_{i}^{(j)}\\!(k)r_{ji}\accentset{\rule{2.79996pt}{0.8pt}}{x}_{i},\leavevmode\nobreak\
\leavevmode\nobreak\ $ (53)

where $\bm{w}_{i}=[w_{i}^{(j)}]$, $j\in\mathcal{V}.$ The following is obtained
from the summation of (51) over all $i\in\mathcal{V}$.

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!0\leq-\sum_{i\in\mathcal{V}}\mathcal{L}_{i}(x_{i}(k\\!+\\!1),u_{i}(k\\!+\\!1),v_{i}(k),\bm{w}_{i}(k))$
(54)
$\displaystyle+\sum_{i\in\mathcal{V}}\mathcal{L}_{i}(\accentset{\rule{2.79996pt}{0.8pt}}{x}_{i},\accentset{\rule{2.79996pt}{0.8pt}}{u}_{i},v_{i}(k),\bm{w}_{i}(k)).$

Substituting (52)-(-B) into (54) yields

$\displaystyle 0$ $\displaystyle\leq$
$\displaystyle\sum_{i\in\mathcal{V}}\\!q_{i}(v_{i}(k),\bm{w}_{i}(k))-\sum_{i\in\mathcal{V}}q_{i}(\accentset{\rule{2.79996pt}{0.8pt}}{\lambda},\accentset{\rule{2.79996pt}{0.8pt}}{\bm{\xi}})$
$\displaystyle\\!+\\!\sum_{i\in\mathcal{V}}v_{i}(k)(\accentset{\rule{2.79996pt}{0.8pt}}{u}_{i}^{2}\\!+\\!d_{i}\\!-\\!\accentset{\rule{2.79996pt}{0.8pt}}{x}_{i})\\!+\\!\sum_{i\in\mathcal{V}}w_{i}^{(i)}(k)(r_{ii}\accentset{\rule{2.79996pt}{0.8pt}}{x}_{i}\\!-\\!\accentset{\rule{2.79996pt}{0.8pt}}{u}_{i})$
$\displaystyle+\\!\sum_{i\in\mathcal{V}}\sum_{j\in\mathcal{V},j\not=i}\\!\\!\\!\\!\\!\\!w_{i}^{(j)}\\!(k)r_{ji}\accentset{\rule{2.79996pt}{0.8pt}}{x}_{i}.$

Since $q_{i}$ is a convex function, it follows that

$\displaystyle 0$ $\displaystyle\leq$
$\displaystyle\sum_{i\in\mathcal{V}}\partial
q_{i}(\accentset{\rule{2.79996pt}{0.8pt}}{\lambda},\accentset{\rule{2.79996pt}{0.8pt}}{\bm{\xi}})^{T}\left(\left[\begin{array}[]{c}\\!\\!\accentset{\rule{2.79996pt}{0.8pt}}{\lambda}\\!\\!\\\
\\!\\!\accentset{\rule{2.79996pt}{0.8pt}}{\bm{\xi}}\\!\\!\end{array}\right]\\!-\\!\left[\begin{array}[]{c}\\!\\!\\!\\!{v}_{i}(k)\\!\\!\\!\\!\\\
\\!\\!\\!\\!\bm{w}_{i}(k)\\!\\!\\!\\!\end{array}\right]\right)\\!$
$\displaystyle+\\!\sum_{i\in\mathcal{V}}v_{i}(k)(\accentset{\rule{2.79996pt}{0.8pt}}{u}_{i}^{2}\\!+\\!d_{i}\\!-\\!\accentset{\rule{2.79996pt}{0.8pt}}{x}_{i})\\!+\\!\sum_{i\i\mathcal{V}}w_{i}^{(i)}(k)(r_{ii}\accentset{\rule{2.79996pt}{0.8pt}}{x}_{i}\\!-\\!\accentset{\rule{2.79996pt}{0.8pt}}{u}_{i})$
$\displaystyle+\\!\sum_{i\i\mathcal{V}}\sum_{j\in\mathcal{V},j\not=i}\\!\\!\\!\\!\\!\\!w_{i}^{(j)}\\!(k)r_{ji}\accentset{\rule{2.79996pt}{0.8pt}}{x}_{i}$
$\displaystyle\leq$
$\displaystyle\sum_{i\in\mathcal{V}}\Gamma\left\|\left[\begin{array}[]{c}\accentset{\rule{2.79996pt}{0.8pt}}{\lambda}\\\
\accentset{\rule{2.79996pt}{0.8pt}}{\bm{\xi}}\end{array}\right]-\left[\begin{array}[]{c}{v}_{i}(k)\\\
\bm{w}_{i}(k)\end{array}\right]\right\|\\!$ (65)
$\displaystyle+\\!\sum_{i\in\mathcal{V}}v_{i}(k)(\accentset{\rule{2.79996pt}{0.8pt}}{u}_{i}^{2}\\!+\\!d_{i}\\!-\\!\accentset{\rule{2.79996pt}{0.8pt}}{x}_{i})\\!+\\!\sum_{i\in\mathcal{V}}w_{i}^{(i)}(k)(r_{ii}\accentset{\rule{2.79996pt}{0.8pt}}{x}_{i}\\!-\\!\accentset{\rule{2.79996pt}{0.8pt}}{u}_{i})$
$\displaystyle+\\!\sum_{i\in\mathcal{V}}\sum_{j\in\mathcal{V},j\not=i}\\!\\!\\!\\!\\!\\!w_{i}^{(j)}\\!(k)r_{ji}\accentset{\rule{2.79996pt}{0.8pt}}{x}_{i}.$

Note that the right hand side of (-B) goes to zero as $k\rightarrow\infty$
because
$(v_{i}(k),\bm{w}_{i}(k))\\!\\!\rightarrow\\!\\!(\accentset{\rule{2.79996pt}{0.8pt}}{\lambda},\accentset{\rule{2.79996pt}{0.8pt}}{\bm{\xi}})$
as $k\rightarrow\infty$,
$\sum_{i\in\mathcal{V}}(\accentset{\rule{2.79996pt}{0.8pt}}{u}_{i}^{2}+d_{i}-\accentset{\rule{2.79996pt}{0.8pt}}{x}_{i})\leq
0$, and
$\sum_{i\in\mathcal{V}}(\sum_{j\in\mathcal{V}}r_{ji}\accentset{\rule{2.79996pt}{0.8pt}}{x}_{i}-\accentset{\rule{2.79996pt}{0.8pt}}{u}_{i})=0$.
This implies
$\lim_{k\rightarrow\infty}\\!\sum_{i\in\mathcal{V}}\mathcal{L}_{i}(\accentset{\rule{2.79996pt}{0.8pt}}{x}_{i},\accentset{\rule{2.79996pt}{0.8pt}}{u}_{i},v_{i}(k),\bm{w}_{i}(k))\\!-\\!\sum_{i\in\mathcal{V}}\mathcal{L}_{i}(x_{i}(k\\!+\\!1),u_{i}(k\\!+\\!1),v_{i}(k),\bm{w}_{i}(k))\\!=\\!0$
which is equivalent to (30).

###  -C Proof of Theorem 3

Given that $y_{i}(k)\in\mathcal{Y}_{i}$ and $z_{i}(k)\in\mathcal{Z}_{i}$,
there exists $\Lambda>0$ such
$\|\bm{\mathsf{g}}_{i}(z_{i}(k),y_{i}(k))\|\leq\Lambda$ for all $k\geq 0$.
Thus,
$\bm{\eta}_{i}(k)\rightarrow\accentset{\rule{2.79996pt}{0.8pt}}{\bm{\eta}}$ as
$k\rightarrow\infty$ according to [31, Proposition 5]. Using the similar
reasoning for obtaining (-B), we have

$\displaystyle 0\\!\\!\\!\\!\\!$ $\displaystyle\\!\\!\\!\\!\leq\\!\\!\\!\\!$
$\displaystyle\\!\\!\\!\\!\sum_{i\in\mathcal{V}}\mathcal{L}_{i}^{ls}(\accentset{\rule{2.79996pt}{0.8pt}}{y}_{i},\accentset{\rule{2.79996pt}{0.8pt}}{z}_{i},\bm{\phi}_{i}(k))\\!-\\!\sum_{i\in\mathcal{V}}\mathcal{L}_{i}^{ls}(y_{i}(k\\!+\\!1),z_{i}(k\\!+\\!1),\bm{\phi}_{i}(k))$
(66) $\displaystyle\leq$
$\displaystyle\sum_{i\in\mathcal{V}}\left(\Lambda\left\|\accentset{\rule{2.79996pt}{0.8pt}}{\bm{\eta}}-\bm{\phi}_{i}(k)\right\|+\sum_{\nu=1}^{m+1}\phi_{i}^{(\nu)}(k)\mathsf{g}_{i}^{(\nu)}(\accentset{\rule{2.79996pt}{0.8pt}}{z}_{i},\accentset{\rule{2.79996pt}{0.8pt}}{y}_{i})\right).$

Since $F_{i}$ is a strongly convex function, it follows that
$\mathcal{L}_{i}^{ls}(y_{i},z_{i},\bm{\phi}_{i}(k))$ is a strongly convex
function for any $\bm{\phi}_{i}(k)$. In other words, there exists
$\gamma_{i}>0$ such that

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\frac{\gamma_{i}}{2}\left\|\left[\begin{array}[]{c}\accentset{\rule{2.79996pt}{0.8pt}}{y}_{i}\\\
\accentset{\rule{2.79996pt}{0.8pt}}{z}_{i}\end{array}\right]-\left[\begin{array}[]{c}{y}_{i}(k+1)\\\
z_{i}(k+1)\end{array}\right]\right\|^{2}\leq\mathcal{L}^{ls}_{i}(\accentset{\rule{2.79996pt}{0.8pt}}{y}_{i},\accentset{\rule{2.79996pt}{0.8pt}}{z}_{i},\accentset{\rule{2.79996pt}{0.8pt}}{\lambda},\bm{\phi}_{i}(k))\\!$
(71)
$\displaystyle\\!\\!\\!\\!\\!\\!-\mathcal{L}_{i}^{ls}(y_{i}(k\\!+\\!1),z_{i}(k\\!+\\!1),\bm{\phi}_{i}(k))$
$\displaystyle\\!\\!\\!\\!\\!\\!-\nabla\mathcal{L}_{i}^{ls}(y_{i}(k\\!+\\!1),z_{i}(k\\!+\\!1),\bm{\phi}_{i}(k))^{T}\\!\\!\left(\left[\begin{array}[]{c}\\!\\!\\!\accentset{\rule{2.79996pt}{0.8pt}}{y}_{i}\\!\\!\\!\\\
\\!\\!\\!\accentset{\rule{2.79996pt}{0.8pt}}{z}_{i}\\!\\!\\!\end{array}\right]\\!-\\!\left[\begin{array}[]{c}\\!\\!\\!\\!y_{i}(k\\!+\\!1)\\!\\!\\!\\!\\\
\\!\\!\\!\\!z_{i}(k\\!+\\!1)\\!\\!\\!\\!\end{array}\right]\right)\\!.$ (76)

By definition,
$\nabla\mathcal{L}_{i}^{ls}(y_{i}(k\\!+\\!1),z_{i}(k\\!+\\!1),\bm{\phi}_{i}(k))^{T}=0$.
Summing (71) over $i\in\mathcal{V}$ and using (66) yields

$\displaystyle\sum_{i\in\mathcal{V}}\frac{\gamma_{i}}{2}\left\|\left[\begin{array}[]{c}\accentset{\rule{2.79996pt}{0.8pt}}{y}_{i}\\\
\accentset{\rule{2.79996pt}{0.8pt}}{z}_{i}\end{array}\right]-\left[\begin{array}[]{c}{y}_{i}(k+1)\\\
z_{i}(k+1)\end{array}\right]\right\|^{2}\\!\\!\\!$
$\displaystyle\\!\\!\\!\leq\\!\\!\\!$
$\displaystyle\\!\\!\\!\sum_{i\in\mathcal{V}}\Lambda\left\|\accentset{\rule{2.79996pt}{0.8pt}}{\bm{\eta}}-\bm{\phi}_{i}(k)\right\|$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!+\sum_{i\in\mathcal{V}}\sum_{\nu=1}^{m+1}\phi_{i}^{(\nu)}(k)\mathsf{g}_{i}^{(\nu)}(\accentset{\rule{2.79996pt}{0.8pt}}{z}_{i},\accentset{\rule{2.79996pt}{0.8pt}}{y}_{i}).$

Given that the right-hand side of this inequality goes to zero as
$k\rightarrow\infty$ according to (66), we have
$(z_{i}(k),y_{i}(k))\rightarrow(\accentset{\rule{2.79996pt}{0.8pt}}{z}_{i},\accentset{\rule{2.79996pt}{0.8pt}}{y}_{i})$
as $k\rightarrow\infty$.