Newton Power Flow Methods for Unbalanced Three-Phase Distribution Networks

Two mismatch functions (power or current) and three coordinates (polar, Cartesian and complex form) result in six versions of the Newton–Raphson method for the solution of power flow problems. In this paper, five new versions of the Newton power flow method developed for single-phase problems in our previous paper are extended to three-phase power flow problems. Mathematical models of the load, load connection, transformer, and distributed generation (DG) are presented. A three-phase power flow formulation is described for both power and current mismatch functions. Extended versions of the Newton power flow method are compared with the backward-forward sweep-based algorithm. Furthermore, the convergence behavior for different loading conditions, R/X ratios, and load models, is investigated by numerical experiments on balanced and unbalanced distribution networks. On the basis of these experiments, we conclude that two versions using the current mismatch function in polar and Cartesian coordinates perform the best for both balanced and unbalanced distribution networks.


Introduction
The electrical power system is one of the most complex system types built by engineers [1]. Traditionally, electricity was generated by a small number of large bulk power plants that use coal, oil, or nuclear fission and was delivered to consumers through the power system in a one-way direction. Due to the modernization of the existing grid, a large number of new grid elements and functions including smart meters, smart appliances, renewable energy resources, and storage devices are being integrated into the grid. Thus, the existing electrical grid is changing rapidly and becoming more and more complex to control. A smart grid (SG) is offered as the solution to this problem [2][3][4].
In a smart grid, most of the new grid elements are directly connected to the distribution network which requires new types of operation and maintenance. The distribution network has been considered as a passive network that totally depends on the transmission network for control and regulation of system parameters. Conventionally, the power flow in the distribution system was one-way traffic (vertical) from the substation (only source) to the end of the feeders. However, the utilization of distributed generation (DG) made the distribution network active in the sense that the distribution network can generate electrical power in the network and transfer the extra power to the transmission network. This changes the direction of the power flow in networks into two-way traffic (horizontal). Therefore, central grid operators or transmission system operators (TSOs) of the power system must have different approaches for maintaining and operating the electrical grid because in this case, the main purpose of the operator has been adjusted to interconnect the various active distribution networks. As the distribution network becomes more active, there is an increasing role of distribution system operators (DSOs). For efficient operation and planning of the power system, it is essential to know the system steady state conditions for various load demands.
A power flow computation that determines the steady state behavior of the network is one of the most important tools for grid operators. The solution of power flow computation can be used to assess whether the power system can function properly for the given generation and consumption. Traditionally, power flow computations were calculated only in the transmission network and the distribution networks were aggregated as buses in the power system model. However, in the new operation and maintenance of the distribution network, the power flow problem computation must be done on the distribution network as well.
A reliable distribution power flow solution method will be required to solve a three-phase power flow problem in unbalanced distribution networks integrated with distributed generations and active resources (i.e., renewable power generations, storage devices, and electric vehicles etc.) [5,6]. There are conventional power flow solution techniques for transmission networks, such as Gauss-Seidel (GS), Newton power flow (NR), and fast decoupled load flow (FDLF) [7][8][9] which are widely used for power system operation, control and planning. However, these conventional power flow methods do not always converge when they are applied to the distribution power flow problem due to some special features of the distribution network: • Radial or weakly meshed (radial network with a few simple loops) structure: In general, a transmission network is operated in a meshed structure, whereas a distribution network is operated in a radial structure where there are no loops in the network and each bus is connected to the source via exactly one path.

•
High R/X ratio: Transmission lines of the distribution network have a wide range of resistance R and reactance X values. Therefore, R/X ratios in the distribution network are relatively high compared to the transmission network.

•
Multi-phase power flow and unbalanced loads: A single-phase representation is used for power flow analysis on transmission network which is assumed to be a balanced network. Unlike the transmission network, a distribution network must use a three-phase power flow analysis due to the unbalanced loads. • Distributed generations: Unlike conventional power plants connected to the transmission network, DGs have fluctuating power output that is difficult to predict and control since it is strongly dependent on weather conditions. Systems with the above features create ill-conditioned systems of nonlinear algebraic equations that cause numerical problems for the conventional methods [10][11][12]. Many methods have been developed on distribution power flow analysis and generally they can be divided into two main categories as:

•
Modification of conventional power flow solution methods : Methods in this category are generally a proper modification of existing methods such as GS, NR and FDLF. • Backward-forward sweep (BFS)-based algorithms : BFS-based algorithms generally take an advantage of the radial network topology. The method is an iterative process in which at each iteration two computational steps are performed, a forward and a backward sweep. The forward sweep is mainly the node voltage calculation and the backward sweep is the branch current or power, or the admittance summation.
In this paper, we focus on the Newton based power flow methods for balanced and unbalanced distribution networks with a general topology. Depending on the problem formulation (power or current mismatch) and specification of the coordinates (polar, Cartesian and complex form), the Newton-Raphson method can be applied in six different ways as a solution method for power flow problems. We refer to [65] for more details on all six versions of the Newton power flow method. In [65], the existing versions of the Newton power flow method [8,18,66] are compared with the newly developed/improved versions of the Newton power flow method (Cartesian power mismatch, complex power mismatch, polar current mismatch, Cartesian current mismatch, and complex current mismatch) for single-phase power flow problems in balanced transmission and distribution networks. It is concluded in [65] that the newly developed/improved versions have better performance than the existing versions of the Newton power flow method.
Therefore, we want to extend the Newton power flow methods developed for a single-phase problem in [65] to three-phase power flow problems. In this paper, only the polar current mismatch version is explained in detail for a three-phase power flow problem and the remaining versions can be derived similarly. Moreover, all six versions are implemented and compared with the BFS algorithm [43] for both balanced and unbalanced distribution networks. Different load models, loading conditions, and R/X ratios are considered in order to analyze the convergence ability of all extended versions. The key contribution of this work is new formulations of the Newton power flow method. Compared to existing versions of the Newton power flow method, our versions use different equations for PV buses in the Jacobian matrix that result in better convergence and robust performance. We present how these versions can be applied to unbalanced distribution networks by studying loads, three-phase load connections, three-phase transformers, and DGs. This paper is structured as follows. In Section 2, mathematical models of the power system, load, three-phase load connection, three-phase transformer, and DG are introduced. Section 3 mathematically describes the three-phase power flow problem. The Newton power flow method, the polar and the current mismatch formulations, and the polar current mismatch version are explained for the three-phase power flow problem in Section 4. The comparison result of all the versions of the Newton power flow method with BFS algorithm in balanced and unbalanced distribution networks is presented in Section 5. Finally, the conclusions are given in Section 6.

Power System Model
Power systems are modeled as a network of buses (nodes) and branches (transmission lines), whereas a network bus represents a system component such as a generator, load, and transmission substation, etc. There are three types of network buses such as a slack bus, a generator bus (PV bus) and a load bus (PQ bus). Each bus in the power network is fully described by the following four electrical quantities: For more details on the power system model we refer to [1].

Load Model
For load buses (PQ buses) in the network, active P and reactive Q power loads must be known in advance. In the power flow analysis, these loads (P and Q) can be modeled as a static or dynamic load. For the power flow computation, the static load models are used, so that active P and reactive Q powers are expressed as a function of the voltages. The following are commonly used models [67]: • Constant power (PQ): The powers (P and Q) are independent of variations in the voltage magnitude |V|: The powers (P and Q) vary directly with the voltage magnitude |V|: The powers (P and Q) vary with the square of the voltage magnitude |V|: The relation between powers (P and Q) and voltage magnitudes |V| is described by a polynomial equation: where a 0 , a 1 , a 2 and b 0 , b 1 , b 2 are constant parameters of the model and satisfy the following equations: The relation between powers (P and Q) and voltage magnitudes |V| is described by an exponential equation: where n is a constant parameter of the model.
Here P 0 , Q 0 , and V 0 are the specified parameters of the each bus in the network.

Load Connection
Three-phase loads can be connected in a grounded Wye (Y) configuration or an ungrounded delta (∆) configuration as shown in Figure 1. Loads are connected phase-to-neutral or phase-to-phase in a four-wire Wye configuration. Similarly, loads are connected phase-to-phase in a three-wire delta configuration. Let us assume that (P p i ) L and (Q p i ) L are the active and reactive power loads, respectively, at bus i for a given phase p and modeled as the exponential load model described in Section 2.1. Then, in the case the Y connection is applied, three-phase loads and currents are given as follows: In the case that the ∆ connection is considered, three-phase loads and currents are given as follows: (2)

Generator Model
Since conventional power plants have controls for the active power P and the voltage magnitude |V|, they are modeled as a PV bus in the power flow analysis. However, most of the DGs do not have both P and |V| controls and therefore they cannot be modeled as a PV bus. Figure 2 shows which type of power converter is employed to which types of renewable energy sources (DGs). Depending on the types of energy sources and energy converters, the DGs are modeled as follows: • The constant power factor model (PQ bus): The active power P output and power factor p f are specified and the reactive power Q is determined by these two variables.

•
The variable reactive power model (PQ bus): The active power P output is specified and the reactive power Q is determined by applying a predetermined polynomial function.

•
The constant voltage model (PV bus): The active power P output and voltage magnitude |V| are specified.
The DGs modeled as PQ buses can be treated as negative PQ loads in power flow analysis.

Transformer Model
Three-phase transformers are modeled by an admittance matrix Y abc T which depends upon the connection of the primary and secondary taps, and the leakage admittance.
where Y abc ps , Y abc sp are a mutual admittance and Y abc pp , Y abc ss are a self admittance of the primary and the secondary taps, respectively. The submatrices of the admittance matrix for different transformer connections are given in Table 2.
In this table, submatrices are given as: and y t is the leakage admittance of the transformer. If the transformer has an off-nominal tap ratio α:β where α and β are tappings on the primary and secondary sides respectively, then the submatrices must be modified as follows: • Divide the self admittance matrix of the primary by α 2 : The admittance matrix (3) for the transformer can be added to the general admittance matrix in (5). For more detailed information, we refer to [15].

Power Flow Problem
The power flow, or load flow problem is the problem of computing the voltage magnitude |V i | and angle δ i in each bus of a power system where the power generation and consumption are given. According to Kirchoff's Current Law (KCL), the relation between the injected current I and the bus voltages V, is described by the admittance matrix Y: In Equation (5), I abc i , V abc j , and Y abc ij are given as: where I p i is the injected current, V p i is the complex voltage at bus i for a given phase p, and Y pq ij is the element of the admittance matrix. The injected current I p i at bus i for a given phase p can be computed from Equation (5) as follows: The mathematical equations for the three-phase power flow problem are given by: where S p i is the injected complex power. Mathematically, the power flow problem is a nonlinear system of equations.

Newton Power Flow Solution Methods
The Newton based power flow methods use the Newton-Raphson method which is used to solve a nonlinear system of equations F( x) = 0. The linearized problem is constructed as the Jacobian matrix equation where J( x) is the square Jacobian matrix and ∆ x is the correction vector. The Jacobian matrix is obtained by ∂x k and is highly sparse in power flow applications [8,62]. Newton power flow methods (NRs) formulate F( x) as power or current-mismatch functions and designate the unknown bus voltages as the problem variables x.

The Power Mismatch Function
The power flow problem (8) is formulated as the power mismatch function F( x) as follows: where (S p i ) sp is the specified complex power at bus i for a given phase p. In general, the specified complex power (S

The Current Mismatch Function
The power flow problem (8) is formulated as the current-mismatch function F( x) as follows: where (S p i ) sp is the specified complex power at bus i for a given phase p. The power mismatch (10) and current mismatch (11) functions given in complex form can be reformulated into real equations and variables using polar and Cartesian coordinates. These two mismatch functions (power and current) and three coordinates (polar, Cartesian and complex form), result in six versions of the Newton-Raphson method for the solution of power flow problems. The detailed explanations of all six versions can be found in [65]. Only the version using the current mismatch functions in polar coordinates is explained in the following section. The remaining versions can be derived similarly.

Polar Current Mismatch Version (NR-c-pol)
In this version, the current mismatch function (11) is rewritten for real and imaginary parts using polar coordinates as follows: The current mismatch function can be written in vector form as follows: where and we look for the solution where the current mismatch function (14) is equal to zero: Application of a first-order Taylor approximation to the current mismatch function (14) results in a linear system of equations that is solved in every Newton iteration: where ∆ x is the correction and J( x) is the Jacobian matrix of the current mismatch function. Here, the Jacobian matrix is obtained by taking all first-order partial derivatives of the current mismatch function with respect to the voltage angles δ p and magnitudes |V| p as: The linear system of Equation (17) that is solved in every Newton iteration can be written in matrix form as follows: where The bus voltage correction in polar coordinate is given by: where h is the iteration counter. Then the complex voltage at bus i can be computed by:

Representation of PV Buses for NR-c-pol
In case of a PV bus j, a voltage magnitude |V abc j | is specified instead of reactive power Q abc j where In the current mismatch formulation, it is possible to choose the reactive power Q p j at the bus j for a given phase p as an unknown variable as a voltage magnitude |V| or an angle δ. Since Q abc j is an unknown variable, all the first-order partial derivatives ∂∆(I r i ) abc ∂Q abc j and ∂(∆I m i ) abc ∂Q abc j must be computed as: When the derivatives ∂∆(I r i ) abc ∂Q abc j and ∂∆(I m i ) abc ∂Q abc j are added into the Jacobian matrix J( x), the Jacobian matrix becomes a rectangular matrix. However, all derivatives of real ∆(I r i ) abc and reactive ∆(I m i ) abc current mismatch functions with respect to |V abc j | cannot be taken since |V abc j | is not an unknown variable. Therefore, we can remove all the ∂∆(I r i ) abc ∂|V abc j | and ∂∆(I m i ) abc ∂|V abc j | from the Jacobian matrix J( x) and the correction ∆|V abc j | can be replaced by ∆Q abc j . Thus, the Jacobian matrix J( x) is again square and therefore we will have a unique solution. The initial reactive power (Q p j ) 0 at PV bus j for a given phase p is computed as follows: The reactive power is updated at each iteration use.
The flow chart of the polar current mismatch version (NR-c-pol) is given in Figure 3. The remaining versions of the Newton power flow method (NR) such as Cartesian current mismatch (NR-c-car), complex current mismatch (NR-c-com), Cartesian power mismatch (NR-p-car), and complex power mismatch (NR-p-com) which are newly developed for a single-phase power flow problem in [65], can be extended to a three-phase power flow problem analogously. let h = 0 and x h = δ h | V| h the initial iterate compute the current-mismatches F( x h ) using (12) and (13) (18) Stop update iterate x h+1 using (21) and (22) yes no h = h + 1 Figure 3. Flow chart of the polar current mismatch version.

Numerical Experiment
We have shown how all versions of the Newton power flow method, originally developed for single-phase power flow problems in [65], can be extended to three-phase power flow problems with unbalanced distribution networks. Depending on the properties of a given network, one version can work better than the other versions. Therefore, it is crucial to study which version is more suitable for which kind of a power network. We use different load models, transformer connections, loading conditions, and R/X ratios in order to analyze the convergence ability and scalability of all versions. Different loading conditions are considered by multiplying each bus's power S by a constant k as S = k * S. Similarly, different R/X ratios are obtained by multiplying each branch resistance by a constant k as Z = k * R + ıX. Finally, the performance of the solution methods is evaluated for constant power and constant polynomial load models as defined in Section 2.1. The most widely used version using power mismatch function in polar coordinates (NR-p-pol [8]) of the Newton power flow method and backward-forward sweep-based algorithm (BFS [43]) are applied for comparison purposes.
Two balanced distribution networks (33-bus [70] and 69-bus [71]) and two unbalanced distribution networks (IEEE 13-bus [72] and IEEE 37-bus [72]) are used for the numerical experiment. All methods are implemented in Matlab and the relative convergence tolerance is set to 10 −5 . The maximum number of iterations is set to 50. Experiments are performed on an Intel computer i5-4690 3.5 GHz CPU with four cores and 64 Gb memory, running a Debian 64-bit Linux 8.7 distribution.

Single-Phase Problems
The convergence results of all solution methods for balanced distribution networks (DCase33 [70] and DCase69 [71]) are shown in Table 3. From Table 3, we see that NR-c-pol and NR-c-car versions have better performance in terms of a number of iterations and the norm of the residual of the mismatch function. Although NR-p-pol [8] and NR-p-car versions converged after the same number of iterations, the value of the residual norm is larger than for the NR-c-pol and NR-c-car versions. This means that if we set the tolerance to 10 −7 , these versions will need extra iterations to converge, whereas NR-c-pol and NR-c-car versions still converge after three iterations. We also see that NR-p-com, NR-c-com and BFS [43] methods need more iterations and have a linear convergence compared to other versions which have a quadratic convergence. These three methods solve the power flow problem in complex form, whereas other versions (NR-p-car, NR-p-pol, NR-c-car, and NR-c-pol) reformulate the problem into real equations using Cartesian and polar coordinates. Figure 4a shows the comparison of the results obtained by proposed solution methods with the well-known result of the existing method [48] for the computed voltage magnitude of DCase69. All the results of the proposed solution methods match the well-known result well with accuracy of 10 −5 as shown in Figure 4b.
A convergence result of all solution methods for the balanced distribution network DCase69 with different loading conditions and different R/X ratios, is shown in Figures 5 and 6, respectively. We see that NR-p-com, NR-c-com, and BFS [43] are more sensitive to the change of loads and R/X ratios compared to other versions that use real variables and values for the problem formulation using polar and Cartesian coordinates. Figure 7 shows the convergence of all solution methods for different load models. It is clear that for all methods, a constant power (PQ) model is more suitable to use on a balanced distribution network.
Thus, we can conclude that NR-c-pol and NR-c-car versions developed in [65] are more suitable for balanced distribution networks than versions using power mismatch functions (NR-p-pol [8], NR-p-car). Furthermore, NR-p-com and NR-c-com versions, as well as BFS [43] are the least preferable methods for balanced distribution networks in terms of convergence and robustness.

Three-Phase Problems
For IEEE 13-bus (DCase13) and IEEE 37-bus (DCase37) test networks, regulators are removed and all three-phase loads are chosen to be connected in a grounded Wye configuration as defined in Section 2.2. For the unbalanced distribution network DCase13, the transformer is connected in Wye-G, whereas DCase37 has the delta-delta transformer connection as defined in Section 2.4. The BFS method [43] is not implemented for three-phase power flow problems since it is not explained in sufficient detail how the three-phase transformer is handled for this method. The convergence result of all solution methods for unbalanced distribution networks (DCase13 and DCase37) is shown in Table 4.
From Table 4, we see that except for the NR-p-com and NR-c-com versions, all methods converged after the same number of iterations for both unbalanced distribution networks. However, NR-c-pol and NR-c-car versions have better performance in terms of both number of iterations and residual norm of the mismatch function, as we had the same result for balanced distribution networks. Again, NR-p-com and NR-c-com versions need more iterations to converge compared to other versions. A convergence result of all solution methods for the unbalanced distribution network DCase13 with different loading conditions and different R/X ratios is shown in Figures 8 and 9, respectively. As in single-phase cases, NR-p-com and NR-c-com versions are more sensitive to the change of loads and R/X ratios compared to other versions. Moreover, NR-c-pol and NR-c-car versions were more stable for the changes and therefore they can be applied to any unbalanced distribution networks with high R/X ratios and loading conditions. All methods performed better when a constant power (PQ) model was applied to three-phase problems as shown in Figure 10. We can conclude that versions using the current mismatch functions (NR-c-pol and NR-c-car) are more suitable than versions using the power-mismatch functions (NR-p-pol [8] and NR-p-car) for unbalanced distribution networks.

Conclusions
In this paper, the Newton power flow methods developed for single-phase problems in [65] are extended to three-phase power flow problems. Mathematical models of the load, three-phase load connection, and three-phase transformer connection are studied and applied in the numerical experiments. The power flow problem is mathematically described for the three-phase problem. The Newton power flow method and its six possible versions are introduced for three-phase power flow problems. As one of the newly developed versions of the Newton power flow method, the polar current mismatch version (NR-c-pol) is explained in detail for a three-phase power flow problem. The existing version (NR-p-pol [8]) and the backward-forward sweep-based algorithm (BFS [43]) are applied for comparison purposes. As a result of the numerical experiment, the polar current mismatch (NR-c-pol) and the Cartesian current mismatch (NR-c-car) versions developed in [65] and extended in this paper perform the best for both balanced and unbalanced distribution networks. We also investigate which version can be applied to what kind of a power network by comparing all versions for distribution networks with different loading conditions, R/X ratios, and load models. We observe that NR-c-pol and NR-c-car versions are more stable to the change of loading conditions and R/X ratios for both balanced and unbalanced networks, whereas the performance of other methods is highly sensitive to them. Therefore, these two versions are the fastest and the most robust methods of other versions that can be applied to single or three-phase power flow problems in any balanced or unbalanced networks. For the subsequent research, these newly developed versions (NR-c-pol and NR-c-car) will be applied to the optimal power flow problem in hybrid networks including DGs as they result in new formulations for equality constraints.