Log(v) 3LPF: A Linear Power Flow Formulation for Unbalanced Three-Phase Distribution Systems

In this work, we introduce Log(v) 3LPF, a linear power flow solver for unbalanced three-phase distribution systems. Log(v) 3LPF uses a logarithmic transform of the voltage phasor to linearize the AC power flow equations around the balanced case. We incorporate the modeling of ZIP loads, transformers, capacitor banks, switches and their corresponding controls and express the network equations in matrix-vector form. With scalability in mind, special attention is given to the computation of the inverse of the system admittance matrix, Ybus. We use the Sherman-Morrison-Woodbury identity for an efficient computation of the inverse of a rank-k corrected matrix and compare the performance of this method with traditional LU decomposition methods using Big-$\mathcal {O}$ notation. We showcase the solver for a variety of network sizes, ranging from tens to thousands of nodes, and compare the Log(v) 3LPF with commercial-grade software, such as OpenDSS.

Backward-Forward Sweep [4], [5], were presented as a solution for modern computers to the power flow problem. Early on, it was noted that convergence is not guaranteed in the Newton-Raphson method [2], [6], [7] and the study of this problem is still an active area of research. The initialization, the problem conditioning or the voltage characteristics may impact the convergence of these algorithms [8]. Multi-phase modelling, introduced in [5] and also studied in [9], [10], is often required in distribution systems, as opposed to single-phase positive sequence, to capture the unbalanced characteristics of distribution networks. Relatively recently, exact convex relaxation techniques for single-phase radial systems [11], [12] have been developed to address the non-convex characteristics of the AC Optimal Power Flow (OPF) problem, originally treated in [13]. These techniques may also work when solving single-phase PF problems (c.f. [11]). However, the proposed convex relaxation relies on the radial assumption and is not necessarily efficient, since the number of variables grows with the square of the number of nodes. The radial assumption is invalid in three-phase systems since the matrix of voltages is indeed not rank 1, unlike in single-phase systems, due to the meshed characteristics of the network. The connection of two multi-phase buses represents a complete subgraph. The caveats of the convex relaxation are discussed in [14].
Since the AC power flow problem was first formulated, multiple linearizations that remove the non-convexities and non-linearities have been proposed. The most widely used linearization relies on the DC assumptions [15], where voltages are nominal and the resistance of the line is neglected. The DC assumption has been proven to provide fast results at the cost of accuracy [16], specially in distribution systems where the voltage may deviate significantly from the nominal values. In the context of radial distribution systems, Baran and Wu [17] proposed a single-phase AC PF formulation, linear in the voltage and current squared. Other single-phase AC PF linearizations were proposed in [18]- [20]. The DistFlow equations were later expanded to three-phase unbalanced networks [21] while considering line losses. However, the linear DistFlow equations are limited to radial networks and the inclusion of losses in the formulation requires training of complex models. Also, the modeling of different network elements presents a challenge since the DistFlow equations are linear in the squared of the voltage and current magnitudes, and they do not consider voltage angles. New linearization techniques for multiphase network are 0885-8950 © 2022 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission.
shown in recent work [22]- [32]. For instance, the work from [22] assumes nominal voltages and no line losses. The authors in [23] use a single iteration of the fixed-point iteration method, and [24] uses curve fitting to approximate the nonlinear voltage ratios assuming small changes. Also, the work from [26] uses a generalized form of the LinDistFlow equations that does not consider voltage angles and where the authors linearize the cross-products of the power flows and remove the dependence of downstream voltages on the upstream voltage. The work from [27] assumes small angle differences, which fails to consider transformers that may introduce 30 or 180°lead or lag phase shift. In [28], the authors compare different linearization methods and concluded that first-order Taylor methods best approximate the voltage phasors. The authors from [29]- [32] used data-driven and regression approaches to derive linear expression for the power flow problem. The work from [22]- [24], [26]- [32] did not consider the modeling of network elements such as multiwinding multi-phase transformers and its multiple connections, or capacitor banks, switches, and the corresponding controls of all these devices. Additionally, some of these methods [22], [23], [26], [27] are tested on small networks and their ability to scale to realistic feeders is a question that remains unanswered. Furthermore, some linearizations rely on the radial assumption, and do not point to fast inversion methods when changes in the boundary conditions of the network tend to be sparse. In this work, we provide an open-source three-phase unbalanced linear power flow solver, the Log(v) 3LPF, that uses the AC PF equations derived using first principles. We take inspiration from previous single-phase work in [33] and expand it to use a logarithmic transform of the voltage phasors as a way of removing other non-linearities that may arise due to different load models or controllers. The methodology to deal with the three-phase model as well as the inclusion of ZIP models and controllers the fast matrix inversion method we propose are all original contributions compared to [33]. Expressing the voltage phasor magnitude in logarithmic form enables a first-order Taylor expansion around zero of log-voltage magnitude and log-voltage angle differences. We work under the assumption that a linear function centered around the balanced case will yield a good approximation. Moreover, the logarithmic transform technique has applications beyond radial networks, and could be generalized [33] for three-phase transmission networks. A linearization provides some convergence guarantees that are well-known issues in iterative solvers [2], [6]- [8]. In addition to linearizing the power flow equations, we provide linear expressions to model constant-impedance (Z), constant-current (I) and constant-power (P) ZIP loads, transformers and its multiple configurations, capacitors banks and switches, as well as control functionality for some of these devices. In general, changes in the controls are modeled as a perturbation of the system admittance matrix, Y bus . This new formulation is presented in a compact matrix-vector form, thus enabling further understanding of the impacts of changing boundary conditions to the network through a simple linear operator. We also provide attention to the efficiency of the solver, specially in the computation of the inverse of Y bus . Since changes in the operating condition of the network tend to be sparse, we leverage the Sherman-Morrison-Woodbury identity for an efficient computation of the inverse. We show results in terms of flops and provide a comparison with traditional Lower-Upper (LU) decomposition methods that are used by commercial software [34]. Finally, we show results from running different IEEE radial networks, ranging from tens to thousands of nodes, and validate our solution against the OpenDSS software [35], that relies on iterative methods to provide an exact solution. A thorough validation of any linearization provides additional value, and is the basis to understand potential use cases. In fact, this work is motivated by the need of algorithms that can solve thousands of PF problems efficiently, to train reinforcement learning (RL) algorithms and learn control policies for the distribution system [36]. Potential applications are broad, ranging from cybersecurity [37] to applications on coupled transmission and distribution studies [38].

II. MODELING OF MULTI-PHASE LINE ELEMENTS
A distribution feeder can be modeled as an undirected graph G = {N , E} with |N | buses, and |E| number of edges. We denote Ω n as the set of phases available in bus n ∈ N , and Ω nm is the set of phases on each edge e = {n, m} ∈ E connecting two buses n, m ∈ N . It should be noted that |Ω nm | = |Ω mn |, and |Ω nm | ≤ min(|Ω n |, |Ω m |) ≤ 3. We further divide the set of edges E into the set of distribution lines E l , the set of transformers E t , and the set of regulators E r so that E = E l ∪ E t ∪ E r . Each edge, e ∈ E, is modeled as a series element using the Π-model representation to include the losses in a long distribution line as shown in Fig. 1. In this section we use the term line or edge interchangeably to refer to any type edges, that is, distribution lines, transformers, or regulators. The equations presented here are derived from first principles, and are applicable to the three cases. We denote v p n as the voltage phasor measured at bus n ∈ N and phase p ∈ Ω n , and i (n) nm,p denotes the current phasor flowing from bus n to bus m, n, m ∈ N on phase p ∈ Ω nm . Y nm ∈ C |Ω nm |×|Ω nm | denotes the admittance matrix of the line, and Y s nm ∈ C |Ω nm |×|Ω nm | is the charging admittance, modeled as a shunt element, half of which is assumed to be located on each end of the line. The vector of current phasors in the line i nm ∈ C |Ω nm | is related to the current flowing from bus n ∈ N to bus m ∈ N , i (n) nm ∈ C |Ω nm | , and shunt current i s nm ∈ C |Ω nm | by Kirchhoff's Current Law (KCL) as follows where the dependency between i nm , i s nm and the vector of voltage phasors v n ∈ C |Ω n | is given by Ohm's law as follows Combining (1)-(2), we obtain for every e = {n, m} ∈ E: and we can write the current flowing from m ∈ N to n ∈ N swapping the two indexes to obtain i (m) mn .
This is not true, however, when modelling transformers and regulators, i.e. Y (n) The derivations that yield our linearization are streamlined by introducing the following matrix, that we refer to as the matrix of power flows going from bus n ∈ N to bus m ∈ N : superscript H stands for complex conjugate transpose, or Hermitian. The following remark is in order: Remark 2: The actual power flows going from bus n ∈ N to bus m ∈ N are the diagonal elements of S (n) nm ; its off-diagonal terms in S (n) nm are redundant. To obtain our final equations we use a matrix operator D(.) such that D(A) is the column vector whose entries are the diagonal entries of A, i.e.: Hence, after several steps, we will use D(S (n) nm ) to obtain the vector of power flows.
In light of (3), the matrix S (n) nm is equal to: The non-linearity in (6) originates from the outer products v n v H n , v n v H m which are also matrices. In a nutshell, our goal in the next section is to linearize these outerproducts. Then, D(S (n) nm ) will bring about the desired linear relationship.

A. Linearizing the Products of the Voltage Phasors Using a Logarithmic Transform of the Voltage Phasor
Every complex number z can be expressed in polar form as z = re jθ where r is the modulus and θ is the angle with respect to the real axis. Similarly, we can express the modulus of the complex number as r = e log(r) . We denote u p n + jθ p n := log(v p n ) as the logarithmic transformation of the voltage phasor v p n on phase p ∈ Ω n , similar to the idea of [33] for the single phase case. We linearize the function that is of the form where | • | denotes voltage magnitude, and θ n is the phase angle at node n ∈ N . If we expand around the balanced case, i.e. |v n | = |v m | = 1 p.u. and θ n = θ m = 0 degrees. Then, the Taylor expansion is of the form and substituting in (8), we obtain the following expression Remark 3: Although we express the product of voltage phasors in polar form, one may also express these products using rectangular coordinates, and then expand the sin and cos functions individually. It should be noted that we use a first-order Taylor series expansion to avoid higher-order terms that would result in non-linear approximations.
The expression obtained in (9) directly applies Taylor to the voltage phasor. However, in this work, we use the logarithmic transform in a similar approach. Specifically we expand as follows which yields a similar expression. The logarithmic transform has inherent benefits when compared to other methods, and these are discussed below.
1) The Logarithmic Transform Preserves the Non-Linearity of the AC Power Flow Equations: Using the logarithmic transform is arguably easier than curve-fitting methods applied to the non-linear terms like the method presented in [24]. Although we take the logarithm of the voltage, i.e. u := log(|v|) to express the voltage as v := |v|e jθ , the non-linear is preserved.
2) Bilinear Terms in the (v, θ) Domain Become the Sum of Linear Terms in the (u, θ) Domain: By using the logarithmic transform, we are able to express the products as follows |v i ||v j |e j(θ i −θ j ) = e u i +u j +j(θ i −θ j ) which results in a linear expression in the exponent, c.f. (10).
3) The Logarithmic Transform Naturally Yields a Maclaurin Series Expansion: The exponent of the product of voltage phasors u i + u j + j(θ i − θ j ) goes to zero in the balanced case, which results in a Maclaurin series expansion. This is arguably easier than the intermediate steps in (8).  Relative error of the voltage phasor when using the logarithmic transform (in red) and the direct method (in black). In total, we draw more than 1000 samples. Expanding using logarithmic transform consistently yields better results.
However, a stronger argument for using the logarithmic transform is that it yields better results. To support this claim, we do a small study. Specifically, we want to understand whether using the approximation from (10) yields better results than (9). We draw samples of the voltage magnitude |v| ∈ [0.8, 1.2] p.u. and voltage angle θ ∈ [− 40,40] in degrees, and compare both methods to the the exact method in (7). The relative error of the methods with respect to the exact case in (7) is shown in Fig. 2. The results using the logarithmic transform consistently outperform those obtained from the direct method, which provides a stronger argument for its use. In the following subsection, we expand the use of the logarithmic transform to three-phase networks.

B. Using the Logarithmic Transform to Linearize the Three-Phase AC Power Flow Equations
With respect to [33] the main challenge is looking at the imbalanced case. For |Ω n | = 3, we can express v n as: v n := e u a n e jθ a n , e u b n e jθ b n , e u c n e jθ c n T (11) Given that the entries of v n are products of e u p n e jθ p n , ∀p ∈ {a, b, c}, ∀n ∈ N in the derivations that follow we use the operator diag(a) on vector a that returns a diagonal matrix whose entries are equal to the entries of a and express element by element products of two vectors a and b as the matrix vector product diag(a)b ≡ diag(b)a. Note that in the rest of the paper diag(.) is used not only to return a square matrix with the elements of the vector in the diagonal when the input is a vector, but it is also used to form a block diagonal matrix if the input is an array of matrices.
In unbalanced distribution systems, the phase angles θ p n , ∀p ∈ {a, b, c}, ∀n ∈ N by design tend to be within the neighborhood of the balanced caseθ n = [0, −2π/3, 2π/3], n ∈ N , and the voltage magnitudes are around 1 p.u. We leverage this domain knowledge to do a first-order expansion of the voltage phasors around 1 p.u. by expanding with respect to u p n ∀p ∈ {a, b, c}, ∀n ∈ N around zero. 1 We also re-center the phase angles prior to applying a first-order Taylor expansion around zero for the phase angles differences. Specifically, by adding and subtracting 120 degrees from phases b and c in v n . We obtain:

Remark 4:
It should be noted that Δ (3) n in (12) is a diagonal matrix that contains the constants as a result of re-centering the voltage angles around the balanced case. As a consequence, all the entries inθ n are close to zero, and not close to those inθ n . In addition, we make a distinction between Δ (3) n and Δ (3) m to accommodate lines that create a phase shift, such as those with transformers other than wye-grounded/wye-grounded. For all other lines, i.
We are now ready to use the expression for v n in (12) to perform our first order expansion for the outer products v n v H n and v n v H m . Recall that we have have defined u p n := log |v p n | where |v p n | is the magnitude of the voltage phasor, and e x ≈ 1 + x since log(1) = 0. Then, in matrix-vector form the first-order Taylor expansion around 0 with respect to u n is: where 1 denotes the all-ones vector. If we substitute (13) in (6), the outer products v n v H m are: m H Also, we know phase angle differences (θ n 1 T − 1θ T m ) are small, and we can further expand e jx ≈ 1 + jx. Dropping the second order terms, we obtain: where J is the all-ones matrix. Now we substitute (14) in (6), reorganize the terms and consider only the diagonal elements in S (n) is the operator that returns the diagonal values of a matrix. If we denote the vector of state variables as u nm [u n , u m ] T , andθ nm [θ n ,θ m ], after some tedious but straightforward algebra, we can express the linear power flow equation in matrix-vector form as follows: where: The matrixỸ nm is defined in the following expression: (17) where the symbol * denotes the complex conjugate and the matrix Y nm is known as the primitive admittance matrix, i.e.
Having expressed in (15) the power flow in the line in terms of the voltage on both sides of the line our next task, which is the subject of Section V is to use the expression as a stepping stone to directly relate the vector of power injections to the nodal voltages in the logarithmic representation. For that, we will leverage the incidence matrix of the network to go from a branch-flow to a bus-injection formulation. But first, to help express the equations in Section V in a more compact form, it is useful to stack all the vectors u e := [u n , u m ],θ e := [θ n ,θ m ], s e := [s nm , s mn ],ỹ e := [ỹ nm ,ỹ mn ], ∀e = {n, m} ∈ E as: and place Y u e , Yθ e , e ∈ E in a block diagonal matrix Using these definitions we can write all branch-flow equations for the lines at once as follows: thus, completing the derivation of our linear affine mapping. Remark 5: Note that the vectors u e andθ e will have the nodal variables repeated a number of times equal to the degree of the node in the graph. For example, in a simple 1 → 2 → 3 network, the voltage at bus 2 will appear twice. These redundant variables are necessary to balance the power flows in the node, so we convert the branch-flow equations into a network bus-injection model, by multiplying by the incidence matrix of the graph.
Before deriving a final expression that relates the vector of power injections and the nodal voltages through the system admittance matrix, i.e. Y bus , we introduce the models that account for the different characteristics of loads. These models, known as ZIP load models, consider the constant impedance, constant current and constant power behavior of a load. We present these models in matrix-vector form to incorporate them into the calculation of the solution.

III. MODELING MULTI-PHASE ZIP LOADS
In general, multi-phase loads can be modeled as phase-tophase in a delta connection or as phase-to-neutral in a wye connection as shown in Fig. 3. Let L denote the set of loads present in the network. We can define the set of wye-connected loads as constant impedance, constant current and constant power wyeconnected loads, respectively. Also, is the set of delta-connected loads. Each load ∈ L can be single phase, |Ω | = 1, two-phase |Ω | = 2 or three-phase |Ω | = 3, but it cannot have more phases than the bus it is connected to. We model loads ∈ L as elements that draw power s ∈ C |Ω | from the network given the voltage v ∈ C |Ω | . Loads can be constant impedance, constant current or constant power loads or a combination of the three. In this section, we introduce the models in matrix-vector form.

A. Modeling Wye-Connected Loads
Next, we provide linear expressions for the power drawn by a wye-connected load.
1) Constant Impedance (Z): Constant impedance loads draw power at different rates to maintain constant admittance of the load, y , ∈ L z Y . As shown in (6), we can calculate the power drawn by a constant-impedance load ∈ L z Y as follows where Y := diag(y ), Y ∈ C |Ω |×|Ω | , and the Taylor expansion is done on e 2u ≈ 1 + 2u .

2) Constant Current (I):
Constant current loads adjust the voltage angle and keep the current of the load, i , ∈ L i Y constant. We obtain an expression for constant-current loads ∈ L i Y as follows where ι := Δ (3) i * , I := diag(ι ) and Δ (3) contains the entries of Δ (3) corresponding to the load phases. The linearization follows from v := e u +jθ n ≈ Δ (3) n (1 + u n + jθ n ) 3) Constant Power (P): The modeling of constant power loads is trivial since the power flow equations are expressed in terms of power injections. Thus, we keep the power drawn by the load, s , ∈ L p Y constant. In the following subsection, we derive analogous expressions for delta-connected loads.

B. Modeling Delta-Connected Loads
Power is drawn by a constant-impedance s , ∈ L z Δ , constant-current s , ∈ L i Δ , constant-power s , ∈ L p Δ deltaconnected load. Delta-connected loads differ from wyeconnected loads in that phase-to-phase readings need to be converted to line-to-neutral. Loads draw line current i , ∈ Ω at a certain line-to-neutral voltage v . We may express the vector of line currents drawn by a delta-connected load i , ∈ L Δ as a function of the vector of phase-to-phase currents, i Δ , ∈ L Δ . When |Ω | = 3, i Δ = [i ab , i bc , i ca ] and we can use Kirchhoff's law to obtain the linear mapping between i and i Δ as follows  (25) is a linear approximation of the power drawn in the terminal by a delta-connected load. (25) assumes that voltage phasors are balanced, and that there is a 30°shift between wye and delta connections.
Expressions for delta-connected ZIP loads are presented next 1) Constant Impedance (Z): The power drawn by a constantimpedance load s , ∈ L z Δ is given as Furthermore, we define the constants: Substituting (14) in (26), taking only the diagonal elements, we obtain

2) Constant Current (I):
The power drawn by a constantcurrent load s , ∈ L i Δ can be obtained as follows

3) Constant Power (P):
The expression for the power drawn from the line can be obtained by using the transformation matrix, i.e. s ≈ Λ s Δ , ∈ L p Δ , completing the derivation of the linear expressions for ZIP loads.

C. Relating the ZIP Constants to the IEEE Test Cases
All the linear expressions from (22)-(29) use the ZIP constants, i.e. the wye-connected constants y , i , s , ∈ L Y and delta-connected constants y Δ , i Δ , s Δ , ∈ L Δ , to model the different characteristics of the load. In general, a load ∈ L is defined by its rated powers =p + jq in VAs and rated voltage readings |v | in Volts, as provided in all the IEEE test cases available in OpenDSS [35]. It should be noted that when |Ω | = 1, ∈ L Y , OpenDSS reports a line-to-neutral voltage. Every other time, the voltage is phase-to-phase. Thus, when we refer to |v |, we are referring to the value provided by OpenDSS in the definition of the load. Also, unbalanced loads must be defined as separate loads, since the voltage and power ratings will differ. Using these readings, we summarize the calculations of the ZIP constants in Table I where 1 |Ω | is the vector of ones with |Ω | entries. In the following section, we model different power delivery elements and provide algorithms to control some of these elements.

IV. MODELING AND CONTROL OF POWER DELIVERY ELEMENTS
In this section, we present the modeling of some standard equipment found in distribution systems. Modeling of highvoltage substation transformers, low voltage residential splitphase transformers, and voltage regulators as auto-transformers with control functionality are presented. Capacitor banks with switching capabilities are represented as negative loads that can inject reactive power. Solar inverters are modeled as negative loads that can control the amount of apparent power injection. We first introduce the modeling of the different components and then incorporate the controls into the models.
A. Modeling of Power Delivery Elements 1) Capacitor Banks: Capacitor banks are modelled as wyeconnected or delta-connected constant impedance loads that inject reactive power to the network. Let C denote the set of the capacitors on the network. We further divide C = C Y ∪ C Δ into wye-connected C Y and delta-connected C Δ capacitors. The impedance matrix of any capacitor is given by y c , c ∈ C Y in a wye connection or y c , c ∈ C Δ in a delta connection. The calculation follows from (22), (28) and Table I. The rated voltage of the capacitor |v c | is analogous to |v | in a load, and the rated power s c = −jq c is analogous to s , where q c is the reactive power of the capacitor bank.
2) Three-Phase Transformers: We model multi-phase multiwinding transformers as line elements with line admittances, such as those shown in (6), which enables the use of (1)- (18) to model transformers. In the Log(v) 3LPF linearization, transformers are represented in the network through the primitive admittance matrix Y nm as shown in (18). However, unlike in the case of a distribution line, the primitive admittance matrix of the transformer Y nm , {n, m} ∈ E t , {n, m} ∈ E l may not be symmetric. We follow the modeling from [40] to obtain the primitive matrix of the transformer. An example of a Wyegrounded Wye-grounded is shown in Fig. 4. In general, each transformer has a core, and its configuration is characterized by the connection of the transformer to both buses n, m ∈ N . If the transformer has more than two windings, e.g. a split-phase 120/240 V transformer, one may apply superposition and treat each winding separately. Let i sc nm ∈ C |Ω nm | denote the shortcircuit current measured on the transformer. The voltage drop v sc nm is related to i sc nm as follows where Y nm = Z −1 nm = diag([1/z 1 , . . . , 1/z |Ω nm | ] T ), and z i , i ∈ Ω nm , {n, m} ∈ E t is the short-circuit impedance in per-unit. Now, let ω nm ∈ C 2|Ω nm | denote the vector of currents in the core of the transformer {n, m} ∈ E t . When |Ω nm | = 1, ω nm = [ω p nm , ω p nm ] T , p ∈ Ω nm . The vector ω nm is related to the short-circuit current as ω nm = βi sc nm where β ∈ {−1, 0, 1} 2|Ω nm |×|Ω nm | is the incidence matrix that relates the two variables. We can also define w nm = [w n , w m ] ∈ C 4|Ω nm | as the vector of currents in the winding. When |Ω nm | = 1, w n = [w n , w n ] and w m = [w m , w m ]. Also, w nm = Nω nm , where N ∈ {0, 1, −1, √ 3, − √ 3} 4|Ω nm |×2|Ω nm | and the scaling factor √ 3 is used to go from line-to-line to line-to-neutral values.
The vector of currents on both terminals i nm = [i n , i m ] T can be expressed as i nm = Xι nm where X ∈ C 2|Ω nm |×4|Ω nm | is the matrix that maps winding current to terminal currents, and depends on the connection. It should be noted that Fig. 4 shows a Wye-grounded Wye-grounded transformer, but any configuration can be implemented through the matrix X. Finally, we multiply (30) by the incidence matrices from the left, and multiplying by their transposes from the right, we can find the relationship between the terminal voltages v nm = [v n , v m ] T and the terminal currents as follows thus, Y nm = XNβY nm β T N T X T , Y nm ∈ C 2|Ω nm |×2|Ω nm | is the primitive admittance matrix of the transformer, used in (17). Refer to [40] for an example on how to calculate the incidence matrices.

3) Distributed Energy Resources (DER):
We model DER resources, in particular solar resources and batteries, as constant power wye-connected loads as in (22). The IEEE 1547 [41] mentions Volt-Var and Volt-Watt as the general method to control the active and reactive power output of the solar panels. The power injection is a function of the voltage, thus, it behaves as a (negative) constant-power load. Also, the power available to the inverter is time dependent, and is subject to fluctuations due to the available power produced in the case of solar panels. This modeling assumption is not a hard constraint, and solar resources could be modeled as any load type in a delta or wye connection, using the constants presented in Table I. In the next section, we introduce solar inverters to control the output of the solar panel, as we will see in the next section.

B. Controlling Power Delivery Elements 1) Voltage Regulators:
We model a regulator r = {n, m} ∈ E r as a line element that connects two buses n, m ∈ N . Voltage regulators are modeled as single-phase auto-transformers, |Ω r | = 1, r ∈ E r . To represent the regulator controls, we use the modeling from [10], and include the controls Γ r in the calculation of the primitive admittance matrix Y nm as follows where τ i ∈ [−16, 16] ∩ Z is the tap number. The coefficients in (33) are due to the fact that most regulators allow voltage to be controlled within ±10% of its nominal value, in 32 steps. In general, the regulator controls the secondary voltage while tapping the voltage in the primary. Thus, γ i = 1, i ∈ {1, . . . , |Ω nm |}. To obtain the tap numbers, the regulator measures the voltage in the relay v r , given (34).
where κ is the ratio between the winding voltage and the regulator control voltage, ξ p is the rated primary current, in Amps, and z v = r v + jx v is the line drop compensator setting, in Volts. The tap number is given as follows where v reg is the regulator setpoint in Volts, v b is the bandwidth, and ν r is the voltage magnitude, v r = diag(ν r )e jθ r , where the values are given in Volts. In general, the voltage regulator setting is 120 V, thus, each step in the tap amounts to 0.75 V.
2) Capacitor Controls: Let s c denote the number of steps in a shunt capacitor, c ∈ C with switches η c,i ∈ {0, 1} that are to be controlled by our algorithm. The impedance values y c , y Δ c to be used in the modeling of a constant impedance load are given as We may consider different types of controls, such as current, voltage, kvar, power factor, or time. The control value in the capacitor η c,i is given as where ϑ is the control input, ϑ c , ϑ c are the maximum and minimum limits, respectively, of the interval outside of which the control operates, and u(.) is the step function. If a voltage override mode is used, then η c,i is given as (38) where v is the voltage measured at some point in the network, and v c , v c are the maximum and minimum voltage settings that allow for the operation of the control. Finally, the different types of controls that use ϑ in (37) may be obtained as follows where ξ r , ρ r are the ratio between the monitored line current and voltage and the corresponding current and voltage control, respectively. In general, the capacitor control measures the voltage v n where the capacitor is connected n ∈ N , or current and power values measured dowstream, i (n) nm , s n nm , {n, m} ∈ E. Also, the operators R(.), I(.) return the real and imaginary part of the phasor, and t denotes time.
3) Control of Smart Inverters: We model the control functionality of smart inverters using Volt-Var (VV) and Volt-Watt (VW) schemes as shown in [36]. The power injected to the system is a function of a piece-wise function, or "droop" curves. This control scheme injects reactive power when voltage decreases and transitions to VAR consumption when the voltage is high. The VV-VW curves are given as follows where η = [η 1 , . . . , η 5 ] are the control settings or parameters that characterize the VV-VW curves,p is the power that can be injected into the system. When modeling solar resources,p in (40) is in general the maximum power produced by the panel p = p unless curtailment is necessary for the safe operation of the system. The calculation of p is out of the scope of this paper, but can be obtained using methods like in [42]. Also,q is the available reactive power, andṽ is a low-pass filtered version of the voltage magnitude, introduced to avoid undesired controls caused by noisy measurements. The quantitiesq andṽ are given as follows.ṽ where τ m i is a time constant, s is the capacity of the inverter, and v i,t is the measured voltage magnitude. Finally, the power injected to the system s i,t = −p i,t − jq i,t , modelled as a negative constant power load, is calculated as follows where τ o i is a time constant. We conclude the modeling of distribution system elements and control functionality. In the next section, we express the power flow equations in terms of the power injections, by leveraging the incidence matrices of the graph. Expressing the equations in matrix-vector form enhances the study of efficient methods to compute updates, e.g. in the controls or in the loads, as we will see later.

V. BUILDING Y BUS AND SOLVING THE PF EQUATIONS
Let N := |N | n=1 |Ω n | and L = {n,m}∈E |Ω nm | denote the total number of nodes and lines in the network, respectively, and by E we denote the N × L incidence matrix of the grid. We can map the vector of power flows s to the vector of power injections s bus ∈ C N given by the incidence matrix, i.e. s bus = Es. Also, we can express the vectors of nodal state variables as u = E T u, u ∈ C N andθ = E Tθ ,θ ∈ C N . Multiplying from the left by E in (21) and substituting these expressions, we obtain where the full expression for the bus injections s bus can be obtained from the ZIP loads. Let us denote by L b A the total number of phases in the set of A−connected constant-b loads, A |Ω | and introduce the adjacency matrices E b A between the load phases and the set of bus-phases, for the set of A−connected constant-b loads: We stack in the vectors denoted by y A , i A and s A all the vectors of constants needed to define the A-connected ZIP load models.
and we can express s bus as follows: (43) includes all regulators when Γ r = I. For the case when Γ r = I, we can add the difference Γ r − I to (43). For that, we define E r ∈ [0, 1] N ×L r where L r = r∈E r |Ω r |, and Y u r , Yθ r are the entries of Y u , Yθ corresponding to the regulators. Including (46) and the action of the regulator controls in (43), we split the matrix into real and imaginary parts, and we obtain Log(V)3LPF :s bus =Ỹ busxbus (47) where:s with: and the matrices Y 1 , Y 2 , Y 3 , Y 4 are defined as (52) The solution to (47) (53) and we may recover the voltage phasors on each bus n ∈ N as follows

A. Discussing the Invertibility ofỸ bus
It is well known that a connected graph with distribution lines and loads is always invertible [43]. This can be proven by showing that the admittance matrices Y nm , {n, m} ∈ E are symmetric and positive definite. In this work, we assume that degree one nodes have loads attached to them, and that voltage regulators are modelled as tap changers installed on a wyegrounded/wye-grounded connected autotransformer, which result in a diagonal admittance matrix (symmetric) that is also positive definite. However, a matrix may become singular when delta-delta, wye-delta or delta-wye transformers are present in the network.
To accommodate these cases, we follow the procedure presented in [43], and add a small shunt admittance to the diagonal of the primitive admittance matrix. This method has also been exploited by tools such as OpenDSS. 2 As a result,Ỹ bus is an invertible but ill-conditioned matrix, dominated by the smallest eigenvectors. An ill-conditioned matrix may be regarded as an advantage, since the vectors u andθ tend to be concentrated in a subspace, and in such case low-order approximations may be used to recover these vectors, resulting in computational gains. However, this is out of the scope of this paper, and we leave it for future work. In the following section, we present an efficient way of updating the inverse ofỸ bus since updates tend to be sparse.

B. Using Sherman-Morrison-Woodbury for an Efficient Computation of Y −1 bus
In its traditional form, the Woodbury identity [44] says that the inverse of a rank-k correction of some matrix of the form (A + UBU T ) −1 can be computed by doing a rank-k correction to the inverse of the original matrix. Initially, this required the matrices A, B, U to be square and A, B non-singular. In our case, A is non-singular for a connected graph. Note that we include all regulators when Γ r = I to avoid disconnecting the graph. In our case, B may not be invertible, and the incidence matrix U is non-square. In [45], new expressions were derived that relax these constraints, namely Thus, the complexity of calculating the inverse in (55) lies in calculating the products and the inverse involving B, since it is the only element that is changing over time. The elements that are tied to the topology of the network A, U remain unchanged, and can be pre-computed.

C. Log(v) vs OpenDSS in Floating-Point Operations
It is known that commercial-grade software such as OpenDSS use LU decomposition methods [34] to calculate the inverse of Y bus ∈ R n×n , where n is also the rank ofỸ bus . In general, the computational complexity of the LU decomposition is 2 3 n 3 , and the inverse takes an additional 2n 2 floating point operations, where n 2 are for the forward substitution and n 2 for the backward substitution. In contrast, using the Woodbury identity to exploit the sparsity of the changes may speed up the process. Let the matrices in (55) be defined as A ∈ R n×n , B ∈ R k×k , U ∈ R n×k , where k is the rank of the correction in (55). The complexity of calculating the inverse using the Woodbury is partitioned as follows r Finally, the subtraction that completes the calculation of the inverse is O(n 2 ) And the complexity of each algorithm is given by the following expressions The performance of the Log(v) 3LPF changes with k which is given as To put these results in perspective, we show in Table II the floating-point  operations needed to compute the inverse of some of the IEEE   TABLE II  COMPUTATION OF THE MATRIX INVERSE IN FLOATING-POINT OPERATIONS and Low Voltage (LV) cases, ranging from the IEEE-13 to the IEEE-8500 to show different network sizes. It is shown that the Log(V) 3LPF scales with the size of the network and provides a significant advantage when computing the inverse, with the exception of the IEEE-34 in which there is an abundance of loads and regulators.
It should be noted that the LV test cases do not have capacitor banks, voltage regulators, constant-impedance or constantcurrent loads. Thus, only constant-power loads are included in these test cases. The Log(v) 3LPF is formulated in terms of power flows and voltages, hence, an update of the system admittance matrix is not necessary since all the changes in boundary conditions are implicitly captured in the vector of power injections. When n k, (56) reduces to Remark 7: In Table II, we use the big-O notation to provide an upper bound for the number of floating point operations required to compute the inverse of the system admittance matrix. However, it is outside the scope of this paper to write a low-level program to compare the computational time needed to parse and solve a test case. Thus, we do not guarantee that the implementation of the Log(v) 3LPF is the most efficient or that the right data structures have been used.

D. Inverting Y bus vs Using the Fixed-Point Iteration Method
It is often argued that, in order to solve the non-linear power flow equations, i.e. i bus = Y bus v, one needs to invert Y bus once. This is not true in general. Loads are modelled by its Norton equivalent that splits the non-linear loads into a linear part included in Y bus and a non-linear part modelled as a compensation current in i bus . Changes to the controls of the network, e.g. voltage regulators, capacitor controls or any switching device induce sparse changes in Y bus that will trigger an update of its inverse. Also, this way of modeling loads works when compensation currents are small [46], which may lead to convergence issues. In addition, we have shown in (57) that for n k, the inverse using the Log(v) 3LPF linearization is O(n 2 ) instead of a regular O(n 3 ). O(n 2 ) is in fact, an upper bound for the number of operations of a single iteration using the fixed-point method.  [21]. The results show the voltage magnitude measured on all phases. Some buses may have less than three phases.The ordering of the buses is based on the hop distance to the feeder head.

VI. SIMULATION CASE STUDY AND RESULTS
We present a study to showcase the numerical results that can be obtained by solving the Log(V) 3LPF equations, shown in (47). We test Log(V) 3LPF on different IEEE test cases, with different sizes, ranging from 41 to 8531 nodes, and different elements and configurations. These tests cases include ZIP loads, substation and split-phase transformers with 120/240 V loads, voltage regulators, capacitor banks, capacitor controls, switches and multi-phase lines.

A. Implementation of the Log(V) 3LPF Solver
We implement the Log(v) 3LPF equations in Python, and use the SciPy sparse linear algebra library to compute the matrix multiplications. SciPy is built on top of Numpy and uses the BLAS routines to perform linear algebra operations such as vector addition, scalar multiplication, dot products, etc. The Log(v) 3LPF is a wrapper built on top of the python package OpenDSSDirect 3 to parse the network model. The solver precomputes all the constant matrices in a fixed topology only once. Updates to the loads or the controls are implicitly modeled in the vector of power injections or as updates to the matrix inverse via the Woodbury identity, respectively. It should be noted that under very particular conditions, e.g. in the IEEE-34 test case, a regular Gauss-Jordan elimination may be more efficient than using the Woodbury identity, given k is large enough. During the initialization process the algorithm will calculate the complexity of calculating the inverse using the big-O notation, and will chose the optimal method. Also, in the LV test cases, the inverse is only computed once during the initialization, since all loads are constant power and there are no voltage regulators. Thus, all changes are modeled in the vector of power injections,s bus . The Log(v) 3LPF solver is available on Github. 4

B. Log(v) 3LPF vs Lossy DistFlow
We first compare the Log(V) 3LPF and Lossy DistFlow [21] linearizations with an OpenDSS solution on the IEEE-123 test case. It should be noted that the Lossy DistFlow linearization cannot handle constant-current loads due to the inherent limitations of the DistFlow equations. Also, the DistFlow linearization does not include the modelling on controls or devices other  than lines or loads. Therefore, we replace the elements for their line equivalent when possible, and treat constant-current loads as constant-power. These simplifications are for comparison purposes, and are not considered in any other analysis in this section. We assume a strong voltage at the feeder-head that is kept constant at 1.02 p.u. We compare the voltage profile on the three phases obtained using the Log(v) 3LPF, Lossy DistFlow and OpenDSS. The results from this comparison are shown in Fig. 5, with Root-Mean-Square Error (RMSE) values reported in Table III. The Log(v) 3LPF provides a better estimation of the voltage on the three phases, with an RMSE less or equal to 0.01 p.u. These results were expected since the DistFlow equations do not include phase angle information. In addition, the Log(v) 3LPF directly models the losses in the line, whereas the lossy Distflow requires the training of complex models to include the losses. This observation also applies to other linearizations such as the one shown in [25]. Nonetheless, the Log (V) 3LPF tends to underestimate the losses, because it tends to overestimate the voltage magnitude and because the phases tend to be around the balanced case. This trend will be seen in the upcoming analysis.

C. Log(v) 3LPF vs OpenDSS
In this subsection, we compare the Log(v) 3LPF with OpenDSS on all IEEE test cases from Table II including the modelling of all the elements available. In all cases, we assume a voltage source that is stiff enough, and that can be treated as an infinite source. All wye-delta or delta-wye transformers assume that the secondary terminal and downstream nodes lag the primary by 30 degrees. Switches, and in general voltage regulators, are modelled as lines with a very high admittance. In most cases, distribution lines are three phase and are modeled using their kron reduction. We first present in Table II the results pertaining the accuracy of the Log(v) 3LPF. The results are not included for the North American LV case due to issues processing the topology and interpreting the geometries. We show the RMSE and Mean Absolute Percent Error (MAPE) of the estimates for the voltage magnitude and angle that characterize the phasors.
In general, the errors are less than 3% of the voltages provided by OpenDSS which we take as the ground truth. We can provide some context for error incurred in the linearization presented here. It has been shown that the earth models used to calculate the impedance of a line, e.g. Carson's equations, may return values for the voltage that differ by more than 0.01 p.u. Therefore, these values are in line with the errors of the Log(v) 3LPF.
1) Solving the Power Flow in the IEEE-8500: The IEEE-8500 is arguably the most complex system available from the IEEE test cases, and perhaps the closest to a real distribution feeder. The IEEE-8500 test case contains multiple triplex lines, that connect the secondary of a center tapped transformer to 120/240 V loads, and are modeled as two-phase lines. In the IEEE-8500 case, loads are defined as two 120 V loads connected line-to-neutral. Also, center tapped transformers are defined as single-phase three-winding transformers and its line equivalent can be obtained by applying superposition in (31). These transformers are connected 7200 V line-to-neutral in one of the phases where the connection of the secondary depends on the polarity. To further showcase the Log(v) 3LPF solver, we present in Fig. 6 the voltage levels on phase A that correspond to 30 buses from the 4875 available in the test case. It should be noted that, as previously shown in Fig. 5, the Log(v) 3LPF tends to overestimate the voltage magnitude and angle. In terms of computational time, even for the IEEE 8500 bus case, our Python implementation took about the same amount of time of OpenDSS. However, we wish to remark that the comparison is not a fair one. In fact, we wrote the Log(v) 3LPF in an interpreted language (Python). One cannot draw conclusions from comparing it with a commercial-grade software like OpenDSS, which is designed to handle the data in an optimized manner, and benefits from being written in a compiled language (Pascal). The development of commercial-grade software is out of the scope of this paper.
2) Using the Log(v) 3LPF to Control the Network: We test the voltage regulator control on the IEEE-13 case. The IEEE-13 contains a single voltage regulator in the network, situated downstream of the substation transformer, keeping a strong voltage at the feeder head. Since there is a single control in the network, the convergence is smooth and controls do not counteract each other. In Fig. 7 we show the actuator applying the controls as well as the resulting impacts of these controls on the downstream bus.

VII. CONCLUSION
In this paper, we have presented Log(v) 3LPF, a linear power flow solver for unbalanced three-phase simulation on distribution systems. The benefit of the proposed linearization we highlighted is in the guarantees of having a solution and in allowing to use a sparse inversion method for the matrix inversion. The Log(v) 3LPF does not rely on any initialization, and it may suited to use under emergency conditions, e.g. a network reconfiguration that makes the network weakly meshed or other abnormal conditions. The Log(v) 3LPF includes the modeling and control of electrical equipment found on distribution systems. Our results indicate that the Log(v) 3LPF provides an accurate estimate of the voltage phasor, and that the error incurred in on par with other assumptions that are usually made to model distribution lines. Furthermore, we compare our results with the OpenDSS software and show that, in general, the Log(v) 3LPF can take advantage of the sparsity of changes in the boundary conditions to reduce the number of computations of the power flow solution. In the context of modeling and control of distribution systems, there are far broader applications of the Log(v) 3LPF that significantly expand its benefits beyond those of a reliable PF solver. The Log(v) 3LPF can be used for the analysis of controllers behavior in distribution systems, e.g. voltage regulator, switches, or capacitor controls. In particular, the inclusion of the phase angle as part of the solution in the Log(v) 3LPF linearization, unlike in the DistFlow equations, facilitates the study of the stability of distributed generators. Other potential applications include the use of the three-phase solver of a distribution OPF environment, by incorporating the power flow constraints as a linear operator, which is similar to using the DC assumption in transmission systems. Relatively recently, some attention has been given to transmission-distribution studies [47], that may use leverage the Log(v) 3LPF to express the power flow equations in both networks as a linear operator. Some of these applications will be further investigated in future research.
ACKNOWLEDGMENT Any opinions, findings, conclusions, or recommendations expressed in this material are those of the authors and do not necessarily reflect those of the sponsors of this work.