Accelerated Voltage Regulation in Multi-Phase Distribution Networks Based on Hierarchical Distributed Algorithm

TL;DR

This paper introduces a hierarchical distributed algorithm for voltage regulation in large multi-phase distribution networks, significantly improving computational efficiency and scalability while maintaining optimality.

Contribution

It presents a novel hierarchical distributed primal-dual gradient algorithm leveraging network structure and linearized power flow models for scalable OPF solutions.

Findings

Achieves over 10-fold acceleration in convergence speed.

Reduces computational loads compared to centralized methods.

Demonstrates effectiveness on a 4,521-node test feeder.

Abstract

We propose a hierarchical distributed algorithm to solve optimal power flow (OPF) problems that aim at dispatching controllable distributed energy resources (DERs) for voltage regulation at minimum cost. The proposed algorithm features unprecedented scalability to large multi-phase distribution networks by jointly exploring the tree/subtrees structure of a large radial distribution network and the structure of the linearized distribution power flow (LinDistFlow) model to derive a hierarchical, distributed implementation of the primal-dual gradient algorithm that solves OPF. The proposed implementation significantly reduces the computation loads compared to the centrally coordinated implementation of the same primal-dual algorithm without compromising optimality. Numerical results on a 4,521-node test feeder show that the designed algorithm achieves more than 10-fold acceleration in the speed of convergence compared to the centrally coordinated primal-dual algorithm through reducing and distributing computational loads.

Figures

Equations

$$P_{ij}$$

$$Q_{ij}$$

$$v_j$$

$${\ell }_{ij}{v}_i$$

$$\mathbf{v}$$

$$R_{ij}:=\sum _{(\zeta ,\xi )\in {\mathcal{E}}_i\cap {\mathcal{E}}_j}2\cdot r_{\zeta \xi },X_{ij}:=\sum _{(\zeta ,\xi )\in {\mathcal{E}}_i\cap {\mathcal{E}}_j}2\cdot x_{\zeta \xi }.$$

$$\displaystyle{\cal Y}_{i}=\left\{(p_{i},q_{i})\big{|}0\leq p_{i}\leq p_{i}^{\text{av}},~{}p_{i}^{2}+q_{i}^{2}\leq\eta_{i}^{2}\right\},$$

$$\displaystyle{\cal Y}_{i}=\left\{(p_{i},q_{i})\big{|}\underline{p}_{i}\leq p_{i}\leq\overline{p}_{i},~{}p_{i}^{2}+q_{i}^{2}\leq\eta_{i}^{2}\right\},$$

$$\displaystyle{\cal Y}_{i}=\left\{(p_{i},q_{i})\big{|}\underline{p}_{i}\leq p_{i}\leq\overline{p}_{i},\underline{q}_{i}\leq q_{i}\leq\overline{q}_{i}\right\}.$$

$$P_0=-P_I-\sum _{i\in \mathcal{N}}{p}_i,$$

$$\underset{\mathbf{p},\mathbf{q}}{\min }$$

$$\underline{\mathbf{v}}\le \mathbf{v}(\mathbf{p},\mathbf{q})\le \overline{\mathbf{v}},$$

$$\mathcal{L}(\mathbf{p},\mathbf{q};\overline{\mathbf{\mu }},\underline{\mathbf{\mu }})=\sum _{i\in \mathcal{N}}{C}_i(p_i,q_i)+C_0(P_0(\mathbf{p}))$$

$$+{\underline{\mathbf{\mu }}}^{\top }(\underline{\mathbf{v}}-\mathbf{v}(\mathbf{p},\mathbf{q}))+{\overline{\mathbf{\mu }}}^{\top }(\mathbf{v}(\mathbf{p},\mathbf{q})-\overline{\mathbf{v}}),$$

$${\mathcal{L}}_{\eta }(\mathbf{p},\mathbf{q};\overline{\mathbf{\mu }},\underline{\mathbf{\mu }})=\sum _{i\in \mathcal{N}}{C}_i(p_i,q_i)+C_0(P_0(\mathbf{p}))$$

$$+{\underline{\mathbf{\mu }}}^{\top }(\underline{\mathbf{v}}-\mathbf{v}(\mathbf{p},\mathbf{q}))+{\overline{\mathbf{\mu }}}^{\top }(\mathbf{v}(\mathbf{p},\mathbf{q})-\overline{\mathbf{v}})-\frac{\eta }{2}\Vert \mathbf{\mu }{\Vert }_2^2.$$

$$\mathbf{p}(t+1)$$

$$\mathbf{q}(t+1)$$

$$\underline{\mathbf{\mu }}(t+1)$$

$$\overline{\mathbf{\mu }}(t+1)$$

$$\mathbf{v}(t+1)$$

$$P_0(t+1)$$

$$\left[\begin{array}{c}\mathbf{y}(t+1)\\ \mathbf{\mu }(t+1)\end{array}\right]={\left[\left[\begin{array}{c}\mathbf{y}(t)\\ \mathbf{\mu }(t)\end{array}\right]-ϵ\left[\begin{array}{c}{\nabla }_{\mathbf{y}}{\mathcal{L}}_{\eta }(\mathbf{y}(t),\mathbf{\mu }(t))\\ -{\nabla }_{\mathbf{\mu }}{\mathcal{L}}_{\eta }(\mathbf{y}(t),\mathbf{\mu }(t))\end{array}\right]\right]}_{\mathcal{Y}\times \mathcal{U}},$$

$$(T(\mathbf{z})-T({\mathbf{z}}^{\prime }){)}^{\top }(\mathbf{z}-{\mathbf{z}}^{\prime })\ge M\Vert \mathbf{z}-{\mathbf{z}}^{\prime }{\Vert }_2^2.$$

$$\Vert T(\mathbf{z})-T({\mathbf{z}}^{\prime }){\Vert }_2^2\le L^2\Vert \mathbf{z}-{\mathbf{z}}^{\prime }{\Vert }_2^2.$$

$$(T(\mathbf{z})-T({\mathbf{z}}^{\prime }){)}^{\top }(\mathbf{z}-{\mathbf{z}}^{\prime })$$

$$0<ϵ\le \overline{ϵ}<2M/L^2$$

$$\displaystyle p_{i}(t+1)=\Big{[}p_{i}(t)-\epsilon\Big{(}\partial_{p_{i}}C_{i}(p_{i}(t),q_{i}(t))-C^{\prime}_{0}(P_{0}(t))$$

$$\displaystyle\hskip 28.45274pt+\sum_{j\in{\cal N}}R_{ij}\big{(}\overline{\mu}_{j}(t)-\underline{\mu}_{j}(t)\big{)}\Big{)}\Big{]}_{{\cal Y}_{i}},\ i\in{\cal N}.$$

$$\sum _{j\in {\mathcal{N}}_k}{R}_{ij}({\overline{\mu }}_j-{\underline{\mu }}_j)+\sum _{j\in \mathcal{N}\backslash {\mathcal{N}}_k}{R}_{ij}({\overline{\mu }}_j-{\underline{\mu }}_j)$$

$$\sum _{k\in \mathcal{K}}{R}_{i{n}_k^0}\sum _{j\in {\mathcal{N}}_k}({\overline{\mu }}_j-{\underline{\mu }}_j)+\sum _{j\in {\mathcal{N}}_0}{R}_{ij}({\overline{\mu }}_j-{\underline{\mu }}_j),$$

$$z_{\zeta \xi }=\left[\begin{array}{ccc}{z}_{\zeta \xi }^{aa}& z_{\zeta \xi }^{ab}& z_{\zeta \xi }^{ac}\\ z_{\zeta \xi }^{ba}& z_{\zeta \xi }^{bb}& z_{\zeta \xi }^{bc}\\ z_{\zeta \xi }^{ca}& z_{\zeta \xi }^{cb}& z_{\zeta \xi }^{cc}\end{array}\right]\in {\mathbb{C}}^{3\times 3},$$

$$Z_{ij}^{\phi \varphi }=\sum _{(\zeta ,\xi )\in {\mathcal{E}}_i\cap {\mathcal{E}}_j}{z}_{\zeta \xi }^{\phi \varphi }\in \mathbb{C}$$

$${\mathbf{v}}_{{}_{\Xi }}=R_{{}_{\Xi }}{\mathbf{p}}_{{}_{\Xi }}+X_{{}_{\Xi }}{\mathbf{q}}_{{}_{\Xi }}+{\tilde{\mathbf{v}}}_{{}_{\Xi }},$$

$${\partial }_{p_j^{\varphi }}{v}_i^{\phi }$$

$${\partial }_{q_j^{\varphi }}{v}_i^{\phi }$$

$$P_0=-\sum _{\varphi \in {\Phi }_0}{P}_I^{\varphi }-\sum _{i\in \mathcal{N}}\sum _{\varphi \in {\Phi }_i}{p}_i^{\varphi },$$

Full text

I Introduction

OPF problems determine the best operating points of dispatchable devices in electric power grids and achieve the optimal system-wide objectives and operational constraints for important applications such as demand response and voltage regulation. However, the increasing penetrations of DERs such as roof-top photovoltaic, electric vehicles, battery energy storage systems, thermostatically controlled loads, and other controllable loads not only provide enormous potential optimization and control flexibility that we can explore [2], but also make OPF more challenging to solve with significantly growing dimensionality. Meanwhile, due to the intermittent nature of the renewable energy resources, their deepening penetration in the distribution networks causes large and rapid fluctuations in power injections and voltages, and calls for fast control paradigms.

Distributed algorithms are developed to facilitate scalable and fast control of large networks of dispatchable DERs by parallel and distributed computation. In literature, such algorithms have been designed and implemented either with a central coordinator (CC), e.g., [3, 4, 5, 6, 7, 8], or among neighboring agents without a CC, e.g., [9, 10, 11, 12, 13]. The latter usually demands that all nodes in the network compute and pass along updated information to their neighbors, which may not be implementable in practice if not all nodes are controllable or equipped with computation and communication capability. This paper will focus on the former communication model, which only requires CC to communicate with controllable nodes. However, what has been overlooked by existing literature on centrally coordinated distributed algorithms is that, while part of the computation loads are distributed among control nodes, the remaining that are executed by CC may surge as the size of the problem increases. Take voltage regulation problems (to be elaborated later) as an example: let $N$ denote the number of networked nodes, and the computational complexity of CC is in the order of $N^{2}$. When $N$ grows, the computational time of CC becomes exceedingly notable, rendering CC the bottleneck and hindering fast algorithm implementation in large systems.

One way to mitigate the computational loads of CC (and therefore accelerate the algorithms) is to consider multiple coordinators among which CC’s loads can be distributed. However, due to the complexity of the coupling term calculated by CC, it is usually challenging to decompose it in an efficient and yet structurally meaningful way. To our best knowledge, little work has been done to explore this area. In particular, hierarchical structures have been under-exploited in developing distributed OPF algorithms.

Moreover, realistic distribution systems are usually featured with multi-phase unbalanced loads. The resultant multi-phase OPF has been studied extensively in literature, e.g., [14, 15, 16, 17, 13, 18]. Indeed, the inter-phase coupling in multi-phase systems adds extra difficulty to the already complex OPF problems and makes computational load reduction and distribution even more challenging.

This work considers a large multi-phase distribution network with tree topology. We consider the linearized distribution flow (LinDistFlow) model [19, 20, 21, 22, 23] as well as its multi-phase extension [13, 17] for the distribution networks. The network is divided into areas featuring subtree topology. A regional coordinator (RC) communicates with all the dispatchable nodes within each subtree, and CC communicates with all RCs. Each RC knows only the topology and line parameters of the subtree that it coordinates, and CC knows only the topology and line parameters of the reduced network which treats each subtree as a node. Given such information availability, we explore the topological structure of the LinDistFlow model to derive a hierarchical, distributed implementation of the primal-dual gradient algorithm that solves an OPF problem. The OPF problem minimizes the total cost over all the controllable DERs and a cost associated with the total network load subject to voltage regulation constraints. The proposed implementation significantly reduces the computation burden of CC compared to the centrally coordinated implementation of the primal-dual algorithm, not only by distributing computational loads among RCs but also by reducing repetitive calculations and information transfers. Convergence of the designed algorithm is guaranteed with unbalanced nonlinear power flow.

Performance of the proposed implementation is verified through numerical simulation of a three-phase unbalanced 4,521-node test feeder with 1,043 controllable nodes. Simulation results show that a 10-fold acceleration in the speed of convergence can be achieved by the hierarchical distributed method compared to the centrally coordinated implementation. This significant improvement in convergence speed makes real-time grid optimization and control possible. Meanwhile, to our best knowledge, the size of the network in our simulation is the largest in distributed OPF studies.

It is worth noting that the model and algorithm design in this work can be readily applied to optimization and control of networked microgrids with each subtree seen as a microgrid. Meanwhile, most existing works on optimization and control of networked microgrids either over-simplify each individual microgrid as a node without considering power flow within it [24, 25, 26], or ignore power flow models among microgrids [27, 28, 29]. Recent work [30] applies a game-theoretic approach to manage a partitioned distribution network based on a noncooperative Nash game, where uniqueness of equilibrium, convergence, and global performance are difficult to characterize. Different from those in the literature, our method models inter- and intra-microgrid power flow dynamics and features provable convergence and global optimality performance.

Also note that, another line of works on voltage regulation are based on local control algorithms, e.g., [31, 32, 22, 23] and references therein. Though local voltage control has the advantages of low computation and communication complexity, it can only solve a limited class of problems and cannot address more general objective functions and constraints as formulated in this work.

The rest of the paper is organized as follows. To provide design intuition, we start with Section II to model the single-phase distribution system and formulate the OPF problem, and follow with Section III to propose a hierarchical distributed implementation of the primal-dual gradient algorithm. We then in Section IV introduce the multi-phase distribution system and its OPF for which we elaborate a hierarchical distributed algorithm and characterize its convergence with nonlinear power flow. Section V presents numerical results, and Section VI concludes this paper.

II Solving OPF in Single-Phase System

II-A Single-Phase Power Flow Model

Consider a radial single-phase power distribution network denoted by ${\cal T}=\{{\cal N}\cup\{0\},{\cal E}\}$ with $N+1$ nodes collected in the set ${\cal N}\cup\{0\}$ where ${\cal N}:=\{1,...,N\}$ and node [math] is the slack bus, and distribution lines collected in the set ${\cal E}$. For each node $i\in\mathcal{N}$, denote by ${\cal E}_{i}\subseteq{\cal E}$ the set of lines on the unique path from node [math] to node $i$, and let $p_{i}$ and $q_{i}$ denote the real and reactive power injected, where negative (resp. positive) power injection means power consumption (resp. generation). Let $v_{i}$ be the squared magnitude of the complex voltage (phasor) at node $i$. For each line $(i,j)\in{\cal E}$, denote by $r_{ij}$ and $x_{ij}$ its resistance and reactance, and $P_{ij}$ and $Q_{ij}$ the real and reactive power from node $i$ to node $j$. Let $\ell_{ij}$ denote the squared magnitude of the complex branch current (phasor) from node $i$ to $j$.

We adopt the following DistFlow model [19, 20] for the radial distribution network:

[TABLE]

Following [21, 22] we assume that the active and reactive power loss $r_{ij}\ell_{ij}$ and $x_{ij}\ell_{ij}$, as well as $r^{2}_{ij}\ell_{ij}$ and $x^{2}_{ij}\ell_{ij}$, are negligible and can thus be ignored. Indeed, the losses are much smaller than power flows $P_{ij}$ and $Q_{ij}$, typically on the order of $1\%$. With the above approximations (1) is simplified to the following linear model:

[TABLE]

where bold symbols $\bm{v}=[v_{1},\ldots,v_{N}]^{\top}$, $\bm{p}=[p_{1},\ldots,p_{N}]^{\top}$, $\bm{q}=[q_{1},\ldots,q_{N}]^{\top}\in\mathbb{R}^{N}$ represent vectors, $\tilde{\bm{v}}$ is a constant vector with every component being the squared voltage magnitude at the slack bus, and the sensitivity matrices $R,X\in\mathbb{R}_{+}^{N\times N}$ respectively contain elements of

[TABLE]

Here, the voltage-to-power-injection sensitivity factors $R_{ij}$ (resp. $X_{ij}$) represents the resistance (resp. reactance) of the common path of node $i$ and $j$ leading back to node 0. Keep in mind that this result serves as the basis for designing the hierarchical distributed algorithm to be introduced later. Fig. 1 (left) illustrates the common path ${\cal E}_{i}\cap{\cal E}_{j}$ for two arbitrary nodes $i$ and $j$ in a radial network and their $R_{ij},X_{ij}$.

II-B OPF and Primal-Dual Gradient Algorithm

Assume node $i\in{\cal N}$ has a dispatchable DER (or aggregation of DERs) whose real and reactive power injections are confined by $(p_{i},q_{i})\in{\cal Y}_{i}$ where ${\cal Y}_{i}$ is a convex and compact set. For example, the feasible region of a PV inverter can be modeled as:

[TABLE]

where $p_{i}^{\text{av}}$ denotes the available active power from a PV system (based on prevailing ambient conditions), and $\eta_{i}$ is the rated apparent capacity. The feasible set for an energy storage system can be modeled as

[TABLE]

for given limits $\underline{p}_{i},\overline{p}_{i}$ and for a given inverter capacity rating $\eta_{i}$. The limits $\underline{p}_{i},\overline{p}_{i}$ are updated during the operation of the battery based on the state of charge. The operating region of small-scale diesel generators can be modeled using constant box constraints

[TABLE]

Let $P_{0}$ denote active power injected into the feeder at node 0, which is linearly approximated by the total active power loads as

[TABLE]

where $P_{I}$ denotes the total uncontrollable (inelastic) power injection in the network. Use $\bm{v}(\bm{p},\bm{q})$ and $P_{0}(\bm{p})$ to represent (2) and (4), respectively, and consider the following OPF problem:

[TABLE]

where the objective $C_{i}(p_{i},q_{i})$ is the cost function for node $i$ and the coupling term $C_{0}(P_{0}(\bm{p}))$ represents the cost associated with the total network load. For example, $C_{0}(P_{0}(\bm{p}))=\alpha(P_{0}(\bm{p})-\tilde{P}_{0})^{2}$ penalizes $P_{0}(\bm{p})$’s deviation from a dispatching signal $\tilde{P}_{0}$ from the bulk system operator, where $\alpha>0$ is a given weighting factor. We make the following assumption for these cost functions.

Assumption 1**.**

$C_{i}(p_{i},q_{i}),\>\forall i\in{\cal N}$* are continuously differentiable and strongly convex in $(p_{i},q_{i})$, with bounded first-order derivative in ${\cal Y}_{i}$; $C_{0}(P_{0})$ is continuously differentiable and convex with bounded first-order derivative. *

Associate dual variables $\underline{\bm{\mu}}$ and $\overline{\bm{\mu}}$ with the left-hand-side and the right-hand-side of (5c), respectively, to write the Lagrangian of (5) as:

[TABLE]

with (5c) treated as the domain of $(\bm{p},\bm{q})$.

In order to design an algorithm with provable convergence, we introduce the following regularized Lagrangian with parameter $\eta>0$ and $\bm{\mu}:=[\underline{\bm{\mu}}^{\top},\overline{\bm{\mu}}^{\top}]^{\top}$:

[TABLE]

Since ${\cal{L}}_{\eta}(\bm{p},\bm{q};\underline{\bm{\mu}},\overline{\bm{\mu}})$ is strongly convex in $\bm{p},\bm{q}$ and strongly concave in $\underline{\bm{\mu}},\overline{\bm{\mu}}$, the next result follows.

Theorem 1**.**

*There exists one unique saddle point $(\bm{p}^{*},\bm{q}^{*};\underline{\bm{\mu}}^{*},\overline{\bm{\mu}}^{*})$ of ${\cal{L}}_{\eta}$. *

Furthermore, the discrepancy due to the regularization term can be bounded and is proportional to $\eta$. We refer the details to [33, 7]. For the rest of this section, we will focus on solving the saddle point of the regularized Lagrangian (7). Specifically, we cast the iterative projected primal-dual gradient algorithm to find the saddle point of (7) as follows:

[TABLE]

where $\epsilon>0$ is a constant stepsize to be determined, $\bm{1}_{N}=[1,\ldots,1]^{\top}\in\mathbb{R}^{N}$, $[\ ]_{\bm{{\cal Y}}}$ is the projection operator onto feasible sets $\bm{{\cal Y}}:=\bigtimes_{i\in{\cal N}}{\cal Y}_{i}$, and $[\ ]_{+}$ is the projection operator onto the positive orthant.

II-C Convergence Analysis

For uncluttered notation, we use $\bm{y}:=[\bm{p}^{\top},\bm{q}^{\top}]^{\top}$ to stack the primal variables and equivalently rewrite (8) as follows:

[TABLE]

where $\bm{{\cal U}}$ denotes the feasible positive orthant for the dual variables. We further let $\bm{z}:=[\bm{y}^{\top},\bm{\mu}^{\top}]^{\top}$ stack all variables, and define the gradient operator $T(\bm{z}):=\begin{bmatrix}\nabla_{\bm{y}}\mathcal{L}_{\eta}(\bm{y},\bm{\mu})\\ -\nabla_{\bm{\mu}}\mathcal{L}_{\eta}(\bm{y},\bm{\mu})\end{bmatrix}$.

Lemma 1**.**

$T(\bm{z})$* is a strongly monotone operator. *

We refer to Appendix for the proof.

By Lemma 1, there exists some constant $M>0$ such that for any $\bm{z},\bm{z^{\prime}}\in\bm{{\cal Y}}\times\bm{{\cal U}}$, one has

[TABLE]

Moreover, based on Assumption 1 the operator $T(\bm{z})$ is also Lipschitz continuous, i.e., there exists some constant $L>0$ such that for any $\bm{z},\bm{z^{\prime}}\in\bm{{\cal Y}}\times\bm{{\cal U}}$, we have

[TABLE]

Lemma 2**.**

*The following relation holds: $M\leq L$. *

Proof.

Apply (11) to have

[TABLE]

Combine (10) and (12) to obtain the result.

Based on the results established so far, we present the next theorem that guarantees the convergence of the primal-dual gradient algorithm (9) with small enough stepsize.

Theorem 2**.**

If the stepsize $\epsilon$ satisfies

[TABLE]

*for some $\overline{\epsilon}$, (9) converges to the unique saddle point of (7) exponentially fast. *

We refer to Appendix for the proof.

II-D Motivation for Hierarchical Design

Note that in (8) the update of any $p_{i}$ (resp. $q_{i}$) involves the knowledge of $\sum_{j\in{\cal N}}R_{ij}(\overline{\mu}_{j}-\underline{\mu}_{j})$ (resp. $\sum_{j\in{\cal N}}X_{ij}(\overline{\mu}_{j}-\underline{\mu}_{j})$). Therefore, at each iteration a coordinator cognizant of the entire network’s sensitivity matrices $R$ and $X$ is required to collect updated dual variables from all nodes, calculate $R^{\top}(\overline{\bm{\mu}}-\underline{\bm{\mu}})$ and $X^{\top}(\overline{\bm{\mu}}-\underline{\bm{\mu}})$, and send back the corresponding result to every node. This becomes computationally more challenging in a larger network containing thousands of or even more controllable endpoints, not to mention the re-calculation of the large $R,X$ matrices in case changes in network topology or regulator taps occur.

This motivates us to design a hierarchical control structure where the large network is partitioned into smaller subtrees, each managed locally by its regional coordinator (RC), and there is a central coordinator (CC) that only manages a reduced network where each subtree is treated as one node. As we will show in next section, the hierarchical control not only distributes but also reduces a large amount of computation overall.

III Hierarchical Distributed Algorithm

In this section, we explore the tree structure of the distribution network as well as the construction of its sensitivity matrices to facilitate the design of an efficient hierarchical distributed algorithm. To this end, we formally define subtree as follows.

Definition 1**.**

*A subtree of a tree ${\cal T}$ is a tree consisting of a node in ${\cal T}$, all its descendants in ${\cal T}$, and their connecting lines. *

We categorize all nodes of distribution network ${\cal T}$ into two groups:

1. $K$ subtrees indexed by ${\cal T}_{k}=\{{\cal N}_{k},{\cal E}_{k}\},\ k\in{\cal K}=\{1,\ldots,K\}$, and 2) a set ${\cal N}_{0}$ collecting all the other “unclustered” nodes in ${\cal N}$. Here, ${\cal N}_{k}$ of size $N_{k}$ is the set of nodes in subtree ${\cal T}_{k}$ and ${\cal E}_{k}$ contains their connecting lines. Thus we have $\cup_{k\in{\cal K}}{\cal N}_{k}\cup{\cal N}_{0}={\cal N}\cup\{0\}$ and ${\cal N}_{j}\cap{\cal N}_{k}=\emptyset,\forall j\neq k$. Assume each subtree ${\cal T}_{k}$ is managed by an RC cognizant of the topology of ${\cal T}_{k}$ and communicating with all the controllable nodes within ${\cal T}_{k}$.

Denote the root node of subtree ${\cal T}_{k}$ by $n_{k}^{0}$, and consider a reduced network ${\cal T}^{r}=\{{\cal N}^{r}\cup\{0\},{\cal E}^{r}\}$ where ${\cal N}^{r}:=\cup_{k\in{\cal K}}\{n_{k}^{0}\}\cup{\cal N}_{0}$ consists of the root nodes of all subtrees and all the unclustered nodes, and ${\cal E}^{r}$ is the set of their connecting lines. We assume there to be a CC cognizant of the topology of the reduced network ${\cal T}^{r}$ and communicating with all the RCs as well as the unclustered nodes.

III-A Hierarchical Distributed Algorithm

For simplicity, we elaborate the algorithm design for real power injections $\bm{p}$ only, and that for $\bm{q}$ follows similarly.

We equivalently rewrite (8a) as:

[TABLE]

Note that in (14) while $\partial_{p_{i}}C_{i}(p_{i}(t),q_{i}(t))$ is local information and the scalar derivative $C^{\prime}_{0}(P_{0}(t))$ can be easily broadcast, the last term $\sum_{j\in{\cal N}}R_{ij}(\overline{\mu}_{j}(t)-\underline{\mu}_{j}(t))$ couples the entire network. How to efficiently compute the coupling term is the key to a scalable algorithm. For that purpose, we introduce the following lemma.

Lemma 3**.**

*Given any two subtrees ${\cal T}_{h}$ and ${\cal T}_{k}$ with their root nodes $n^{0}_{h}$ and $n^{0}_{k}$, we have $R_{ij}=R_{n^{0}_{h}n^{0}_{k}},X_{ij}=X_{n^{0}_{h}n^{0}_{k}}$, for any $i\in{\cal N}_{h}$, and any $j\in{\cal N}_{k}$. Similarly, given any unclustered node $i\in{\cal N}_{0}$ and a subtree ${\cal T}_{k}$ with its root node $n^{0}_{k}$, we have $R_{ij}=R_{in^{0}_{k}},X_{ij}=X_{in^{0}_{k}}$, for any $j\in{\cal N}_{k}$. *

Proof.

We have the two following facts, which are also illustrated in Fig. 1. First, by (3), $R_{ij}$ (resp. $X_{ij}$) is the summed resistance (resp. reactance) on the common path of node $i$ and $j$ leading back to node 0. Second, any node in one subtree and any node in another subtree (or any node in one subtree and one unclustered node) share the same common path back to node 0. The result follows immediately.

Lemma 3 indicates that the sensitivity matrices $R$ and $X$ possess block structure such that all nodes within one subtree ${\cal N}_{h}$ share the same sensitivity with respect to all nodes within another subtree ${\cal N}_{k}$; see Fig. 2 for illustration. This permits a hierarchical distributed way to recalculate the coupling terms.

For clustered node $i\in{\cal N}_{k}$, we decompose $\sum_{j\in{\cal N}}R_{ij}(\overline{\mu}_{j}-\underline{\mu}_{j})$ as:

[TABLE]

where the first term of (15) consists of information within subtree $k$ (together with the line parameter from $n_{k}^{0}$ to bus 0, i.e., $R_{n_{k}^{0}n_{k}^{0}}$ and $X_{n_{k}^{0}n_{k}^{0}}$, which can be informed by CC), the second from all the other subtrees, and the third from the unclustered nodes. For convenience, denote by $\alpha_{i}=\sum_{j\in{\cal N}}R_{ij}(\overline{\mu}_{j}-\underline{\mu}_{j})=\alpha_{k,i}^{\text{in}}+\alpha_{k}^{\text{out}}$ with $\alpha_{k,i}^{\text{in}}$ denoting the first term in (15) and $\alpha_{k,i}^{\text{out}}$ the summation of the second and third terms in (15). Note that for $i\in{\cal N}_{k}$, $\alpha_{k,i}^{\text{in}}$ is accessible by RC $k$ cognizant of the topology of subtree $k$, and that $\alpha_{k}^{\text{out}}$ does not involve the network structure of any other subtrees but only that of the reduced network ${\cal T}^{r}$ known by CC.

For unclustered node $i\in{\cal N}_{0}$, similarly we can decompose $\sum_{j\in{\cal N}}R_{ij}(\overline{\mu}_{j}-\underline{\mu}_{j})$ as

[TABLE]

whose computation only requires the topology and line parameters of the reduced network ${\cal T}^{r}$ coordinated by CC.

Eqs. (15)–(16) motivate us to design a hierarchical distributed implementation of the primal-dual gradient algorithm (8), where CC and RCs are in charge of only a portion of the whole system: CC manages the reduced network ${\cal T}^{r}$ and coordinates RCs as well as unclustered nodes without knowing any structural or node-wise information within subtrees, and each RC $k$ manages its own subtrees ${\cal T}_{k}$ without knowing structural or node-wise information of the other subtrees or the reduced network. This also leads to a secure and privacy-preserving design such that regional topology and node-wise information are both protected.

We refer the detailed presentation of the single-phase hierarchical algorithm to [1] as its structure is similar to the multi-phase algorithm to be introduced later.

III-B Complexity Reduction

The hierarchical implementation not only enables parallel computation of the coupling terms in the gradient algorithm, but also largely reduces computational loads and communication overhead due to the following two reasons:

• •

The term $\sum_{h\in{\cal K},h\neq k}R_{n^{0}_{h}n^{0}_{k}}\sum_{j\in{\cal N}_{h}}(\overline{\mu}_{j}-\underline{\mu}_{j})$ requires less computation than the original $\sum_{j\in{\cal N}_{h},h\neq k}R_{ij}(\overline{\mu}_{j}-\underline{\mu}_{j})$;

• •

The intermediate computation result $\alpha_{k}^{\text{out}}$ is the same for all the nodes in ${\cal N}_{k}$, reducing a lot of repetitive computations.

III-B1 Computational Complexity Characterization

We compare the computational complexity in calculating the coupling term $R^{\top}(\overline{\bm{\mu}}-\underline{\bm{\mu}})$ by centrally coordinated algorithm and hierarchical distributed algorithm as follows.

• •

Centrally coordinated algorithm takes $N^{2}$ multiplications and $N(N-1)$ additions leading to total complexity of $\mathcal{O}(N^{2})$.

• •

For hierarchical distributed algorithm, computational complexity for general network and clustering strategy is challenging to characterize. We instead make simplification and provide intuitions on the computational reduction under the following ideal topological and clustering scenario: Assume that all $N$ nodes with voltage constraints are equally clustered into $K$ subtrees, each with $N_{k}=N/K$ nodes. This also ignores the third term in (15).

1. (a)

To calculate (15) for all $i\in{\cal N}_{k}$, each RC performs $N_{k}^{2}$ multiplications and $N_{k}(N_{k}-1)$ additions for the first term, $N_{k}-1$ additions for the second term, and $N_{k}$ additions to add the two terms. This results in $2N_{k}^{2}+N_{k}-1$ operations for each RC and $K(2N_{k}^{2}+N_{k}-1)$ for all RCs. 2. (b)

CC performs $K^{2}$ multiplications and $(K-1)K$ additions to put the second term of (15) together. 3. (c)

Adding up the results of last two steps and using $N=KN_{k}$ lead to a total of $2N^{2}/K+N+2K^{2}$ operations (ignoring lower-order terms in $K$). 4. (d)

Apparently the computation depends on the choice of $K$, e.g., when $K=1$ or $K=N$, the clustering is trivial and there is no computational reduction. When we choose $K=(N^{2}/2)^{1/3}$, the total computation burden is minimized to $\mathcal{O}(N^{4/3})$, largely decreased from $\mathcal{O}(N^{2})$.

Though in practice such ideal topological and clustering structure is uncommon, a simple geographical or administrative clustering strategy will nevertheless reduce a significant amount of computation. As will be shown in Section V, a 4-subtree clustering of the original network gives us a 4-fold acceleration in convergence speed by reducing computational loads. How to optimally cluster a distribution network remains an open question for our future work.

III-C Multi-Level Distributed Control

In large distribution systems, one subtree may still contain too many nodes to be handled efficiently. Meanwhile, some smaller areas within subtrees may want to shelter their topology and node-wise information from RC or CC, e.g., for security concerns.

This as well as the fractal structure of tree topology motivates us to consider deeper clustering and multi-level control, i.e., to apply similar approaches to cluster nodes within a subtree into smaller “sub-subtrees”, as illustrated in Fig. 3. Mathematically, this is done by decomposing $\alpha_{k,i}^{\text{in}}$ in the same as we do with $\alpha_{i}$ in (15). Even deeper clustering can be done likewise if necessary. We omit further details here.

IV Multi-Phase System Hierarchical Control

Given the intuitions built up from previous sections, we are now ready to design hierarchical distributed algorithms for multi-phase systems.

IV-A Multi-Phase System Modeling and OPF Formulation

Define $\mathfrak{i}:=\sqrt{-1}$. Let $a,b,c$ denote the three phases, and $\Phi_{i}$ the set of phase(s) of node $i\in{\cal N}$, e.g., $\Phi_{i}=\{a,b,c\}$ for a three-phase node $i$, and $\Phi_{j}=\{b\}$ for a single b-phase node $j$. Also, in a three-phase system one usually has $\Phi_{0}=\{a,b,c\}$. Define ${\cal N}^{\phi}\subseteq{\cal N}$ the subset of ${\cal N}$ collecting nodes that have phase $\phi$. Denote by $p_{i}^{\phi}$, $q_{i}^{\phi}$, $V_{i}^{\phi}$ and $v_{i}^{\phi}$ the real power injection, the reactive power injection, the complex voltage phasor, and the squared voltage magnitude, respectively, of node $i\in{\cal N}$ at phase $\phi\in\Phi_{i}$. Denote by $N_{\Xi}:=\sum_{i\in{\cal N}}|\Phi_{i}|=\sum_{\phi\in\Phi_{0}}|{\cal N}^{\phi}|$ the total cardinality of the multi-phase system, where $|\cdot|$ calculates the cardinality of a set. We make the following assumptions to obtain a linearized power flow model for the multi-phase system.

Assumption 2**.**

We consider a multi-phase distribution network where

1. 1)

the line losses are small and ignored, and 2. 2)

the magnitudes of three-phase voltages are approximately equal and the phase differences among three-phase voltages are close to $2\pi/3$, i.e., $\frac{V_{i}^{a}}{V_{i}^{b}}\approx\frac{V_{i}^{b}}{V_{i}^{c}}\approx\frac{V_{i}^{c}}{V_{i}^{a}}\approx e^{\mathfrak{i}2\pi/3}$.

Let $z_{\zeta\xi}$ be the phase impedance matrix of line $(\zeta,\xi)\in{\cal E}$. For example, if line $(\zeta,\xi)$ has three phases,

[TABLE]

where the off-diagonal elements represent mutual impedance between two phases. Similar to $R_{ij}$ and $X_{ij}$ defined in Eq. (3), we define

[TABLE]

the summarized impedance (if $\varphi=\phi$) or the mutual impedance (if $\varphi\neq\phi$) of the common path of node $i$ and $j$ leading back to node 0, and $\overline{Z}^{\varphi\phi}_{ij}$ its conjugate.

We denote by $\bm{v}_{{}_{\Xi}}=[[v_{1}^{\phi}]^{\top}_{\phi\in\Phi_{1}},\ldots,[v_{N}^{\phi}]^{\top}_{\phi\in\Phi_{N}}]^{\top}\in\mathbb{R}^{N_{\Xi}}$ the multi-phase squared voltage magnitude vector, and $\bm{p}_{{}_{\Xi}}=[[p_{1}^{\phi}]_{\phi\in\Phi_{1}}^{\top},\ldots,[p_{N}^{\phi}]^{\top}_{\phi\in\Phi_{N}}]^{\top}\in\mathbb{R}^{N_{\Xi}}$ and $\bm{q}_{{}_{\Xi}}=[[q_{1}^{\phi}]^{\top}_{\phi\in\Phi_{1}},\ldots,[q_{N}^{\phi}]^{\top}_{\phi\in\Phi_{N}}]^{\top}\in\mathbb{R}^{N_{\Xi}}$ the multi-phase power injection vectors. We then extend the linearization (2) to its multi-phase counterpart written as

[TABLE]

where $\tilde{\bm{v}}_{{}_{\Xi}}\in\mathbb{R}^{N_{\Xi}}$ is a constant vector depending on squared voltage magnitudes at all phases of the slack bus, and the voltage-to-power sensitivity matrices $R_{{}_{\Xi}},X_{{}_{\Xi}}\in\mathbb{R}^{N_{\Xi}\times N_{\Xi}}$ are determined by the linear approximation method developed for multi-phase system [17, 13]. The elements of $R_{{}_{\Xi}},X_{{}_{\Xi}}$ are calculated as follows (we use $a=0$, $b=1$, and $c=2$ when calculating $\phi-\varphi$):

[TABLE]

for any $\varphi\in\Phi_{i}$, $\phi\in\Phi_{j}$, $i,j\in{\cal N}$, with $\omega=e^{-\mathfrak{i}2\pi/3}$ and $\mathfrak{Re}\{\cdot\}$ and $\mathfrak{Im}\{\cdot\}$ denoting the real part and the imaginary part of a complex number. Note that when $\varphi=\phi$, Eqs. (19) coincide with $R_{ij}$ and $X_{ij}$ in Eqs. (3) for any nodes $i,j\in{\cal N}$; otherwise, Eqs. (19) calculate the summarized mutual impedance—rotated by phase difference $\pm 2\pi/3$—of the common path of nodes $i,j\in{\cal N}$ leading back to node 0.

We further extend Eq. (4) to its multi-phase counterpart as

[TABLE]

where $P_{I}^{\phi}$ is the total uncontrollable (inelastic) power injection at phase $\phi\in\Phi_{0}$. We use $\bm{v}_{{}_{\Xi}}(\bm{p}_{{}_{\Xi}},\bm{q}_{{}_{\Xi}})$ and $P_{0}(\bm{p}_{{}_{\Xi}})$ to represent Eqs. (18) and (20), and formulate the OPF problem for the multi-phase system as follows:

[TABLE]

Associate dual variables $\overline{\bm{\mu}}_{{}_{\Xi}}$ and $\underline{\bm{\mu}}_{{}_{\Xi}}$ with (21c) and we write the regularized Lagrangian of (21) as:

[TABLE]

Note that the multi-phase sensitivity matrices $R_{{}_{\Xi}}$ and $X_{{}_{\Xi}}$ in Eqs. (19) have the similar structure as their single-phase counterparts $R$ and $X$ defined in Eqs. (3), i.e., the values of $\partial_{p_{j}^{\varphi}}v_{i}^{\phi}$ and $\partial_{q_{j}^{\varphi}}v_{i}^{\phi}$ for any $i,j\in{\cal N}$ only depend on the common path of $i$ and $j$ leading back to node 0, as well as their angle difference $\phi-\varphi$. This motivates us to design a similar hierarchical control structure for the multi-phase system.

IV-B Multi-Phase Hierarchical Distributed Algorithm

The primal-dual gradient algorithm for solving the regularized Lagrangian of the convex optimization problem (21) reads:

[TABLE]

where (22a)–(22d) are for all $\phi\in\Phi_{i}$ and all $i\in{\cal N}$. One can obtain similar results as in Theorem 1–2 for the Lagrangian ${\cal{L}}^{\Phi}_{\eta}$ and the primal-dual gradient algorithm (22). We omit the details here to avoid repetition.

Note that the last terms in the primal update steps (22a)–(22b) not only couple all nodes but also multiple phases together. Nevertheless, similar to Eq. (15), we can decompose the coupling terms based on Eqs. (17)–(19) along with the radial topology of the network. To this end, we consider the same subtree structure and notations defined in Section III. We use real power updates for illustration.

For clustered node $i\in{\cal N}_{k}$, we have the following decomposition for any $\phi\in\Phi_{i}$:

[TABLE]

[TABLE]

with $Z^{\varphi\phi}_{n_{h}^{0}n_{k}^{0}}=\sum_{(\zeta,\xi)\in{\cal E}_{n_{k}^{0}}\cap{\cal E}_{n_{h}^{0}}}z^{\varphi\phi}_{\zeta\xi}$ summarizing the (mutual) impedance on the common path of subtree $h$ and subtree $k$ and $Z^{\varphi\phi}_{jn_{k}^{0}}=\sum_{(\zeta,\xi)\in{\cal E}_{j}\cap{\cal E}_{n_{k}^{0}}}z^{\varphi\phi}_{\zeta\xi}$ summarizing the (mutual) impedance on the common path of unclustered node $j$ and subtree $k$.

For unclustered node $i\in{\cal N}_{0},\forall\phi\in\Phi_{i}$, similarly we have:

[TABLE]

One can apply similar approaches to obtain the decomposed results to update reactive power injections. Then the resultant equivalent form of Eqs. (22) can be implemented in a hierarchical distributed way with the collaboration of RCs and CC. We summarize the resultant algorithm for multi-phase systems based on the primal-dual gradient algorithm (22) and the decomposition strategies (24)–(25) in Algorithm 1 where we use $s^{\phi}_{i}(t)=\alpha^{\phi}_{i}(t)+\mathfrak{i}\beta^{\phi}_{i}(t)$ for notational simplicity. We also use $s_{k}^{\phi_{\text{out}}}$ and $s^{\phi_{\text{in}}}_{k,i}$ to denote the parts of $s^{\phi}_{i}(t)$ that are calculated outside of subtree $k$ and within subtree $k$, respectively. See Fig. 4 for an illustration of clustering and information flows in Algorithm 1. Due to their mathematical equivalence, Algorithm 1 and Eqs. (22) share the same dynamic properties, which will also be illustrated in Section V.

IV-C Nonlinear Power Flow Analysis

In this part, we implement Algorithm 1 with nonlinear power flow, which can be either balanced or unbalanced, and characterize its performance. To this end, we denote the voltage magnitudes and the total active power loads updated by the nonlinear multi-phase power flow equations by $(\hat{\bm{v}}_{{}_{\Xi}},\hat{P}_{0})=\mathcal{F}(\bm{p}_{{}_{\Xi}},\bm{q}_{{}_{\Xi}})$, where the nonlinear power flow $\mathcal{F}$ can be updated based on some specific model from literature or can be measured from physical power flow.

Assumption 3**.**

The discrepancy between the linearized model and the nonlinear model is bounded, i.e., there exists some constants $e_{1},e_{2}>0$ such that given any $(\bm{p}_{{}_{\Xi}},\bm{q}_{{}_{\Xi}})\in\bm{{\cal Y}}$ we have

[TABLE]

Similar to Section II-C, we define $\bm{\mu}_{{}_{\Xi}}:=[\underline{\bm{\mu}}_{{}_{\Xi}}^{\top},\overline{\bm{\mu}}_{{}_{\Xi}}^{\top}]^{\top}$ and $\bm{z}_{{}_{\Xi}}:=[\bm{p}_{{}_{\Xi}}^{\top},\bm{q}_{{}_{\Xi}}^{\top},\bm{\mu}_{{}_{\Xi}}^{\top}]^{\top}$ to equivalently rewrite the primal-dual gradient dynamics (22) as

[TABLE]

with its gradient operators denoted by

[TABLE]

Note that Lipschitz continuity and strong monotonicity hold for $T_{{}_{\Xi}}(\bm{z}_{{}_{\Xi}})$ with positive constants $L_{{}_{\Xi}}$ and $M_{{}_{\Xi}}$ for any feasible $\bm{z}_{{}_{\Xi}},\bm{z}^{\prime}_{{}_{\Xi}}$ as follows:

[TABLE]

To implement the gradient algorithm with nonlinear power flow, we replace $C^{\prime}_{0}(P_{0}(t))$ in (22a) with $C^{\prime}_{0}(\hat{P}_{0}(t))$, $v^{\phi}_{i}(t)$ in (22c)–(22d) with $\hat{v}^{\phi}_{i}(t)$, and (22e)–(22f) with $(\hat{\bm{v}}_{{}_{\Xi}}(t+1),\hat{P}_{0}(t+1))=\mathcal{F}(\bm{p}_{{}_{\Xi}}(t+1),\bm{q}_{{}_{\Xi}}(t+1))$. In other words, the values of $\hat{\bm{v}}_{{}_{\Xi}},\hat{P}_{0}$ are updated by the nonlinear power flow while the gradient of $\hat{\bm{v}}_{{}_{\Xi}},\hat{P}_{0}$ with respect to decision variables are calculated based on the linearized model (18)–(20). Similar model-based feedback control has been applied and characterized with provable performance in recent literature, e.g., [34, 3, 7].

We denote the gradient operator with nonlinear power flow by $\hat{T}_{{}_{\Xi}}(\bm{z}_{{}_{\Xi}})$, and the resultant projected gradient algorithm as:

[TABLE]

Based on Assumption 3 and Lipschitz continuity of our gradient operator, we must have

[TABLE]

for some constant $\rho>0$. We characterize the convergence of dynamics (36) next.

Theorem 3**.**

Under Assumptions 1–3 and a stepsize $\epsilon$ chosen according to

[TABLE]

for some $\bar{\epsilon}$, dynamics (36) converges to the saddle point $\bm{z}_{{}_{\Xi}}^{*}$ of ${\cal{L}}^{\Phi}_{\eta}$ as

[TABLE]

Proof.

The distance between $\bm{z}_{{}_{\Xi}}(t+1)$ and $\bm{z}_{{}_{\Xi}}^{*}$ can be calculated as:

[TABLE]

[TABLE]

with $\Delta=1+\epsilon^{2}L_{{}_{\Xi}}^{2}-2\epsilon M_{{}_{\Xi}}$. In the above, the first inequality comes from non-expansiveness of projection operator, the second from (37), the third from (35a)–(35b), and the last inequality is obtained by recursively executing previous steps. Then based on (35c) and condition (38) one has $0<\Delta\leq\bar{\Delta}<1$ for some $\bar{\Delta}$ such that when $t\rightarrow\infty$, (39) follows.

Theorem 3 indicates that, despite the discrepancy introduced by linear and nonlinear power flows, dynamics (36) converges to the saddle point of ${\cal{L}}^{\Phi}_{\eta}$ with bounded distance $\frac{\rho}{2M_{{}_{\Xi}}/\epsilon-L_{{}_{\Xi}}^{2}}$. Therefore, given any accuracy requirement $\delta>0$, if stepsize is chosen such that $\epsilon<\min\{\frac{2M_{{}_{\Xi}}\delta}{\rho+L^{2}_{{}_{\Xi}}\delta},\frac{2M_{{}_{\Xi}}}{L^{2}_{{}_{\Xi}}}\}$, we have

[TABLE]

V Numerical Results

A three-phase unbalanced, 11,000-node test feeder is constructed by connecting an IEEE 8,500-node test feeder and a modified EPRI Ckt7 test feeder at the substation. Fig. 4 shows the single-line diagram of the feeder with its line width proportional to the nominal power flow on it. The primary side of the feeder is modeled in detail, while the loads on the secondary side are lumped into corresponding distribution transformers, resulting in a 4,521-node network with 1,043 controllable (aggregated) loads. We group all the nodes into unclustered nodes and four subtrees marked in Fig. 4. Subtrees 1–4 contain 357, 222, 310, and 154 nodes with controllable loads, respectively. We fix the loads on all 292 unclustered nodes for simplicity.

The three-phase unbalanced nonlinear power flow model is simulated in OpenDSS. With default control of capacitors and regulators in OpenDSS [35], one achieves the voltage profile shown in Fig. 6 with orange dots, where under-voltages are observed. We next disable the control of all the capacitors and regulators to obtain the heavily under-voltage scenario marked with blue dots in Fig. 6. We implement Algorithm 1 with this scenario as the initial condition.

The simulation is conducted on a laptop with Intel Core i7-7600U CPU @ 2.80GHz 2.90GHz, 8.00GB RAM, running Python 3.6 on Windows 10 Enterprise Version.

V-A Numerical Performance Evaluation

For each controllable node $i$ of phase $\phi$, we consider minimizing the cost of its deviation from its nominal (most preferred) load level $(p_{i}^{\phi}(0),q_{i}^{\phi}(0))$, i.e., $C^{\phi}_{i}(p^{\phi}_{i},q^{\phi}_{i})=(p^{\phi}_{i}-p_{i}^{\phi}(0))^{2}+(q^{\phi}_{i}-q_{i}^{\phi}(0))^{2}$. We focus on voltage regulation here. So, we set $C_{0}(P_{0})$ to $0.0005(P_{0}-\tilde{P}_{0})^{2}$ with a small weight and $\tilde{P}_{0}=P_{0}(0)$. $\underline{v}_{i}$ and $\overline{v}_{i}$ are uniformly set to $0.95$ p.u. and $1.05$ p.u., respectively. We implement Algorithm 1 with a constant stepsize $3.5\times 10^{-4}$ for the primal update, and $3.5\times 10^{-3}$ for the dual.

V-A1 Convergence and Computational Speed Improvement

It takes about 1730 iterations to reach 1% of the optimal value and 3000 iterations to reach the optimal; see the red curve in Fig. 6 (right). This result is also comparable to those in literature, e.g., [9].

The convergence dynamics of the centrally coordinated algorithm is identical to the hierarchical distributed implementation as sampled and illustrated in Fig. 5. However, due to the overall complexity reduction, the computational time per iteration is reduced from more than 4s to 1s. Meanwhile, if we consider parallel computation111We imitate the performance of parallel computation on one PC by timing the computation for CC and each RC. If the algorithm is to be implemented on multiple devices, we also need consider communication delay which is ignored here. and takes the slowest cluster to estimate the overall time consumed, each iteration only takes 0.37s. Therefore, the overall computational speed improvement is more than 10 folds without compromising any accuracy of the OPF solution. Noticing that the computational time at each cluster is approximately proportional to the square of node number, we expect faster performance if more clusters, each with a smaller number of nodes, are divided.

V-A2 Voltage Regulation

We plot the regulated voltages obtained by Algorithm 1 with green dots in Fig. 6. Note that the voltage magnitudes of all nodes are strictly constrained within the $[0.95,1.05]$ p.u. bound. In contrast, the default control of regulators and capacitors cannot guarantee that all the voltages are within this bound.

V-A3 Comparison with Single-Phase Algorithm

We apply the single-phase algorithm with the same setup for comparison. The detailed results of the single-phase algorithm are referred to [1]. The results show that it takes multi-phase algorithm more time to execute each iteration (1s v.s 0.45s), which is expected since the multi-phase algorithm executes more computation by considering inter-phase sensitivity. On the other hand, thanks to more accurate linearization, the minimal cost obtained by the multi-phase algorithm is 17971, 4.6% smaller than the value of 18840 obtained by the single-phase algorithm in [1], as illustrated in Fig. 7. This result also echoes with Eq. (39), which indicates that accuracy of linearization model $\rho$ affects algorithm performance.

VI Conclusion

We proposed a hierarchical distributed implementation of the primal-dual gradient algorithm to solve an OPF problem. The objective of OPF is to minimize the total cost over all the controllable DERs and a cost associated with the total network load, subject to voltage regulation constraints. By utilizing the information structure of tree/subtrees to reduce and distribute computational loads the proposed implementation is scalable to large multi-phase distribution networks. Performance of our design is analytically characterized and numerically corroborated. The significant improvement in convergence speed shows the great potential of the proposed method for grid optimization and control in real time.

Acknowledgments

This work was authored in part by the National Renewable Energy Laboratory, operated by Alliance for Sustainable Energy, LLC, for the U.S. Department of Energy (DOE) under Contract No. DE-EE-0007998. Funding provided by U.S. Department of Energy Office of Energy Efficiency and Renewable Energy Solar Energy Technologies Office. The views expressed in the article do not necessarily represent the views of the DOE or the U.S. Government. The U.S. Government retains and the publisher, by accepting the article for publication, acknowledges that the U.S. Government retains a nonexclusive, paid-up, irrevocable, worldwide license to publish or reproduce the published form of this work, or allow others to do so, for U.S. Government purposes.

-A Proof of Lemma 1

Proof.

Denote by $f(\bm{y})=\sum_{i\in{\cal N}}C_{i}(p_{i},q_{i})+C_{0}(P_{0}(\bm{p}))$ and $\bm{\mu}^{\top}\bm{g}(\bm{y})=\underline{\bm{\mu}}^{\top}(\underline{\bm{v}}-\bm{v}(\bm{p},\bm{q}))+\overline{\bm{\mu}}^{\top}(\bm{v}(\bm{p},\bm{q})-\overline{\bm{v}})$ for simplicity. Then $T(\bm{z})$ can be equivalently decomposed into the following operators:

[TABLE]

We can verify that the first operator $\begin{bmatrix}\nabla_{\bm{y}}f(\bm{y})\\ \nabla_{\bm{\mu}}\frac{\phi}{2}\|\bm{\mu}\|_{2}^{2}\end{bmatrix}$ is strongly monotone since $f(\bm{y})$ and $\frac{\phi}{2}\|\bm{\mu}\|_{2}^{2}$ are strongly convex in $\bm{y}$ and $\bm{\mu}$, respectively. The second (linear) operator is monotone since

[TABLE]

Therefore, $T(\bm{z})$ is a strongly monotone operator as the result of combining a strongly monotone operator and a monotone operator.

-B Proof of Theorem 2

Proof.

We have

[TABLE]

where the first inequality comes from non-expansiveness of projection operator, and the second from (10) and (11). Then by Lemma 2 and condition (13), we have $0\leq 1+\epsilon^{2}L^{2}-2\epsilon M<1$, i.e., $\|\bm{z}(t+1)-\bm{z}^{*}\|_{2}^{2}\leq\Delta\|\bm{z}(t)-\bm{z}^{*}\|_{2}^{2}$ for some constant $0<\Delta<1$.

Therefore, dynamics (9) converges to the unique saddle point of (7) exponentially fast.