A Theory of Solvability for Lossless Power Flow Equations -- Part I: Fixed-Point Power Flow

TL;DR

This paper introduces a new fixed-point formulation of lossless power flow equations, providing explicit approximations and convergence guarantees, advancing understanding of power flow solvability in both meshed and radial networks.

Contribution

It presents a novel fixed-point power flow model parameterized by graph-theoretic matrices, with proven convergence and solvability conditions for lossless power networks.

Findings

Iterates rapidly converge to high-voltage solutions

Explicit approximation accurately predicts power flow near base case

Derived conditions guarantee solution existence and uniqueness

Abstract

This two-part paper details a theory of solvability for the power flow equations in lossless power networks. In Part I, we derive a new formulation of the lossless power flow equations, which we term the fixed-point power flow. The model is stated for both meshed and radial networks, and is parameterized by several graph-theoretic matrices -- the power network stiffness matrices -- which quantify the internal coupling strength of the network. The model leads immediately to an explicit approximation of the high-voltage power flow solution. For standard test cases, we find that iterates of the fixed-point power flow converge rapidly to the high-voltage power flow solution, with the approximate solution yielding accurate predictions near base case loading. In Part II, we leverage the fixed-point power flow to study power flow solvability, and for radial networks we derive conditions guaranteeing the existence and uniqueness of a high-voltage power flow solution. These conditions (i) imply exponential convergence of the fixed-point power flow iteration, and (ii) properly generalize the textbook two-bus system results.

Tables

Table 1. TABLE I : Testing of Fixed-Point Power Flow

	Base Load			High Load
Test Case	FPPF	${\delta }_{\max }$	${\delta }_{\mathrm{avg}}$	FPPF	${\delta }_{\max }$
	Iters.	(p.u.)	(p.u.)	Iters.	(p.u.)
14 bus system	4	0.001	0.000	8	0.090
RTS 24	4	0.003	0.001	8	0.081
30 bus system	4	0.003	0.002	8	0.104
New England 39	4	0.006	0.004	8	0.086
57 bus system	5	0.011	0.003	8	0.118
RTS ’96 (3 area)	4	0.003	0.001	8	0.084
118 bus system	3	0.001	0.000	7	0.054
300 bus system	6	0.022	0.004	8	0.059
PEGASE 1,354	5	0.011	0.001	8	0.070
Polish 2,383 wp	4	0.003	0.000	8	0.078
PEGASE 2,869	5	0.015	0.002	8	0.098
PEGASE 9,241	6	0.063	0.003	9	0.133

Table 2. TABLE II : Solver success rates under random initializations for the IEEE 118 bus system.

IC Spread ($\alpha$)	NR	FDLF	FPPF
0.05	0.98	0.98	1.00
0.10	0.53	0.53	1.00
0.15	0.18	0.18	1.00
0.2	0.03	0.03	1.00
0.3	0.00	0.00	1.00
0.5	0.00	0.00	1.00
0.7	0.00	0.00	0.99
0.9	0.00	0.00	0.99

Equations

$$\Delta +4{\Gamma }_{g\ell }^2$$

$$\mathcal{E}={\mathcal{E}}^{\ell \ell }\cup {\mathcal{E}}^{g\ell }\cup {\mathcal{E}}^{gg},$$

$$A=\left(\begin{array}[]{c}A_{L}\\[2.15277pt] \hline\cr\\[-8.61108pt] A_{G}\end{array}\right)=\left(\begin{array}[]{c|c|c}A_{L}^{\ell\ell}&A_{L}^{g\ell}&\mathbbold{0}\\[2.15277pt] \hline\cr\\[-8.61108pt] \mathbbold{0}&A_{G}^{g\ell}&A_{G}^{gg}\end{array}\right)\,.$$

$$V=\left(\begin{array}{c}{V}_L\\ V_G\end{array}\right),Q=\left(\begin{array}{c}{Q}_L\\ Q_G\end{array}\right),B=\left(\begin{array}{cc}{B}_{LL}& B_{LG}\\ B_{GL}& B_{GG}\end{array}\right).$$

$$P_i$$

$$Q_i$$

$$V_L^{\ast }\triangleq -B_{LL}^{-1}{B}_{LG}{V}_G.$$

$$\mathsf{D}\triangleq [V_i^{\ast }{V}_j^{\ast }{B}_{ij}{]}_{(i,j)\in \mathcal{E}}.$$

$$\mathsf{L}\triangleq A\mathsf{D}{A}^{\mathsf{T}}.$$

$$\mathsf{S}\triangleq \frac{1}{4}[V_L^{\ast }]B_{LL}[V_L^{\ast }].$$

$$\mathsf{D}=\mathrm{blkdiag}({\mathsf{D}}_{\ell \ell },{\mathsf{D}}_{g\ell },{\mathsf{D}}_{gg}).$$

$$A$$

$$-\left(\begin{array}{ccc}{A}_L^{\ell \ell }(-)& A_L^{g\ell }(-)& \text{mathbbold}0\\ \text{mathbbold}0& \text{mathbbold}0& A_G^{gg}(-)\end{array}\right).$$

$$A(+{)}^{\mathsf{T}}V=A(+{)}^{\mathsf{T}}\left(\begin{array}{c}{V}_L\\ V_G\end{array}\right)$$

$$A(-{)}^{\mathsf{T}}V=A(-{)}^{\mathsf{T}}\left(\begin{array}{c}{V}_L\\ V_G\end{array}\right)$$

$$[V_i{V}_j{B}_{ij}{]}_{(i,j)\in \mathcal{E}}$$

$$=\left[\left(\begin{array}{c}{A}_L^{\ell \ell }(+{)}^{\mathsf{T}}{V}_L\\ A_G^{g\ell }(+{)}^{\mathsf{T}}{V}_G\\ A_G^{gg}(+{)}^{\mathsf{T}}{V}_G\end{array}\right)\right]\left[\left(\begin{array}{c}{A}_L^{\ell \ell }(-{)}^{\mathsf{T}}{V}_L\\ A_L^{g\ell }(-{)}^{\mathsf{T}}{V}_L\\ A_G^{gg}(-{)}^{\mathsf{T}}{V}_G\end{array}\right)\right][B_{ij}{]}_{(i,j)\in \mathcal{E}}.$$

$$\mathsf{D}=\left[\left(\begin{array}{c}{A}_L^{\ell \ell }(+{)}^{\mathsf{T}}{V}_L^{\ast }\\ A_G^{g\ell }(+{)}^{\mathsf{T}}{V}_G\\ A_G^{gg}(+{)}^{\mathsf{T}}{V}_G\end{array}\right)\right]\left[\left(\begin{array}{c}{A}_L^{\ell \ell }(-{)}^{\mathsf{T}}{V}_L^{\ast }\\ A_L^{g\ell }(-{)}^{\mathsf{T}}{V}_L^{\ast }\\ A_G^{gg}(-{)}^{\mathsf{T}}{V}_G\end{array}\right)\right][B_{ij}],$$

$$V_L=[V_L^{\ast }]v⟺v=[V_L^{\ast }{]}^{-1}{V}_L.$$

$$[V_i{V}_j{B}_{ij}{]}_{(i,j)\in \mathcal{E}}=\mathsf{D}\cdot [h(v)],$$

$$h(v)=\left(\begin{array}{c}{h}_{\ell \ell }(v)\\ h_{g\ell }(v)\\ h_{gg}(v)\end{array}\right)$$

$$h_e(v)=\{\begin{array}{ll}{v}_i{v}_j& \text{if }e=(i,j)\in {\mathcal{E}}^{\ell \ell }\\ v_j& \text{if }e=(i,j)\in {\mathcal{E}}^{g\ell }\\ 1& \text{if }e=(i,j)\in {\mathcal{E}}^{gg}\end{array}.$$

$$P=A[V_i{V}_j{B}_{ij}{]}_{(i,j)\in \mathcal{E}}\sin (A^{\mathsf{T}}\theta ),$$

$$P=A\mathsf{D}\cdot [h(v)]\cdot \sin (A^{\mathsf{T}}\theta ).$$

$$\psi (v,y)\triangleq [h(v){]}^{-1}\left(A^{\mathsf{T}}{\mathsf{L}}^{†}P+{\mathsf{D}}^{-1}Cy\right).$$

$$\mathbbold{0}_{c}=C^{\sf T}\boldsymbol{\arcsin}(\psi(v,y))\qquad{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}(\boldsymbol{\mathrm{mod}}\,2\pi)}\,.$$

$$\sin (\eta )=\psi (v,y).$$

$$P$$

$$\tilde{\psi }$$

$$P$$

$$=\underset{=\mathsf{L}}{\underbrace{A\mathsf{D}{A}^{\mathsf{T}}}}{\mathsf{L}}^{†}P+A\mathsf{D}{\mathsf{D}}^{-1}Cy=\Pi P+\underset{=\text{mathbbold}0}{\underbrace{AC}}y=P,$$

$$\sin (A^{\mathsf{T}}\theta )=\psi (v,y).$$

$$C^{\mathsf{T}}\arcsin (\psi (v,y))=2\pi C^{\mathsf{T}}k.$$

Full text

I Introduction

The power flow equations describe the balance and flow of power in a synchronous AC power system. The solutions of these equations (also called operating points of the network) describe the configurations of voltages and currents which (i) are physically consistent with Kichhoff’s and Ohm’s laws, and (ii) meet the prescribed boundary conditions, specified in terms of fixed power injections or fixed voltage levels at particular nodes in the network. Knowledge of the current system operating point is crucial, as is understanding how the current operating point will change as control actions are taken or as unexpected contingencies occur. As such, the power flow equations are embedded in nearly every power system analysis or control problem, including optimal power flow and its security-constrained variants, transient and voltage stability assessment, contingency screening, short-circuit analysis, and wide-area monitoring/control [1].

As the equations are nonlinear, the existence of real-valued solutions is not guaranteed: lightly loaded networks typically possess many solutions [2], while a network which is loaded sufficiently will possess none. Despite this potential for both multiple reasonable solutions and infeasibility, typically there is a single desirable solution, characterized by high voltage magnitudes at buses and small inter-bus current flows. This solution is often termed stable, as it behaves in a manner consistent with the intuition of operators, and moreover, is a locally exponentially stable equilibrium point for some simplified dynamic grid models [3, 4]. The ability to accurately and consistently calculate this high-voltage solution is incredibly important, and fairly reliable numerical techniques are available for this purpose [2, 5, 6]. While our results have computational implications, our main interest and motivation is the question of power flow feasibility/solvability: for what classes of networks and loading scenarios can we guarantee that the power flow equations are solvable for a useful solution, and what can be rigorously said about this solution?

Aside from intellectual merit, there are at least two important engineering motivations for understanding solvability. The first is to better understand the convergence of iterative numerical algorithms for solving power flow equations. When a power flow solver diverges, it may be because of a numerical instability in the algorithm, an initialization issue [7, 8], or it may be because no power flow solution exists to be found [9]. Without a coherent theory of power flow solvability, it is difficult to distinguish between these cases. Our proposed algorithm in Part I is based on a carefully chosen fixed-point iteration. For some restricted classes of networks, our theoretical results in Part II provide a certificate that a unique power flow solution exists, and specify a large set of initializations from which our fixed-point iteration converges exponentially to this solution.

The second motivation comes from the desire to operate power systems safely yet non-conservatively. Due to the large capital costs of transmission infrastructure investment, system operators are incentivized to operate power networks close to their maximum power transfer limits. The present work is an additional step towards characterizing these nonlinear transfer limits, and understanding in a precise mathematical way how the transfer limits depend on the internal structure of the grid. In this context, our results in Part II provide a topology-dependent loading margin for the grid. This loading margin can serve as a solvability certificate, or as a lower bound on the distance to the maximum power transfer boundary.

I-A Contributions of Part I and Preview of Part II Results

This two-part paper presents a new model of power flow in lossless networks, and then leverages this model to obtain (i) a new iterative power flow algorithm, (ii) an approximation of the high-voltage solution, and (iii) new theoretical results on power flow solvability. Our new model is inspired by the way that phase angles are eliminated in the standard textbook analysis of the two-bus PV ​​-PQ power system [10, Chapter 2]. We begin with the lossless power flow equations in polar form, with voltage variables $(V,\theta)$ and power variables $(P,Q)$, and proceed to eliminate the phase angles from the model. The state variables of our new model are (i) the normalized voltages $v_{i}=V_{i}/V_{i}^{*}$ at PQ buses, where $V_{i}^{*}$ denotes the open-circuit voltage at the $i$th PQ bus, and (ii) a set of slack variables which enforce Kichhoff’s voltage law around cycles in meshed networks. Voltage phase angles are uniquely recovered as functions of these state variables.111In that phase angles are absent, our model is conceptually similar to the Baran-Wu branch flow model [11, 5].

For networks without cycles, these slack variables are discarded, and the model can quickly be manipulated into the fixed-point form $v=f(v)$, where the function $f$ depends on the grid topology, parameters and loading. Motivated by this important radial case, we call our reformulation the fixed-point power flow (FPPF). The FPPF model is the main result of Part I, and is summarized in Theorem III.5. In Section III-G we show how the FPPF model naturally leads to an explicit approximate solution to the power flow equations, yielding voltage magnitudes at PQ buses and phase angles at all buses as explicit functions of active and reactive power injections.

In Section IV we numerically study the FPPF and the accompanying approximate solution using standard IEEE test cases. We show that the lossless power flow equations can be quickly and reliably solved by iterating the FPPF, for both lightly and heavily loaded systems. The convergence is exponential, at a rate comparable to the fast-decoupled power flow approach, and is extremely insensitive to the choice of initialization. We also show that our approximate solution is quite accurate in these same standard test cases.

Throughout Part I and Part II we restrict our attention to lossless networks, for two main reasons. Firstly, many high-voltage transmission networks are approximately lossless, with resistance/reactance ratios below $0.2$ (see Section IV). For such networks, practice has shown that the lossless assumption is not restrictive. Indeed, the standard power flow model used for dispatch — the linearized “DC Power Flow” — explicitly relies on this lossless assumption [12], and is widely used in industry. Secondly, power flow solvability remains poorly understood for transmission networks, and the lossless case should be understood before attempting an analysis of the lossy case. We comment further on resistances in Section IV and in our conclusions in Part II.

As an informal preview of our main result in Part II, the pair of existence and uniqueness conditions we derive for radial networks with no connections between PQ buses are

[TABLE]

The first inequality captures voltage stability of PQ load buses: $\Delta\in{[0,1)}$ is related to reactive power loading, while $\Gamma_{g\ell}\in{[0,\frac{1}{2})}$ is related to active power flow between generators and loads. Roughly speaking, this inequality ensures that voltage magnitudes at PQ buses stay high. The second inequality on $\Gamma_{gg}$ is an angle stability condition between generator PV buses, and ensures that phase angle differences between PV buses stay relatively small. The quantities $\Delta,\Gamma_{g\ell}$ and $\Gamma_{gg}$ depend only on the data of the power flow problem. These conditions are an exact generalization of conditions found in standard textbooks [10, Chapter 2], [1, Section 8.1.1] for the canonical two-bus PV ​​-PQ power system, generalizing the so-called circle diagram [1, Figure 8.3]. They also unify and extend recent solvability conditions developed for decoupled active power flow [4] and for decoupled reactive power flow [13], which in the above notation read as $\Gamma_{gg}<1$ and $\Delta<1$, respectively. We also present weaker results for networks with connections between PQ buses, which guarantee only the existence of a solution.

I-B Structure of Paper

Section II formally states our modeling assumptions leading to the standard model of coupled, lossless power flow used in the remainder of both papers.

Section III contains the main results of Part I. We introduce the stiffness matrices (Section III-A), derive the fixed-point power flow model in a step-by-step fashion (Sections III-B–III-E), discuss the derivation (Section III-F), and derive an approximate power flow solution based on the FPPF (Section III-G. The derivation is presented for meshed networks, with the result for radial networks stated as a corollary.

Section IV validates our results numerically on standard test cases, while Section V summarizes and concludes. The remainder of this section summarizes some vector and matrix notation used extensively throughout the paper, some of which is non-standard but convenient.

I-C Preliminaries and Notation

Sets, vectors and functions: For a finite set $\mathcal{N}$, $|\mathcal{N}|$ is its cardinality. The set (resp. ≥0, >0) is the real (resp. nonnegative, strictly positive) numbers, and $\boldsymbol{\mathrm{j}}$ is the imaginary unit. For $x\in^{n}$ and an index set $\mathcal{I}\subset\{1,\dots,n\}$, $[x_{i}]_{i\in\mathcal{I}}\in^{|\mathcal{I}|\times|\mathcal{I}|}$ is the diagonal matrix containing the appropriate elements of $x$. When no confusion can arise, we will simply write $[x]\in^{n\times n}$ for the diagonal matrix with $x$ on the diagonal. We let $\mathbbold{1}_{n}$ and $\mathbbold{0}_{n}$ denote the $n$-dimensional vectors of unit and zero entries, and $\mathbbold{0}$ is a matrix of all zeros of appropriate dimensions. The $n\times n$ identity matrix is $I_{n}$. The subspace $\mathbbold{1}_{n}^{\perp}\triangleq\{x\in^{n}\;|\;\mathbbold{1}_{n}^{\sf T}x=0\}$ is the subspace of n perpendicular to $\mathbbold{1}_{n}$. For $x\in^{n}$, we define the vector functions $\boldsymbol{\sin}(x)\triangleq(\sin(x_{1}),\ldots,\sin(x_{n}))^{\sf T}$, with $\boldsymbol{\arcsin}(x)$ and $\boldsymbol{\cos}(x)$ defined similarly, and for $x\in^{n}_{\geq 0}$, $\sqrt{x}=(\sqrt{x_{1}},\ldots,\sqrt{x_{n}})^{\sf T}$.

*Graphs and graph matrices : * A graph is a pair $(\mathcal{N},\mathcal{E})$, where $\mathcal{N}$ is the set of nodes and $\mathcal{E}\subseteq\mathcal{N}\times\mathcal{N}$ is the set of edges. If a label $e\in\{1,\dots,|\mathcal{E}|\}$ and an arbitrary direction is assigned to each edge $e=(i,j)\in\mathcal{E}$, the node-edge incidence matrix $A\in^{|\mathcal{N}|\times|\mathcal{E}|}$ is defined component-wise as $A_{ke}=1$ if node $k$ is the source node of edge $e$ and as $A_{ke}=-1$ if node $k$ is the sink node of edge $e$, with all other elements being zero. A graph is radial (acyclic, a tree) if it contains no cycles. For $x\in^{|\mathcal{N}|}$, $A^{\sf T}x\in^{|\mathcal{E}|}$ is the vector with components $x_{i}-x_{j}$, with $(i,j)\in\mathcal{E}$. We call $\mathrm{ker}(A)$ the cycle space of the graph. If the graph is connected, then $\mathrm{ker}(A^{\sf T})=\mathrm{im}(\mathbbold{1}_{|\mathcal{N}|})$, with $\mathrm{ker}(A)=\emptyset$ for acyclic graphs. In this radial case, for every $P\in\mathbbold{1}_{|\mathcal{N}|}^{\perp}$, there exists a unique $p\in^{|\mathcal{E}|}$ satisfying Kirchoff’s Current Law (KCL) $P=Ap$ [14, 15]. The vector $P$ is interpreted as nodal injections, with $p$ being the associated edge flows. If a weight $w_{ij}>0$ is assigned to each edge $(i,j)\in\mathcal{E}$, then ${\sf L}={\sf L}^{\sf T}=A[w]A^{\sf T}$ is the weighted Laplacian matrix for the graph, which is positive semidefinite with $\mathrm{ker}({\sf L})=\mathrm{im}(\mathbbold{1}_{|\mathcal{N}|})$.

Torus geometry: The set $\mathbb{S}^{1}$ denotes the unit circle, an angle is a point $\theta\in\mathbb{S}^{1}$, and $|\theta_{1}-\theta_{2}|$ denotes the geodesic distance between two angles $\theta_{1},\theta_{2}\in\mathbb{S}^{1}$. The $n$-torus $\mathbb{T}^{n}=\mathbb{S}^{1}\times\cdots\times\mathbb{S}^{1}$ is the Cartesian product of $n$ unit circles. For $\gamma\in[0,\pi/2)$ and a given graph $(\mathcal{N},\mathcal{E})$, let $\Theta(\gamma)=\{\theta\in\mathbb{T}^{|\mathcal{N}|}:\,\max_{(i,j)\in\mathcal{E}}|\theta_{i}-\theta_{j}|\leq\gamma\}$ be the closed set of angles $\theta=(\theta_{1},\dots,\theta_{n})$ with neighboring angles $\theta_{i}$ and $\theta_{j}$ no further than $\gamma$ apart.

II The Lossless Power Flow Model

This section introduces the power flow model used throughout the paper and states all modeling assumptions.

II-A Network and Branch Models

We model a steady-state synchronous power network as a connected, weighted, and undirected graph $(\mathcal{N},\mathcal{E})$ with nodes (buses) $\mathcal{N}$ and edges (branches) $\mathcal{E}\subseteq\mathcal{N}\times\mathcal{N}$. To each bus $i\in\mathcal{N}$ we associate a phasor voltage $V_{i}\angle\theta_{i}$, where $V_{i}>0$ is the bus voltage magnitude and $\theta_{i}\in\mathbb{S}^{1}$ is the voltage phase angle, and a complex power injection $P_{i}+\boldsymbol{\mathrm{j}}Q_{i}$. The real part $P_{i}$ is the active power injection, while $Q_{i}$ is the reactive power injection. There will be two types of buses in our network: we will have $n\geq 1$ load buses buses, denoted by the set $\mathcal{N}_{L}$, and $m\geq 1$ generator buses, denoted by the set $\mathcal{N}_{G}$. Without loss of generality, we order these buses as $\mathcal{N}_{L}=\{1,\ldots,n\}$ and $\mathcal{N}_{G}=\{n+1,\ldots,n+m\}$, with $\mathcal{N}=\mathcal{N}_{L}\cup\mathcal{N}_{G}$. Models for these buses are stated in Section II-B. This partitioning of buses $\mathcal{N}=\mathcal{N}_{L}\cup\mathcal{N}_{G}$ induces a partitioning of the branches222To keep notation under control, we will abbreviate $\mathbbold{1}_{|\mathcal{E}^{\ell\ell}|}$, $\mathbbold{0}_{|\mathcal{E}^{\ell\ell}|},I_{|\mathcal{E}^{\ell\ell}|}$ by $\mathbbold{1}_{\ell\ell},\mathbbold{0}_{\ell\ell},I_{\ell\ell}$, and similarly for the other edge sets.

[TABLE]

where, $\mathcal{E}^{\ell\ell}$ contains all branches between nodes $i,j\in\mathcal{N}_{L}$, $\mathcal{E}^{g\ell}$ contains all branches between generators and loads, etc. The incidence matrix $A\in^{(n+m)\times|\mathcal{E}|}$ of the graph $(\mathcal{N},\mathcal{E})$ inherits both the nodal and branch partitions, and may be written as the 2 $\times$ 3 block matrix

[TABLE]

For example, $A_{L}^{\ell\ell}$ is a matrix of size $n\times|\mathcal{E}^{\ell\ell}|$, mapping variables defined on the branches between load buses to variables defined at the load buses incident to those branches. The zero submatrices in (2) indicate that branches between generators cannot be incident to any load buses, and vice versa. Since the incidence matrix assigns an arbitrary orientation to each branch, we assume without loss of generality that all branches in $\mathcal{E}^{g\ell}$ are oriented from generators to loads: for each $(i,j)\in\mathcal{E}^{g\ell}$, $i\in\mathcal{N}_{G}$ and $j\in\mathcal{N}_{L}$. It follows that all non-zero elements of ${A}^{g\ell}_{L}$ are equal to $-1$, while all non-zero elements of $A^{g\ell}_{G}$ are equal to $+1$. Figure 1 illustrates these conventions.

We assume all transmission lines are inductive, and model them with the standard lumped parameter ${\Pi}$-model, which allows for the inclusion of inductive/capactive shunt loads, off-nomial tap ratios, and line-charging capacitors [16]. With this model, each branch $(i,j)\in\mathcal{E}$ is weighted by a purely imaginary admittance $y_{ij}=\boldsymbol{\mathrm{j}}b_{ij}$ where $b_{ij}\leq 0$.

We encode the admittances and grid topology in the bus admittance matrix $Y\in^{(n+m)\times(n+m)}$, with elements $Y_{ij}=-y_{ij}$ and $Y_{ii}=\sum_{j\neq i}^{n+m}y_{ij}+y_{\mathrm{shunt},i}$, where $y_{\mathrm{shunt},i}=\boldsymbol{\mathrm{j}}b_{\mathrm{shunt},i}$ is the shunt element at bus $i$. The susceptance $b_{\mathrm{shunt},i}$ and can be positive (capacitive) or negative (inductive). The susceptance matrix $B\in^{(n+m)\times(n+m)}$ is defined as the imaginary part of $Y$, and is characterized as follows [17].

Fact II.1** (Properties of the Susceptance Matrix)**

If the network contains no phase-shifting transformers, then

1. (i)

$B_{ij}=B_{ji}\geq 0$, with $B_{ij}>0$ if and only if $(i,j)\in\mathcal{E}$ or $(j,i)\in\mathcal{E}$; 2. (ii)

$B_{ii}=-\sum_{j=1,j\neq i}^{n+m}B_{ij}+b_{\mathrm{shunt},i}$* for all $i\in\mathcal{N}$.*

We vectorized the relevant variables as $V=(V_{1},\ldots,V_{n+m})^{\sf T}$, $\theta=(\theta_{1},\ldots,\theta_{n+m})^{\sf T}$, $P=(P_{1},\ldots,P_{n+m})^{\sf T}$, and $Q=(Q_{1},\ldots,Q_{n+m})^{\sf T}$. Like the incidence matrix, these vectors and the susceptance matrix $B$ inherit the partitioning of buses $\mathcal{N}=\mathcal{N}_{L}\cup\mathcal{N}_{G}$ as

[TABLE]

The principal submatrix $B_{LL}\in^{n\times n}$ describes the weighted interconnections between PQ buses and is central to our analysis; we refer to it as the grounded susceptance matrix, and impose the following standing assumption on it.

Assumption II.2** (Grounded Susceptance Matrix)**

The grounded susceptance matrix $B_{LL}$ in (3) is negative definite.

Assumption II.2 is usually satisfied in practical networks [18, Section III], and always satisfied in the absence of line-charging and shunt capacitors due to irreducibility and diagonal dominance of the susceptance matrix [19]. This assumption does not disallow shunt capacitors, but merely limits their size so that they do not “overcompensate” the network.

II-B Bus Models

Load Models: Each load bus $i\in\mathcal{N}_{L}$ is modeled as a PQ bus, with a fixed active power injection $P_{i}\in$ and a fixed reactive power injection $Q_{i}\in$. At PQ buses, voltage magnitudes $V_{i}>0$ and phases $\theta_{i}\in\mathbb{S}^{1}$ are free variables.

Generator Models: Each generator bus $i\in\mathcal{N}_{G}$ is modeled as a PV bus, with an active power injection $P_{i}\in$ fixed by the prime mover and constant voltage magnitude $V_{i}>0$ regulated on the network-side by an Automatic Voltage Regulator (AVR) system. At PV buses, phase angles $\theta_{i}\in\mathbb{S}^{1}$ and reactive power injections $Q_{i}\in$ are free variables. In the present work we do not consider generator reactive power constraints. We will therefore generally be unconcerned with reactive power injections at PV buses, as they are uniquely determined by (4b) given the other state variables. In other words, reactive power injections at PV buses may be considered as “outputs” of the power flow problem, rather than state variables.

The above bus models are uncontroversial when considering power flow solvability, and are the standard models [1] used by industry [20]. Since the network is lossless, no “slack bus” is required; all generators are treated as PV buses.333Equivalently, select PV bus $(n+m)$ as the slack bus, and its power injection $P_{n+m}=-\sum_{i=1}^{n+m-1}P_{i}$ is determined a priori.

II-C The Power Flow Equations

Denoting the voltage phasor at bus $i\in\mathcal{N}$ by $U_{i}=V_{i}e^{\boldsymbol{\mathrm{j}}\theta_{i}}$, the bus current injections $I\in\mathbb{C}^{n+m}$ in the network are given by the so-called nodal equations $I=YU=\boldsymbol{\mathrm{j}}BU$. The equations for balance of power may then be written as $P+\boldsymbol{\mathrm{j}}Q=[U]\mathrm{conj}(I^{*})$, where $\mathrm{conj}(\cdot)$ denotes complex conjugation. After expanding the complex exponentials into trigonometric terms and equating real and imaginary parts, one arrives at the celebrated lossless power flow equations

[TABLE]

written here in polar form, and where we have suppressed the reactive power equations for the PV buses (Section II-B). The equations (4a)–(4b) are a set of $2n+m$ nonlinear equations in the $n+m$ angles $\theta$ and the $n$ PQ bus voltage magnitudes $V_{L}$. As the angles enter only as differences, there are effectively only $n+m-1$ angles, and the problem appears to be overdetermined. However, note from (4a) that $\sum_{i=1}^{n+m}P_{i}=0$ for all choices of angles and voltage magnitudes; this equation states the balance of active power in a lossless network. In other words, (4a) is not an independent set of equations. Rather than eliminate one angle and one active power flow equation, we will retain both and impose the following assumption on the active power injections $P=(P_{1},\ldots,P_{n+m})^{\sf T}$.

Assumption II.3** (Balance Assumption)**

The necessary condition $P\in\mathbbold{1}_{n+m}^{\perp}$ is satisfied.

III Reformulation of the Power Flow Equations to Fixed-Point Power Flow

We now begin the formulation of our new fixed-point power flow model. Several key matrices — which we term the power network stiffness matrices — will appear very naturally during this reformulation process, and we now introduce the reader to them. These matrices quantify the strength of the grid, and are analogous to the stiffness matrices encountered in the theories of mechanical statics and vibrations.

III-A Open-Circuit Voltages and Stiffness Matrices

Letting $V_{G}=(V_{n+1},\ldots,V_{n+m})^{\sf T}$ denote the vector of generator voltage magnitudes, we introduce the following definition and characterization.

Definition 1** (Open-Circuit Load Voltages)**

The open-circuit load voltages $V_{L}^{*}\in^{n}$ are defined using the susceptance matrix (3) as

[TABLE]

Proposition III.1** (Open-Circuit Solution)**

Each component of $V_{L}^{*}$ is strictly positive. Moreover, when $P=\mathbbold{0}_{n+m}$ and $Q_{L}=\mathbbold{0}_{n}$, $(\theta,V_{L})=(\mathbbold{0}_{n+m},V_{L}^{*})$ is a solution of (4a)–(4b).

Proof:

. See appendix. ∎

When there are no shunt elements attached at PQ buses, the matrix $-B_{LL}^{-1}B_{LG}$ is a nonnegative row-stochastic matrix [21, Lemma II.1], and each open-circuit load voltage $V_{i}^{*}$ is therefore a weighted average of PV bus voltage set points. Capacitive shunt elements tend to push these open-circuit voltages up, while inductive shunt elements pull them down.

We now use these open-circuit voltages to define three matrices of interest, which have units of power.

Definition 2

(Laplacian, Branch, and Nodal Stiffness Matrices): The branch stiffness matrix $\mathsf{D}\in^{|\mathcal{E}|\times|\mathcal{E}|}$ is the diagonal positive definite matrix444For convenience we abuse our notation a bit, and here use $V_{i}^{*}$ also for the fixed PV bus voltages.

[TABLE]

The Laplacian stiffness matrix ${\sf L}={\sf L}^{\sf T}\in^{(n+m)\times(n+m)}$ is the symmetric positive semidefinite matrix

[TABLE]

With the partitioning of the susceptance matrix $B$ in (3), the nodal stiffness matrix $\mathsf{S}\in^{n\times n}$ is the symmetric, negative definite Metzler matrix

[TABLE]

To interpret ${\sf D}$, note from (4a) that $V_{i}V_{j}B_{ij}$ is the maximum active power transfer along the branch $(i,j)\in\mathcal{E}$. Thus, the branch stiffness matrix (6) captures the maximum branch-wise power transfers when voltages $V_{i}$ take their open-circuit values $V_{i}^{*}$. The Laplacian stiffness matrix (7), first introduced in [4], is a generalization of the branch matrix $\mathsf{D}$ for meshed networks. The nodal stiffness matrix was introduced by the author in [13]. Roughly speaking, $\mathsf{D}$ and ${\sf L}$ quantify how strong the branches of the network are, while $\mathsf{S}$ quantifies the interconnection strength between PQ buses. All three matrices depend on the open-circuit voltages (5), which in turn depend quite densely on the interconnection structure of the network, including shunt elements at PQ buses. The branch stiffness matrix $\mathsf{D}$ inherits the partitioning of the branches $\mathcal{E}$ in (1) as

[TABLE]

III-B Procedure for Deriving Fixed-Point Power Flow

The first major challenge one encounters in analyzing (4a)-(4b) is that variables appear as complicated products of trigonometric nonlinearities $\sin(\cdot)$ and $\cos(\cdot)$ with quadratic nonlinearities $V_{i}V_{j}$. More precisely, for any branch $(i,j)\in\mathcal{E}$ the product $V_{i}V_{j}$ will be a quadratic nonlinearity if $i$ and $j$ are PQ buses, linear if $i$ and $j$ are a PQ ​​-PV pair, or a constant if both buses are PV buses. We must therefore first figure out how to isolate these different combinations in a useful way. After this initial hurdle, our reformulation procedure will consist of three main steps:

1. Step 1: Isolate the phases in the active power flow (4a);
2. Step 2: Reformulate the $n$ reactive power flow equations (4b) into a form which partially isolates the phase angles;
3. Step 3: Combine the two reformulations by eliminating the phase angles from the reformulated reactive power flow, resulting in a single equation in the voltage variables, and rearrange the result into a fixed-point equation.

III-C Reformulation Step 1: Active Power Flow

We begin by introducing some additional notation associated with the incidence matrix (2). Due to the way (see Section II-A) that edge directions were assigned for the graph $(\mathcal{N},\mathcal{E})$, we may write $A$ as the difference between two nonnegative matrices $A=A(+)-A(-)$, i.e.,

[TABLE]

The matrix $A(+)$ indexes the buses at the sending end of each branch, while $A(-)$ indexes the corresponding receiving end buses. It follows that $A(+)^{\sf T}V$ is vector of sending-end voltages, while $A(-)^{\sf T}V$ is the vector of receiving end voltages. In vector form, we compute that

[TABLE]

For $(i,j)\in\mathcal{E}$, the quadratic products $V_{i}V_{j}B_{ij}$ can therefore be constructed from (11a)–(11b) through the formula

[TABLE]

It follows from (12) that the branch stiffness matrix $\mathsf{D}=[V_{i}^{*}V_{j}^{*}B_{ij}]_{(i,j)\in\mathcal{E}}$ in (6) has the vector representation

[TABLE]

and the submatrices in (9) may be easily identified by comparing to (13). We now find it useful to introduce a change of voltage variables. Using the open-circuit load voltages $V_{L}^{*}$ defined in (5), consider the bijective change of variables

[TABLE]

Thus, $v_{i}$ is the load bus voltage $V_{i}$ normalized by its open-circuit value $V_{i}^{*}$. Inserting the coordinate change (14) into (12), we find that

[TABLE]

where the map $h:^{n}\rightarrow^{|\mathcal{E}|}$ satisfies $h(\mathbbold{1}_{n})=\mathbbold{1}_{|\mathcal{E}|}$ and is defined by

[TABLE]

More explicitly, for any edge $e=(i,j)\in\mathcal{E}$, $h_{e}(v)$ satisfies

[TABLE]

In equation (15), we have decomposed the products $V_{i}V_{j}B_{ij}$ into a product of the branch stiffness matrix $\mathsf{D}$ and a nonlinear function $h(v)$. We now return to the active power flow (4a). In vector form, (4a) is written as

[TABLE]

as may be verified by applying the definition of the incidence matrix. Substituting our result (15) for $[V_{i}V_{j}B_{ij}]_{(i,j)\in\mathcal{E}}$ into the active power flow (18) we obtain

[TABLE]

Let $c\triangleq\mathrm{dim}(\mathrm{ker}(A))$ be the dimension of the cycle space for the graph $(\mathcal{N},\mathcal{E})$. In other words, $c$ is the number of independent cycles. Let $C\in^{|\mathcal{E}|\times c}$ be a matrix whose columns span $\mathrm{ker}(A)$.555One may always find a so-called totally unimodular basis for the cycle space $\mathrm{ker}(A)$, in which case one may take $C\in\{-1,0,1\}^{|\mathcal{E}|\times c}$ as the corresponding (oriented) edge-cycle incidence matrix [22, Section 3]. In this case, $C$ is totally unimodular [22, Theorem 3.4]; we proceed with this choice. Inspired by [4, SI Theorem 1], the following result shows that by using the Laplacian stiffness matrix ${\sf L}$, equation (19) can be solved for $\boldsymbol{\sin}(A^{\sf T}\theta)$ in terms of the voltages $v$ and a set of slack variables $y\in^{c}$ which account for loop flows of active power.

Lemma III.2

(Active Power Reformulation for Meshed Networks): Consider the active power flow equation (4a), and let $\psi:_{>0}^{n}\times^{c}\rightarrow^{|\mathcal{E}|}$ be defined by666Here, ${\sf L}^{\dagger}$ denotes the Moore-Penrose pseudoinverse of ${\sf L}$.

[TABLE]

The following two statements are equivalent:

1. (i)

$(\theta,V_{L})\in\Theta(\pi/2)\times^{n}_{>0}$* is a solution of (4a);* 2. (ii)

$(v,y)\in_{>0}^{n}\times^{c}$* is a solution of*

[TABLE]

with branch-wise angle differences \eta=A^{\sf T}\theta\,\,{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}(\boldsymbol{\mathrm{mod}}\,2\pi)}\in[-\frac{\pi}{2},\frac{\pi}{2}]^{|\mathcal{E}|} recovered via

[TABLE]

Proof:

. Our development so far shows that (4a) is equivalent to (19), so we show equivalence of (19) and (21)–(22). First, simply make a change of variables $\tilde{\psi}=\boldsymbol{\sin}(A^{\sf T}\theta)$. Then (19) may be written as the pair of equations

[TABLE]

The equation (23a) is linear in $\tilde{\psi}$, and we claim that $\psi(v,y)$ as defined in (20) is the general solution. Substituting (20) into (23a), we indeed find that

[TABLE]

where ${\sf L}{\sf L}^{\dagger}=\Pi$ is the projection matrix onto $\mathbbold{1}_{n+m}^{\perp}$, and we used the facts that $P\in\mathbbold{1}_{n+m}^{\perp}$ and that $AC=\mathbbold{0}$ by construction. Thus (20) is indeed the general solution to (23a), with the first term being a particular solution and the second term parameterizing the homogeneous solution in the slack variable $y\in^{c}$. Thus, (19) is equivalent to

[TABLE]

It remains only to show that (21) and (24) are equivalent. If $(\theta,V_{L})\in\Theta(\frac{\pi}{2})\times^{n}_{>0}$ is a solution of (4a), it follows from (24) that $\boldsymbol{\arcsin}(\psi(v,y))=A^{\sf T}\theta+2\pi k$ for some integer vector $k=(k_{1},\ldots,k_{|\mathcal{E}|})^{\sf T}\in\mathbb{Z}^{|\mathcal{E}|}$. Left-multiplying this equality by $C^{\sf T}$ and using that $C^{\sf T}A^{\sf T}=(AC)^{\sf T}=\mathbbold{0}$, we have

[TABLE]

Since $C\in\{-1,0,1\}^{|\mathcal{E}|\times c}$, $C^{\sf T}k\in\mathbb{Z}^{c}$, and therefore each component of the right-hand side of (25) is an integer multiple of $2\pi$. It follows that (25) is equivalent to (21). Conversely, if (21) holds, then $C^{\sf T}\boldsymbol{\arcsin}(\psi(v,y))=2\pi k^{\prime}$ for some integer vector $k^{\prime}=(k^{\prime}_{1},\ldots,k^{\prime}_{c})^{\sf T}\in\mathbb{Z}^{c}$. Since $C$ has full rank and is totally unimodular, we can always find a $k\in\mathbb{Z}^{|\mathcal{E}|}$ such that $C^{\sf T}k=k^{\prime}$ [23, Theorem 5.20]. The general solution to $C^{\sf T}\boldsymbol{\arcsin}(\psi(v,y))=2\pi k^{\prime}$ may therefore be written as $\boldsymbol{\arcsin}(\psi(v,y))=2\pi k+A^{\sf T}\theta$ for some $\theta\in^{n+m}$, with $2\pi k$ being the particular solution and $A^{\sf T}\theta$ parameterizing the homogeneous solution. Taking the $\boldsymbol{\sin}(\cdot)$ of both sides yields (24), which completes the proof.

∎

The modulo operation in Lemma III.2 is subtle, but is required to capture so-called loop flows; see [24, Remark 5.3.2] for details. Equations (21) and (22) are our desired reformulation of the active power flow equations (4a); (20) and (22) are essentially the explicit solution, while (21) enforces that Kichhoff’s voltage law holds true around the cycles of the network. The following corollary shows that in radial networks, the formula (20) which uses the Laplacian stiffness matrix ${\sf L}$ may be replaced by a simpler formula using the branch stiffness matrix $\mathsf{D}$.

Corollary III.3

(Active Power Reformulation for Radial Networks): If the graph $(\mathcal{N},\mathcal{E})$ describing the network is radial, then the cycle constraints (21) are discarded, and the explicit solution $\psi(v,y)$ from Lemma III.2 reduces to

[TABLE]

where

[TABLE]

are the unique branch-wise active power flows, satisfying Kirchhoff’s current law $P=Ap$.

Proof:

. That the cycle constraints (21) can be discarded is clear, as there are no cycles by assumption, so we may set $Cy=\mathbbold{0}_{|\mathcal{E}|}$ in (20). Comparing (20) and (26) then, we need only prove that $A^{\sf T}{\sf L}^{\dagger}P=\mathsf{D}^{-1}p$. To do so, set $Z\triangleq A^{\sf T}{\sf L}^{\dagger}P$, and note that

[TABLE]

Since $\mathrm{ker}(A)=\emptyset$, we may left-multiply by the left-inverse $(A^{\sf T}A)^{-1}A^{\sf T}$ of $A$ and insert (27) to obtain

[TABLE]

which shows the result. ∎

For later use, we use $h(v)$ from (16) to expand out the acyclic solution (26) in terms of the branch partitions as

[TABLE]

where

[TABLE]

are the branch-wise phase angle differences.

III-D Reformulation Step 2: Reactive Power Flow

It is convenient to define an “unoriented” version of the incidence matrix $A$, denoted by $|A|\in^{(n+m)\times|\mathcal{E}|}$, with all non-zero elements set to $+1$. From (10), it is clear that $|A|=A(+)+A(-)$. As we did for the incidence matrix (2), we partition $|A|$ according to bus and branch types as

[TABLE]

We now focus on the reactive power flow equations (4b). Using Lemma A.3, (4b) can be written in vector form as

[TABLE]

Adding and subtracting $|A|_{L}[V_{i}V_{j}B_{ij}]_{(i,j)\in\mathcal{E}}\mathbbold{1}_{|\mathcal{E}|}$ from the right-hand side, (31) becomes

[TABLE]

Applying Lemma A.3 (i)–(iii) to the first two terms, we obtain

[TABLE]

Working on the first term in (LABEL:Eq:ReactiveLOL), by substituting for $V_{L}=[V_{L}^{*}]v$ from the coordinate change (14) and simplifying using the nodal stiffness matrix $\mathsf{S}$ from (8), we find that

[TABLE]

For the second term in (LABEL:Eq:ReactiveLOL), again introduce $\eta=A^{\sf T}\theta$ and insert the formula (15) for $[V_{i}V_{j}B_{ij}]_{(i,j)\in\mathcal{E}}$ to obtain

[TABLE]

which is the main result of Step 2. The following lemma summarizes our reformulations.

Lemma III.4** (Reactive Power Reformulation)**

Consider the reactive power flow equation (4b). The following two statements are equivalent:

1. (i)

$(\theta,V_{L})\in\Theta(\pi/2)\times^{n}_{>0}$* is a solution of (4b);* 2. (ii)

$(\eta,v)\in[-\frac{\pi}{2},\frac{\pi}{2}]^{|\mathcal{E}|}\times_{>0}^{n}$* is a solution of (33) .*

III-E Reformulation Step 3: Eliminate the Phase Angles

We now combine the two reformulations (22) and (33). To begin, we apply $\sin^{2}\eta_{ij}+\cos^{2}\eta_{ij}=1$ component-wise to (22) and to obtain

[TABLE]

where we have selected the positive square root for all components, as we are interested in phase angle vectors $\theta\in\Theta(\pi/2)$. Next, left-multiply both sides of (33) by $-\frac{1}{4}\mathsf{S}^{-1}[v]^{-1}$ and rearrange to obtain

[TABLE]

Substituting (34) into (35), we may combine Lemma III.2 and Lemma III.4 to obtain our most general main modeling result.

Theorem III.5

(Fixed-Point Power Flow for Meshed Networks): Consider the coupled power flow equations (4a)–(4b). The following two statements are equivalent:

1. (i)

$(\theta,V_{L})\in\Theta(\pi/2)\times^{n}_{>0}$* is a solution of (4a)–(4b);* 2. (ii)

$(v,y)\in_{>0}^{n}\times^{c}$* is a solution of the fixed-point power flow*

[TABLE]

where

[TABLE]

where $h(v)$ is as in (16), with the angular differences \eta=A^{\sf T}\theta\,\,\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}(\boldsymbol{\mathrm{mod}}\,2\pi)\in[-\frac{\pi}{2},\frac{\pi}{2}]^{|\mathcal{E}|} recovered via $\eta=\boldsymbol{\arcsin}(\psi(v,y))$ .

The model (36a)–(36b) depends only on the scaled PQ bus voltage magnitudes $v\in^{n}_{>0}$ and on the slack variables $y\in^{c}$ which enforce the cycle constraints arising from Kirchhoff’s voltage law. Phase angles are completely absent, and are recovered as “outputs” by applying $\arcsin(\cdot)$ component-wise to (37b). The terminology fixed-point power flow comes from the special form that the model takes for radial networks.

Corollary III.6

(Fixed-Point Power Flow for Radial Networks I): Consider the coupled power flow equations (4a)–(4b) and assume that the graph $(\mathcal{N},\mathcal{E})$ describing the network is radial. The following two statements are equivalent:

1. (i)

$(\theta,V_{L})\in\Theta(\pi/2)\times^{n}_{>0}$* is a solution of (4a)–(4b);* 2. (ii)

$v\in^{n}_{>0}$* is a fixed point777A fixed point of $f_{\rm radial}$ is a point $v\in^{n}_{>0}$ satisfying $f_{\rm radial}(v)=v$. of the mapping $f_{\rm radial}:^{n}_{>0}\rightarrow^{n}$ defined by*

[TABLE]

where

[TABLE]

with $h(v)$ as in (16) and the angular differences $\eta=A^{\sf T}\theta\in[-\frac{\pi}{2},\frac{\pi}{2}]^{|\mathcal{E}|}$ recovered via $\eta=\boldsymbol{\arcsin}(\psi(v))$ .

Proof:

. By Corollary III.3, $\psi(v,y)$ from (37b) reduces to $\psi(v)$ in (39b), and hence $u(v,y)$ from (37a) reduces to $u(v)$ in (39a). Again by Corollary III.3, we may omit the cycle constraint (36b). It follows that $f_{\rm mesh}(v,y)$ as defined in (36a) reduces to $f_{\rm radial}(v)$ as in (38). ∎

The function $f_{\rm radial}(v)$ may be written in an alternative form which emphasizes the importance the branch partitioning $\mathcal{E}=\mathcal{E}^{\ell\ell}\cup\mathcal{E}^{g\ell}\cup\mathcal{E}^{gg}$. This form is less compact, but more explicit, and will be used for analysis purposes in Part II.

Corollary III.7

(Fixed-Point Power Flow for Radial Networks II): With the notation $u(v)=(u_{\ell\ell}(v),u_{g\ell}(v),u_{gg}(v))^{\sf T}$, an equivalent expression for the fixed point map $f_{\rm radial}(v)$ in Corollary III.6 is

[TABLE]

Proof:

. See appendix. ∎

III-F Discussion of FPPF Derivation

The fixed-point power flow model (36) is parameterized by the stiffness matrices $\mathsf{D},\mathsf{L}$, and $\mathsf{S}$ and the power demands $P$ and $Q_{L}$. As discussed in Section III-A, the stiffness matrices encode both the topology of the network and the values of relevant parameters such as line susceptances and generator voltages. A particularly important observation is that the active and reactive power loads $P$ and $Q_{L}$ each enter the model multiplied by the inverse of a stiffness matrix ($A^{\sf T}{\sf L}^{\dagger}P$ in the first case, $\mathsf{S}^{-1}[Q_{L}]$ in the second case). The FPPF model reveals that these very specific combinations of network parameters and loads are the important quantities to focus on.

In the radial case, the map $f_{\rm radial}$ in (40) consists of five distinct, easily interpretable terms:

• •

the first term is the open-circuit voltage profile of the network in the scaled variables $v_{i}=V_{i}/V_{i}^{*}$;

• •

the second term is a reformulation of the decoupled reactive power flow equation; see [13].

• •

the third term proportional to $u_{g\ell}(v)$ accounts for the influence of active power branch flows between PV and PQ buses on the voltage magnitudes at PQ buses.

• •

the fourth and fifth terms proportional to $u_{\ell\ell}(v)$ account for the influence of active power branch flows between PQ buses on PQ bus voltage magnitudes. There are two terms to account for this because a branch-wise flow of active power between two PQ buses affects the voltage magnitude at both ends of the branch.

A more detailed examination of the internal structure of these terms is deferred to Part II [25]. Conceptually, the FPPF model is similar Baran-Wu branch flow model [5] in that phase angles are absent. In contrast to the branch flow model, where sending-end branch power flows and currents are state variables, the voltages $v$ are the only state variables in (40); this difference in complexity can perhaps be attributed to the lossless character of the network considered here.

III-G Approximate Solution to Lossless Power Flow

We now leverage the FPPF model to derive an explicit approximate solution to the power flow equations (4a)–(4b). By construction, the FPPF model (36) has the property that under open-circuit conditions (when $Q_{L}=\mathbbold{0}_{n}$ and $P=\mathbbold{0}_{n+m}$), it holds that $(v,y,\eta)=(\mathbbold{1}_{n},\mathbbold{0}_{c},\mathbbold{0}_{|\mathcal{E}|})$ is a solution:

[TABLE]

This is the solution we desire in practice: a high-voltage solution ($v\approx\mathbbold{1}_{n}$) with small angular differences ($\eta\approx\mathbbold{0}_{|\mathcal{E}|}$). By Taylor expanding each side of (36) around this solution, we can obtain an explicit expression for the “linearized” power flow solution $(\eta_{\rm lin},v_{\rm lin})$ to lowest order in both $Q_{L}$ and $P$. Expanding (36b) around the open-circuit solution, to first order we find that

[TABLE]

By construction, $C$ is full column rank with $\mathrm{Im}(C)=\mathrm{ker}(A)$ (Section III-C). It follows that $AC=\mathbbold{0}$ and that $C^{\sf T}{\sf D}^{-1}C$ is invertible. We therefore conclude that $y_{\rm lin}=\mathbbold{0}_{c}$ is the approximate solution for the cyclic slack variables. While we omit the details, expanding both (36a) and the expression $\eta=\boldsymbol{\arcsin}(\psi)$ to lowest order in $P$ and $Q_{L}$ yields

[TABLE]

The first two terms in (41a) are the approximate solution for decoupled reactive power flow, as obtained in [26, 13]. The third term is novel, and captures how — to lowest order — active power injections affect voltage magnitudes quadratically. The linearization (41b) is the standard DC Power Flow approximation [12]. Taken together, (41) is an approximate solution to the lossless power flow equations. Intuitively, this approximation will be accurate when (i) loop flows of power are insignificant, and (ii) the quantities $\|{\sf S}^{-1}Q_{L}\|$ and $\|A^{\sf T}{\sf L}^{\dagger}P\|$ are both small.888In Part II we will show for radial networks that appropriately restricting these two quantities guarantees the existence of a power flow solution. We note that higher-order explicit approximations can be obtained by retaining higher-order terms in the Taylor expansion.

IV Numerical Experiments

We now present simulation studies on standard test cases to test the computational performance of the FPPF model as well as the accuracy of the approximate solution (41). Algorithm 1 describes the implementation of the FPPF used for the tests described here; many variations and computational refinements are of course possible. The basic approach is to iterate the mapping $v_{k+1}=f_{\rm mesh}(v_{k},y_{k})$ to update the scaled voltages $v$. For meshed networks the slack variables $y$ must also be updated. While many options are possible, we will use a Newton step

[TABLE]

where

[TABLE]

is the Jacobian matrix of (36b). After a desired relative tolerance $\epsilon$ is reached, the algorithm terminates.

Remark 1

(Computational Remarks): Various combinations of constants such as $A^{\sf T}{\sf L}^{\dagger}P$, ${\sf D}^{-1}C$, and $|A|_{L}{\sf D}$ appear in the iterations of Algorithm 1; these can be precomputed then stored for future use. For example, one may solve the sparse linear equation ${\sf L}z=P$, calculate $A^{\sf T}z$, then store the result. The iteration $v\leftarrow f_{\rm mesh}(v,y)$ should be computed by right-multiplying (36a) by ${\sf S}$ and solving the corresponding sparse linear system. A factorization (e.g., Cholesky) of ${\sf S}$ can be stored, and each update will then require only one forward-backward substitution; other sparsity-exploiting techniques could also be applied. For the Newton step (42), the Jacobian (43) is sparse and the variable portion $\left(I_{|\mathcal{E}|}-[\psi_{k}]^{2}\right)^{-\frac{1}{2}}[h(v_{k})]^{-1}$ is diagonal. In summary, each iteration of Algorithm 1 requires the solution of an $n\times n$ system of sparse equations (with a constant coefficient matrix) and the solution of a $c\times c$ system of sparse equations for the slack Newton step. $\square$

Remark 2

(Dimensional Comparison To Standard Power Flow Methods): Classic Newton-Raphson or Fast Decoupled Load Flow implementations iterate on the $n+m$ phase angles and $n$ voltage magnitudes. In contrast, the Algorithm 1 iterates on the $n$ voltage magnitudes and $c$ slack variables. Typically $c<n+m$, and therefore the FPPF algorithm requires the solution of (sometimes substantially) smaller systems of linear equations than these implementations. For example, the 39 bus New England system has $c=8$, the 2383 Polish system has $c=514$, and the 9241 PEGASE system has $c=6809$. We also note that many different solution techniques (including Newton’s method) could be applied to the system of $n+c$ nonlinear equations (36). Algorithm 1 is just one option, but is well-motivated by our development in Part II. $\square$

We apply Algorithm 1 and the approximate solution (41) to the standard MATPOWER test cases [27, 28]. In all cases, branch and shunt conductances were set to zero, in line with our main theoretical assumption.999For context on this assumption, the mean branch $R/X$ ratios of the networks in Table I are approximately 0.3, 0.2, 0.4, 0.15, 0.4, 0.1, 0.3, 0.2, 0.1, 0.3, 0.2, and 0.2. The mean differences between the voltage solutions computed with and without resistances are 0.007, 0.006, 0.010, 0.004, 0.024, 0.006, 0.003, 0.010, 0.005, 0.005, 0.005, and 0.007 per unit, respectively. The cycle-space matrix $C\in\{-1,0,1\}^{|\mathcal{E}|\times c}$ was generated using C = null(A,"r") in MATLAB. Simulation results are shown in Table I. The second column shows the number of iterations required by Algorithm 1 to reach a relative solution tolerance of $\epsilon=0.001$ p.u. on the voltage magnitudes. As can be seen, the iterations converge to the high-voltage solution quickly. The number of iterations is essentially independent of system size, and in all cases the solution agrees with the one determined by MATPOWER. Columns 3 and 4 show the max and mean prediction errors

[TABLE]

between the exact voltage magnitude solution $V_{L}$ and the approximate solution determined by (41), in per unit. In most of the cases the maximum error is quite small, with larger test cases showing larger maximum errors. Across all cases however, the mean error is consistently quite small, indicating a good overall approximation.

To examine performance in more heavily loaded networks, each case was loaded along the base case direction 90% of the way to the power flow insolvability boundary, as determined by continuation power flow cpf in MATPOWER. The previous experiments were repeated, and the results are shown in columns four and five of Table I. As would be expected, convergence to the power flow solution now takes more iterations, and the approximate solution (41) yields less accurate predictions; mean error $\delta_{\rm avg}$ is smaller, but omitted.

For base case loading, the voltage profile results for the 39 bus system are plotted in Figure 2. As stated, the profile obtained by iterating the FPPF (blue crosses) coincides with the solution determined by MATPOWER (black circles). The approximate solution (41) is plotted in red. The approximate solution is quite accurate, but systematically overestimates the voltage values; this behaviour is consistent across all test cases.

For the base case IEEE 300 bus test system, Figure 3 compares the convergence of the FPPF to the stock implementations of Newton-Raphson (NR) and Fast-Decoupled Load Flow (FDLF) from MATPOWER. All three solvers were initialized from a flat start; the vertical axis is the $2$-norm of the difference between current iterate of voltage magnitudes and the exact voltage magnitude solution vector. Both the FPPF and the FDLF show linear convergence (with slightly different rates) while the NR iterates converge quadratically. Figure 4 repeats the comparison for heavy loading. In this case the NR requires only one additional iteration to reach machine precision, while the FPPF and FDLF iteration counts double; this is consistent across all cases. We conclude that the FPPF algorithm convergence is similar or slightly favourable compared to FDLF, but does not approach that of NR.

The linear convergence of the FPPF is consistent with our results in Part II, where we will show — for a subclass of radial networks — that the FPPF is a contraction mapping whenever there exists a high-voltage solution. Those theoretical results do not generalize to the cases considered numerically here. However, they do suggest that contractivity of the FPPF and the existence of a high-voltage solution will go hand-in-hand.

As a final test, we examine the sensitivity of Algorithm 1 to initialization. If the FPPF is a contraction mapping on a large subset of voltage-space, we would expect the convergence of Algorithm 1 to be insensitive to the choice of initialization. We consider the IEEE 118 bus test system at base case loading. We fix the angle initial condition at $\theta=\mathbbold{0}_{n+m}$, and generate 1000 random voltage magnitude initial conditions, each with components pulled uniformly from the interval $[1-\alpha,1+\alpha]$ for various values of $\alpha$. For each initial condition, we run the NR, FDLF, and FPPF algorithms. If the iterates converge to the known high-voltage solution, we mark the test successful; otherwise, we say the solver has failed. Table II shows the fraction of successful tests for various values of $\alpha$.

For very small values of $\alpha$ (i.e., initial conditions close to a flat start) all three solvers behave similarly. As $\alpha$ increases, both the NR and FDLF solvers increasingly struggle to calculate the high-voltage solution; for some initializations these solvers converge to a low-voltage solution, for other initializations they diverge. On the larger test systems, the NR and FDLF are even more sensitive than suggested by Table II. The results show that for NR and FDLF, even one mildly perturbed component of an otherwise tolerable initial condition causes failure to converge to the high-voltage solution, In contrast, the FPPF iteration recovers the high-voltage solution from nearly every constructed initial condition.

V Conclusions

We have developed a new formulation of lossless power flow equations, which we term the fixed-point power flow. The model is naturally parameterized in terms of the power network stiffness matrices, which concisely encode the topology and parameters of the network, and leads immediately to an approximate solution of the power flow equations. We then proposed an algorithm for solving the power flow equations based on the FPPF. Numerical testing shows that this algorithm converges quickly for base case and stressed conditions. The convergence is linear, at a rate comparable or favorable compared to fast-decoupled methods, and is very insensitive to initialization.

In Part II [25] of this paper, we restrict ourselves to the case of radial networks, and leverage the fixed-point power flow to derive sufficient and tight conditions for the existence and uniqueness of a stable power flow solution. The analysis presented there will also provide guarantees for the convergence of the FPPF iteration $v_{k+1}=f_{\rm radial}(v_{k})$ for a subclass of radial networks. Future directions and open problems are deferred to the conclusions of Part II.

Appendix A Supporting Results and Proofs

Proof of Proposition III.1:   That the definition (5) is well-posed follows from Assumption II.2. From Assumption II.2 and Fact II.1, we conclude that $-B_{LL}$ is a symmetric positive definite $M$-matrix.101010A matrix $A\in^{n\times n}$ is a $Z$-matrix if $A_{ij}\leq 0$ for all $i\neq j$. A $Z$-matrix $A$ is a nonsingular $M$-matrix if it can be expressed as $A=sI_{n}-B$, where $B\in^{n\times n}$ has nonnegative elements and $s>\rho(B)$, where $\rho(B)$ is the spectral radius of $B$ [29, Chapter 6]. If $A$ is a nonsingular $M$-matrix, then the elements of $A^{-1}$ are nonnegative [29, Chapter 6, Theorem 2.3, $\mathbf{N}_{38}$]. Hence $-B_{LL}^{-1}$ is nonnegative. The submatrix $B_{LG}$ is nonnegative by Fact II.1, and $V_{G}$ is strictly positive. Since the network is connected, $B_{LG}$ contains at least one non-zero positive element in each row and column, and it follows that $V_{L}^{*}$ is strictly positive. Now suppose that $P=\mathbbold{0}_{m+m}$ in (4a) and $Q_{L}=\mathbbold{0}_{n}$ in (4b). Using Lemma (iii), substituting $\theta=\mathbbold{0}_{n+m}$ into (44), and applying Lemma (iii)(iii), it follows that $(\theta,V_{L})=(\mathbbold{0}_{n+m},V_{L}^{*})$ is a solution of (4b). Similarly, substituting $\theta=\mathbbold{0}_{n+m}$ into (4a) shows that $(\mathbbold{0}_{n+m},V_{L}^{*})$ also solves (4a). This completes the proof. $\square$

Lemma A.1** (Partitioned Incidence Matrix II)**

Consider the incidence matrix (2) along with its plus/minus decomposition (10). For $v\in_{>0}^{n}$, let $v^{-1}=(v_{1}^{-1},\ldots,v_{n}^{-1})^{\sf T}$. The following identities hold:

1. (i)

$[A_{L}^{g\ell}(-)^{\sf T}v]^{-1}\mathbbold{1}_{g\ell}=A_{L}^{g\ell}(-)^{\sf T}v^{-1}=(|A_{L}^{g\ell}|)^{\sf T}v^{-1}$** 2. (ii)

$[A_{L}^{\ell\ell}(-)^{\sf T}v]^{-1}\mathbbold{1}_{\ell\ell}=A_{L}^{\ell\ell}(-)^{\sf T}v^{-1}$** 3. (iii)

$[A_{L}^{\ell\ell}(+)^{\sf T}v]^{-1}\mathbbold{1}_{\ell\ell}=A_{L}^{\ell\ell}(+)^{\sf T}v^{-1}$** 4. (iv)

$[v]^{-1}A_{L}^{\ell\ell}(+)=A_{L}^{\ell\ell}(+)[A_{L}^{\ell\ell}(+)^{\sf T}v^{-1}]$** 5. (v)

$[v]^{-1}A_{L}^{\ell\ell}(-)=A_{L}^{\ell\ell}(-)[A_{L}^{\ell\ell}(-)^{\sf T}v^{-1}]$** 6. (vi)

$[v]^{-1}A_{L}^{g\ell}(-)=A_{L}^{g\ell}(-)[A_{L}^{g\ell}(-)^{\sf T}v^{-1}]$* .*

Proof:

. (i) By construction $A_{L}^{g\ell}(-)\in\{0,1\}^{n\times|\mathcal{E}^{g\ell}|}$ has exactly one element equal to one in each column; if the column corresponds to branch $(i,j)\in\mathcal{E}^{g\ell}$, then the non-zero element is in row $j$. Thus, $A_{L}^{g\ell}(-)^{\sf T}v^{-1}$ produces a branch vector with entry $v_{j}^{-1}$ in the entry corresponding to branch $(i,j)\in\mathcal{E}^{g\ell}$. It follows that $[A_{L}^{g\ell}(-)^{\sf T}v][A_{L}^{g\ell}(-)^{\sf T}v^{-1}]=I_{g\ell}$, and the result follows by right-multiplying by $\mathbbold{1}_{\ell\ell}$ and rearranging. The third equality is trivial as $A_{L}^{g\ell}(-)=|A|_{L}^{g\ell}$ by definition. The proofs of (ii)-(v) are analogous.

∎

Lemma A.2** (Identities)**

The following identities hold:

1. (i)

$\!\begin{aligned} L_{1}(v)&\triangleq[v]^{-1}|A|_{L}^{\ell\ell}[h_{\ell\ell}(v)]\\ &=A_{L}^{\ell\ell}(+)[A_{L}^{\ell\ell}(-)^{\sf T}v]+A_{L}^{\ell\ell}(-)[A_{L}^{\ell\ell}(+)^{\sf T}v]\,,\end{aligned}$** 2. (ii)

$L_{2}(v)\triangleq[v]^{-1}|A|_{L}^{g\ell}[h_{g\ell}(v)]=|A|_{L}^{g\ell}$* .*

Proof:

. (i): Using the decomposition $|A|_{L}^{\ell\ell}=A_{L}^{\ell\ell}(+)+A_{L}^{\ell\ell}(-)$ from (30) and inserting the expression for $h_{\ell\ell}(v)$ from (16), we find that

[TABLE]

Using Lemma A.1 (iv) and (v) to substitute for the first factor in each term, then reshuffling the diagonal matrices, we find that

[TABLE]

from which the result follows.

(ii): Similar to (i), first substitute for $|A|_{L}^{g\ell}$ from (30) and for $h_{g\ell}(v)$ from (16) to find

[TABLE]

Applying Lemma A.2 (vi) to the first two terms in the product, we obtain

[TABLE]

which completes the proof. ∎

Lemma A.3

(Alternate Expressions for the Reactive Power Flow Equation): The reactive power flow (4b) can be written in vector form as

[TABLE]

Moreover, the following expressions are all equal:

1. (i)

$-[V_{L}][B_{ii}]_{i\in\mathcal{N}_{L}}V_{L}-|A|_{L}[V_{i}V_{j}B_{ij}]_{(i,j)\in\mathcal{E}}\mathbbold{1}_{|\mathcal{E}|}$* ,* 2. (ii)

$-[V_{L}]\left(B_{LL}V_{L}+B_{LG}V_{G}\right)$* ,* 3. (iii)

$-[V_{L}]B_{LL}\left(V_{L}-V_{L}^{*}\right)$* .*

Proof:

. For $i\in\mathcal{N}_{L}$, the $i$th component of (44) is given by

[TABLE]

which is exactly (4b). To obtain (i), simply set $A^{\sf T}\theta=\mathbbold{0}_{|\mathcal{E}|}$ in (44). The equivalence of (i) and (ii) follows most easily by observing that both are equivalent to (4b) by setting $\cos(\theta_{i}-\theta_{j})=1$ in (4b) and separating the sum in the appropriate fashion to obtain (i) or (ii). The equivalence of (ii) and (iii) is immediate by invoking Assumption II.2 and inserting (5).

∎

Proof of Corollary III.7:   Let $k(v)$ denote the final term in (38). Expanding $k(v)$ by substituting the block incidence matrix (2), the block $\mathsf{D}$ matrix (9), and the partitioned $h(v)$ from (16), we obtain

[TABLE]

where we have rearranged some diagonal matrices, and where by combining (39a) and (39b),

[TABLE]

Lemma A.2(ii) shows that $L_{2}(v)=|A|_{L}^{g\ell}$ independent of $v$, so the previous simplifies to

[TABLE]

Applying Lemma A.2(i) to explicitly expand $L_{1}(v)$ leads immediately to (40). $\square$