Hierarchical Optimal Power Flow with Improved Gradient Evaluation

We model a single-phase equivalent power distribution network as a directed tree graph $\mathcal{T}:=\{\mathcal{N}^{+},\mathcal{E}\}$, where $\mathcal{N}^{+}=\mathcal{N}\cup\{0\}$, with $0$ indexing the root node (the substation, also known as the slack bus), and $\mathcal{N}=\{1,...,N\}$ indexing other nodes. The set $\mathcal{E}$ collects the ordered pairs of nodes representing the lines in the network. We define the directions of lines as pointing away from the root, e.g., in line $(i,j)$, node $i$ is closer to the root than node $j$. Let $p_{i}$ and $q_{i}$ denote the net active and reactive power injections (i.e., power supply minus consumption) and $v_{i}$ denote the squared voltage magnitude at each node $i\in\mathcal{N}^{+}$. Let $\ell_{ij}$ denote the squared current magnitude, $z_{ij}=r_{ij}+\boldsymbol{i}x_{ij}$ denote the series impedance, and $P_{ij}$ and $Q_{ij}$ denote the sending-end active and reactive power on each line $(i,j)\in\mathcal{E}$.

We adopt the distribution power flow model [12, 13]:

Assume the power injections $(p_{i},q_{i})$ to be controllable within a given operating region:

Define vectors $\boldsymbol{p}:=\left[p_{1},\ldots,p_{N}\right]^{\top}$, $\boldsymbol{q}:=\left[q_{1},\ldots,q_{N}\right]^{\top}$, $\boldsymbol{v}:=\left[v_{1},\ldots,v_{N}\right]^{\top}\in\mathbb{R}^{N}$. Following [8], we express the voltage vector as a function of power injections, i.e., $\boldsymbol{v}(\boldsymbol{p},\boldsymbol{q})$, implicitly defined by (1). Consider the following OPF problem:

where $f_{i}$ is a strongly convex cost function for the control of power injections at node $i\in\mathcal{N}$, inequality (3b) is element-wise, given constant squared voltage limits $\underline{\boldsymbol{v}}$ and $\overline{\boldsymbol{v}}$.

Let $\underline{\boldsymbol{\mu}}$ and $\overline{\boldsymbol{\mu}}$ be the dual variables associated with (3b) to penalize the violation of voltage lower and upper bounds, respectively. To design a convergent primal-dual algorithm, we consider the regularized Lagrangian of (3):

where $\boldsymbol{\mu}=\left(\underline{\boldsymbol{\mu}},\overline{\boldsymbol{\mu}}\right)$ and $\epsilon>0$ is a constant factor for the regularization term. Function (II-B) is strongly convex in primal variables $(\boldsymbol{p},\boldsymbol{q})$ and strongly concave in dual variables $\boldsymbol{\mu}$. The saddle point of (II-B) serves as an approximate solution to (3), with a bounded error that is linear in $\epsilon$ [10, 14]. The following primal-dual gradient algorithm has been commonly used to iteratively approach such a saddle point:

where $\boldsymbol{u}:=(\boldsymbol{p},\boldsymbol{q})$ collects the controllable power injections; $f(\boldsymbol{u}):=\sum_{i\in\mathcal{N}}f_{i}(p_{i},q_{i})$ is the objective (3a); $\boldsymbol{v}(t)=\boldsymbol{v}\left(\boldsymbol{u}(t)\right)$ both denote the voltage profile under power injections $\boldsymbol{u}(t)$ at time step $t$. Note that $\boldsymbol{v}\left(\boldsymbol{u}(t)\right)$ does not have to be calculated by (1) or any other mathematical model of power flow; instead, it can be measured from the power network after actuating the network with the intermediate decisions $\boldsymbol{u}(t)$ during the process (5). This feedback-based implementation has been widely adopted to compensate for likely modeling errors [7, 8, 9, 11]. The subscripts $\mathcal{Y}$ and $+$ represent the projections onto $\prod_{i\in\mathcal{N}}\mathcal{Y}_{i}$ and the nonnegative orthant, respectively. Without loss of generality, we use constant step sizes $\sigma_{u},\sigma_{\mu}>0$ to update the primal and dual variables, respectively.

The gradient $\partial{f}/\partial{\boldsymbol{u}}$ in (5a) can be easily computed, and the voltage $\boldsymbol{v}(t)$ in (5b)–(5c) can be conveniently measured, both locally at each node $i\in\mathcal{N}$. Therefore, the remaining key challenge is to compute the gradient $\partial{\boldsymbol{v}}/\partial{\boldsymbol{u}}$ that potentially couples all the nodes in the power network. This will be addressed in the next section.

In power flow equations (1), the nonlinear and implicit dependence of voltage $\boldsymbol{v}$ on power injections $\boldsymbol{u}$ makes it difficult to quickly and precisely compute the gradient $\partial{\boldsymbol{v}}/\partial{\boldsymbol{u}}$. Reference [8] proposed a backward-forward sweep method for gradient calculation, which is computationally inefficient in large networks. In contrast, prior work [10] considered the linearized power flow model in [12, 13]:

which yields the following approximation of $\partial{\boldsymbol{v}}/\partial{\boldsymbol{u}}$:

where $\mathbb{P}_{j\wedge h}$ denotes the common part of the unique paths from nodes $j$ and $h$ back to the root node.

We illustrate the error of linear model (6) by introducing the following Lemma.

Consider an arbitrary node $h\in\mathcal{N}$, and $u_{h}:=(p_{h},q_{h})$. We take the partial derivatives of the variables in the nonlinear power flow model (1):

The exact partial derivatives are hard to solve from (9) due to their complex interdependence. To facilitate computation, we first simplify (9d) by replacing the partial derivatives with their approximations under the linear model (6):

We recall the following results from [8] for all the nodes $h\in\mathcal{N}$ and lines $(i,j)\in\mathcal{E}$:

where $\mathds{1}(j\in\mathbb{P}_{h})$ is an indicator that equals 1 if node $j$ lies on the unique path from node $h$ to the root, and 0 otherwise. By (11) and (7), we can calculate $\partial\hat{\boldsymbol{\ell}}/\partial{\boldsymbol{u}}$ in (10) as follows:

The results in (7), (11), (12) can replace the partial derivatives on the right-hand side (RHS) of (9c) to obtain an improved gradient evaluation $\partial{\boldsymbol{v}}/\partial{\boldsymbol{u}}$ as:

where node $i$ is the unique parent of node $j$ in the directed tree network. By (7a), there is $R_{jh}=R_{ih}+2r_{ij}\mathds{1}(j\in\mathbb{P}_{h})$, which converts the result in (13a) into:

The first term on the RHS of (14) is the same as (7a) derived from the linear model (6), while the second and third terms compensate for the effects of the quadratic terms that were neglected from (6). The improved partial derivatives (13b) with respect to reactive power injections $q_{h}$ have the same structure as (13a), and thus we will not repeat.

We consider the single-phase equivalent models of the IEEE 37-node and 123-node test networks. They are clustered into subtrees as shown by Figure 1 (left panel) and Figure 2, respectively.

For the ease of illustration, we make the following modifications to the original network models on the IEEE PES website (https://cmte.ieee.org/pes-testfeeders/resources/):

1. 1.

   We model each multi-phase line as a single-phase line with the average impedance of multiple phases. We also convert the multi-phase wye- and delta-connected loads into single-phase loads by taking their average over multiple phases.

2. 2.

   We multiply the original load data of the 37-node network by six, and those of the 123-node network by two, to create scenarios with serious voltage issues.

3. 3.

   All the loads are treated as constant-power loads. Detailed models of capacitors, regulators, and breakers are not simulated.

For each network, let $\mathcal{N}_{L}$ denote the set of nodes that host nonzero loads. We denote the negative net power injections from the IEEE load data by $(\underline{p}_{h},\underline{q}_{h})$ for all load nodes $h\in\mathcal{N}_{L}$, and treat them as the nominal injections or the most preferred injections. The feasible power injection regions in (2) are defined as $\mathcal{Y}_{h}=\left\{(p_{h},q_{h})\mid\underline{p}_{h}\leq p_{h}\leq 0.3\underline{p}_{h}<0,\underline{q}_{h}\leq q_{h}\leq 0.3\underline{q}_{h}<0\right\}$. The OPF objective function in (3a) is defined as $f_{h}(p_{h},q_{h})=(p_{h}-\overline{p}_{h})^{2}+(q_{h}-\overline{q}_{h})^{2}$ for all $h\in\mathcal{N}_{L}$, which aims to minimize the disutility caused by the deviation from the nominal loads.

For each network, the voltage magnitude of the root (slack) node is fixed at $1.05$ per unit (p.u.). The lower and upper limits for safe voltage are set at $0.95$ p.u. and $1.05$ p.u., respectively. The step sizes for primal and dual variable updates are empirically chosen as $\sigma_{u}=2\times 10^{-3}$ and $\sigma_{\mu}=1\times 10^{-3}$, respectively.

For each given power injection $\boldsymbol{u}=\left(\boldsymbol{p},\boldsymbol{q}\right)$, an OpenDSS power flow simulation is performed to obtain the corresponding voltage $\boldsymbol{v}(\boldsymbol{u})$ and the associated line-flow quantities $(\boldsymbol{P},\boldsymbol{Q},\boldsymbol{\ell})$, which are treated as the feedback signals measured from the power network. The Python 3.7 programs for OPF algorithms and the OpenDSS simulations are run on a laptop equipped with Intel Core i7-9750H CPU @ 2.6 GHz, 16 GB RAM, and Windows 10 Professional OS.

Due to space limitation, we only show the simulation results in the 123-node network model. The results in the 37-node network are similar.

Voltage safety: The voltage magnitudes at different nodes of the 123-node network are plotted in Figure 3, for three cases. The first case, referred to as “no control”, just takes the nominal power injections $(\underline{\boldsymbol{p}},\underline{\boldsymbol{q}})$ from the IEEE data, which causes severe violation of the voltage lower bound. The second case, referred to as “linear control”, applies the primal-dual gradient algorithm in [10] based on the linearized power flow model (6) to determine power injections, which enhances the overall voltage level but still leaves most nodes below the voltage lower bound. The third case, referred to as “improved control (ours)”, solves the OPF with the proposed Algorithm 1 based on our improved gradient evaluation method, which effectively lifts all the voltages above the lower bound. This result verifies our analysis in Section III that the improved gradient evaluation can prevent voltage violation caused by the linear model that optimistically estimates the voltages to be safe.

Computational efficiency: Figure 4 shows the change of voltage at a node of the 123-node network (similar trends are displayed for the voltages at other nodes). This result verifies again that the voltage converges to a safe level under the proposed Algorithm 1, compared to the previous algorithm [10] that leads to a violation of voltage lower bound. Moreover, Algorithm 1 converges in fewer iterations. Indeed, Algorithm 1 spends 205 seconds to complete 2,000 iterations, compared to the previous algorithm in [10] which spends 166 seconds. This would be an acceptable increase, considering the improved voltage safety achieved by Algorithm 1.