Matrix-Based Generalization for Power-Mismatch Newton-Raphson Load Flow Computations With Arbitrary Number of Phases

The standard power-mismatch Newton method is still frequently used for computing load flow due to its simplicity and generality. In this paper, a matrix-based generalization for the usual power flow equations to an arbitrary number of phases is derived. The proposed equations enable computing power injections and the Jacobian matrix in terms of submatrices that compose the network admittance matrix. Besides the more compact representation, another advantage of the proposed generalization is execution time reduction compared to the standard scalar formulation. Simulations are carried out to demonstrate the time reduction achieved via the proposed equations.

N , P Numbers of system buses and phases. Set of system phases. Hadamard product. ∂ x y Partial derivative of y with respect to x. V Vector of bus voltage magnitudes.

I. INTRODUCTION
Load flow algorithms are widely used in the analysis of power systems via simulation [1]- [5]. Due to the ample variety of possible applications and power system topologies, multiple methods have been proposed for solving load flow [6]- [18]. In general, it is acknowledged that Newton-type methods are the most general, since no restrictive assumptions regarding network topology are required. Furthermore, such methods are frequently favored due to their precision and quadratic convergence properties [9].
In what follows, we review representative algorithms that have been considered for solving load flow problems. The current-summation method [6]- [8] provides good performance when solving radial and highly-loaded networks, avoiding ill-conditioning problems that may be encountered when applying gradient-based methods. Its main limitation consists in not being applicable to meshed grids; considering that such topologies are increasingly common in the distribution scenario, this disadvantage is significant. Newton-type methods can be formulated either in terms of power-mismatch or current-mismatch functions [9]. The power-based approach corresponds to the traditional Newton-Raphson load flow, whereas the current-based formulation has been more recently proposed as the current injection method [10], [11]. It has been shown that the latter method tends to perform better than traditional textpowermismatch load flow due to the linearity of the currentmismatch equations with respect to voltage phasors. In [12], an additional contribution to current-mismatch load flow is given by incorporating an acceleration parameter to the current injection method, further improving convergence properties. In [13], the authors propose an approximation-based Jacobian decomposition for power-mismatch load flow which accelerates Jacobian computation. However, it is based applicable only to radial feeders with short lengths.
Among more recent methods, it is of special interest to comment those based on optimal power flow (OPF) and non-iterative approaches [14]- [18]. Many problems concerning power system operation require simultaneous optimization of system parameters (e.g. voltage and power losses) and other objectives such as, for instance, economical power dispatch and allocation of reactive support units. Hence, solving OPF problems is significantly harder than ordinary load flow due to the system parameters being simultaneously subjected to the load flow equations and objective-related constraints. The authors in [16] propose the usage of semi-definite programming for solving OPF and demonstrate that, under sufficient conditions, the method is always capable of obtaining a globally optimal solution. In [15], it is shown that relaxation of quadratically constrained quadratic programs can provide exact solutions in polynomial time for OPF in distribution systems. In [14], conditions under which relaxed optimal power flow achieves the global optimum for meshed networks are derived.
A non-iterative method for solving load flow, named holomorphic embedding load flow (HELM), was proposed in [18]. The main advantages yielded by this method are its independence from starting point selection and guaranteed convergence, given that at least one solution exists. Despite its attractive features, extensive comparisons with iterative algorithms carried out in [19] have highlighted its main disadvantages: relatively low precision and no guarantee of convergence to the correct solution (from a system operation viewpoint) in problems with multiple solutions. The deterministic nature of HELM makes it a strong candidate for real-time applications, in contrast to iterative methods that may pose a disadvantage to such applications. As an example, the authors in [17] have integrated HELM to genetic algorithm optimization in order to achieve load scheduling for minimizing voltage unbalance.
It must be noted that, despite the availability of all methods previously reviewed, Newton-Raphson power-mismatch load flow is still often favored due to its intuitiveness and convenient nature of direct specification of bus power injections [9]. Hence, any contribution that enhances its implementation and execution time is desirable. In this work, we present a generalized formulation of the power-mismatch equations that can reduce processing time for multi-phase systems.
To the best of our knowledge, no work in the literature has presented a matrix-form generalization to the powermismatch equations for an arbitrary number of phases. At most, scalar equations for power injection and Jacobian elements have been derived for three-phase systems [9]. This fact represents a further motivation for achieving a matrix representation of power-mismatch load flow, since it will permit abstracting from the number of phases when handling the power injection and Jacobian matrix computations. Although the usual single-phase and three-phase systems may be simple enough for using scalar notation, this is not the case with higher-order systems, whose practical implementation tends to expand due to the increasing power density being required for transmission [20].
The following are examples of high phase order power systems, whose representation in power-mismatch load flow is contemplated by the equations derived in this work: six and twelve-phase transmission lines [21]- [23], double-circuit three-phase transmission lines [24], [25] and pairs of pole-sharing three-phase distribution lines [26], [27].
In this paper, matrix equations are derived for powermismatch load flow, which can be seen as a direct generalization of single-phase equations where scalars are replaced by corresponding vectors and matrices. This result is of theoretical and practical interest because it: • Enables the processing of load flow directly in terms of bus mutual and self-admittance matrices, which can be previously computed from the specified network elements via their primitive admittance matrices; • Provides a concise, matrix-based representation of power and Jacobian equations for multi-phase systems, independently of the number of system phases; • Reduces execution time due to the multithreaded processing of matrix operations [28].

II. PRELIMINARIES
Consider a power system with N buses and P phases.
In each iteration of the power-mismatch Newton method, the values of voltage angles and magnitudes are iteratively updated according to the following equations: where the superscript n denotes iteration number, ∂ x y denotes the partial derivate of y with respect to x and the NPx1 vectors P and Q are the active and reactive power errors [29], respectively. The sensitivity matrices that compose the Jacobian in (2) are all NPxNP and given by H rs The Jacobian may be seen as being composed by PxP submatrices. The relation between such submatrices and the sensitivity submatrices is, for instance, where r, s ∈ . Identical considerations apply for M, N and L. The previously defined bus vectors can be interpreted in a similar manner, e.g.

III. DERIVATION OF POWER EQUATIONS
Consider an arbitrary P-phase circuit element connecting buses i and j. The relation between bus voltage and current phasors (denoted with a circumflex) is given by the following equations [30]:Î whereÎ ij is the Px1 vector of currents from bus i to bus j and the matrices Y The network admittance matrix Y is square of order NP. It can be expressed by means of PxP submatrices Y = [Y ij ] N xN , i, j = 1, 2, . . . , N , which can be computed via primitive matrices as [30]: (6) where i is the set of buses connected to bus i and Y (sh) ii is the primitive self-admittance matrix of bus i. Henceforth, the primitive matrices are considered as known via the specification of all circuit elements.
Theorem 1: Consider a power system with voltage magnitude and angle vectors V = [V i ] N x1 and θ = [θ i ] N x1 , respectively. The active power P i and reactive power Q i injection vectors for each power system bus i may be expressed by the following matrix equations: , r ∈ , the admittance submatrices Y ij = G ij + jB ij and the symbol denotes the Hadamard product.
Proof: Applying Kirchoff's current law to bus i yields:Î (9) whereÎ i is the vector of bus i nodal currents. Let S i be a column vector whose elements are the complex power injections in bus i. It can be computed as S i =V i Î * i , where (·) * denotes complex conjugation. Using (9) in this expression for complex power, we have: where Y ij = G ij + jB ij was used. The voltage phasor vector V i can be decomposed in polar coordinates asV , with r ∈ . Using this decomposition in (10), the following is obtained: VOLUME 8, 2020 (12) Recall that the vectors θ c It is straightforward to demonstrate that the following identities hold: Applying (13)-(15) to (11) and (12), the matrix-form active and reactive power equations given by (7) and (8) are obtained, as we wanted to demonstrate.

IV. DERIVATION OF JACOBIAN EQUATIONS
Consider that the Jacobian sensitivity matrices are composed by PxP submatrices, according to our previously proposed representation, e.g. H = [H ij ] N xN , i, j = 1, 2, . . . , N and likewise for the other sensitivity matrices. Furthermore, each submatrix can be decomposed in P column vectors of P-th order, where each vector is the derivative of a P i or Q i with respect to a given phase voltage magnitude or angle. In other words, all submatrices can be expressed as follows: where r ∈ . In Theorem 2, equations for the component vectors in terms of admittance submatrices and voltage magnitude and angle subvectors are presented and demonstrated.
where the index k is used to avoid confusion when i = j and 1 r is a P-element column vector whose elements are zero, expect for the r-th element which is equal to 1.
Proof: Consider at first the computation of H r ij . Assume i = j and apply the derivate with respect to θ r j to (7): where the sum over i disappears because the derivative is with respect to a single j ∈ i . In a similar manner for i = j, it follows that: Proceeding to M r ij , the following is obtained by deriving (7) with respect to V r j , i = j: For the case i = j, the following is obtained: The derivations for N r ij and L r ij are analogous, respectively, to those for H r ij and M r ij . Due to this reason, they are not detailed here.
Remark: In our computations, it has been assumed that the Jacobian matrix is of order 2NP and thus comprises all network buses, which are taken as PQ buses. However, this does not imply any loss of generality; for the cases of PV and V θ buses, the corresponding rows and columns of the Jacobian may be omitted to obtain fixed values of voltage magnitude and/or angle during all iterations.

V. PROCESSING TIME REDUCTION
Let O(S m ) be the asymptotic time complexity for solving the load flow problem S via the power-mismatch Newton-Raphson load flow, using one of the methods considered in this work (scalar or matrix formulation). That is, m can be scalar or matrix. It can be noted that the asymptotic complexity of the load flow problem is not altered when choosing the scalar or matrix processing method, since the matrix method is equivalent to applying the scalar method P times, and P is independent of N . For this reason, we can conclude that both methods exhibit, asymptotically (that is, N → ∞) and for every S, equal computational complexities, i.e. O(S scalar ) = O(S matrix ). Hence, the proposed method is identical to scalar power-mismatch Newton-Raphson in terms of convergence properties, scalability to larger systems and asymptotic complexity.
The objective of our proposed matrix formulation is to reduce execution time with respect to the scalar formulation, not asymptotic complexity. In what follows, we derive an expression for roughly estimating T matrix − T scalar and show that it may be expected to be negative, where T m denotes the average execution time of power-mismatch Newton-Raphson, using method m, for a given load flow problem.
Denote by t matrix and t scalar the average times required for computing any of (7), (8), (20)(21)(22)(23) via matrix and scalar formulations, respectively. Additionally, let t o be the average time required by the remaining computations (which is dominated by Jacobian construction and inversion) that are common to both methods. During a single iteration of method m, (7) and (8) are computed N times each, amounting to 2N computations. Furthermore, (20-23) yield a total of 4N 2 computations. Thus, a simple estimate of T m is: Using (28), we obtain T matrix − T scalar = 2 · (2N 2 + N ) · (t matrix − t scalar ). Denote by t p the average execution time required for computing a single scalar component of (7), (8), (20)(21)(22)(23). In the scalar formulation, t scalar = P · t p . However, due to multithreaded processing of the matrix formulation, it can be expected that t matrix < P · t p . For this reason, T matrix − T scalar < 0. It is of interest to highlight that this derivation suggests absolute time savings proportional to N 2 .

VI. PRELIMINARY VALIDATION -SMALL FEEDERS
Preliminary computer simulations were carried out in order to compare execution times of power-mismatch load flow achieved via standard scalar equations [9] and the previously derived matrix equations. The systems considered for testing were the IEEE reference 4-bus and 13-bus distribution feeders; in the 4-bus feeder, the wye-wye step-down transformer and unbalanced load configuration was selected. The 4-bus feeder is considered ideal for preliminary tests of three-phase modeling, whereas the 13-bus feeder is often used as a reference for power system studies and benchmarking of load flow methods. Such systems and their respective documentations are described in more detail in [31].
Both feeders were tested for different loading levels. To modulate system loading, a multiplicative factor α ∈ A = {0.5 + 0.1 · i; i = 1, 2, . . . , 10} is applied to all specified active and reactive power injections. The scalar and matrix-based load flow routines were implemented in MATLAB R2018a. For each combination of feeder, loading level and method (scalar or matrix equations), load flow was executed 1000 times in a dual-core 2.7GHz laptop. The selected convergence criterion was max i,r { P r i , Q r i } < , i = 1, 2, . . . , N and r ∈ , with = 10 −6 p.u. in the power system bases [31]. Execution time results obtained for the 4-bus and 13-bus feeders are given in Figs. 1(a) and 1(b), respectively. In Figs. 1(c) and 1(d), we plot the number of load flow iterations until convergence and the average time required for Jacobian inversion and auxiliary computations (t o , as defined in Section V).
The results show that, despite both methods converging in the same number of iterations for all loading levels, processing time was smaller for the matrix formulation. This corroborates what was discussed in Section V, in which it was claimed that the proposed method can reduce execution time even though it remains with the same asymptotic complexity as that of the scalar equations.
As could be reasonably expected, average execution time tends to be non-decreasing function of loading. Averaging over all values of α, percentage reductions in execution time were equal to 39.4% and 42.1% for the 4-bus and 13-bus feeders, respectively. The fact that average percentage time reduction was similar for both feeders can be explained by analyzing Figs. 1(c), 1(d) and (28). In Figs. 1(c) and 1(d), it can be seen that t o is relatively small with respect to execution time and, for this reason, does not dominate (28). Hence, we can apply the approximation t o ≈ 0 and obtain the following expression for evaluating relative time savings: which is approximately independent of N .

VII. VALIDATION CONSIDERING A LARGER SYSTEM
It has been pointed out in Section V that the scalability and quadratic convergence properties of Newton-Raphson load flow are mantained in the proposed method; additionally, the results from Section VI demonstrate favorable time reductions for the 4-bus and 13-bus test feeders via application of the derived equations. However, it may still be objected that time gains may not necessarily scale well for large N . In fact, the approximate analysis in (29) is based on the assumption of small t o , which is not guaranteed for large systems. Hence, in order to provide further validation of the proposed method, computer simulations were carried out for the largest of the original IEEE test feeders, namely the 123-bus system. Initially, experiments analogous to those of Section VI are carried out: for each α ∈ A, average execution times and t o are computed from 100 iterations of load flow, both for standard scalar equations and matrix formulation.
To extend the results of this case study, analysis of execution time as a function of system R X ratio is also carried out. To this end, a resistance factor γ is applied to the real part of the impedance matrix Z = R + jX of each distribution line of the 123-bus system, yielding modified matrices Z = γ R+jX which are then used to compute Y . Similarly to what was done for variable α, load flow is iterated 100 times for each γ in order to compute the averages of execution time and t o . In this particular case study, the adopted convergence criterion and γ domain are, respectively, = 10 −6 p.u. and γ ∈ = {1.0 + 0.1 · i; i = 0, 1, . . . , 30}. All the obtained results are given in Figs. 2(a) to 2(d).
The results suggest that time savings yielded by the matrix formulation scale well to larger N . In fact, the obtained average percentage reductions in execution time for variable α and γ were, respectively, 47.1% and 46.4%; such values are similar to those obtained in the previous case studies for the smaller feeders. It can also be observed that t o did not dominate execution time, as was the case for small N ; thus, the approximate analysis which yielded (29) still holds.
An additional point of interest is the fact that behaviors of the average time plots for scalar and matrix formulations are very similar, aside from the obvious difference in scale. Among the features of such similarity, it is worth noting that, for variable γ , both methods diverge for γ > 3.7, which was indicated in Figs. 2(c) and 2(d) with the notation T m → ∞.
The above-mentioned features of the obtained results suggest validity of the discussions in Section V regarding convergence properties of the matrix formulation. In fact, it was asserted that the proposed method reduces execution time of the Newton-Raphson method, despite having the same asymptotic complexity. This is precisely what the results indicate: the average time plots are reduced in scale, whereas behavior of execution time (including divergence)  as a function of system ill-conditioning (i.e. increasing α or γ ) are essentially the same for both scalar and matrix approaches.

VIII. CONCLUSION
In this paper, general equations in matrix form have been presented for power-mismatch Newton load flow with an arbitrary number of phases; a comprehensive derivation of such equations was carried out. The obtained equations allow the computation of injected power and the Jacobian matrix as functions of P-th order admittance submatrices and bus voltage magnitude/angle subvectors. Simulations involving the IEEE 4-bus, 13-bus and 123-bus test feeders were carried out and successfully confirmed time reduction via the proposed formulation. The results suggest that the proposed matrix generalization for the power-mismatch equations can be applied for reducing the processing time in applications of multi-phase load flow. The following characteristics and results could be observed related to the proposed approach: • Compared to ordinary scalar formulation, the obtained equations were shown to be of greater convenience notation-wise because the number of phases is abstracted from notation, which becomes more compact; • It was shown by means of theoretical considerations and computer experiments that usage of the derived equations yields execution time savings when compared to the scalar version of load flow. Combined with the fact that convergence properties are unchanged (since the Newton-Raphson method itself is not altered), this suggests that the obtained matrix notation should be preferred in order to reduce execution time; • The results suggest that relative execution time reduction is roughly independent of system size, which means that time gains due to the matrix equations have good scalability. Hence, the derived equations can in principle be recommended for systems of arbitrary size.
In future work, the authors intend to evaluate processing time benefits yielded by applying the derived equations in VOLUME 8, 2020 non-deterministic load flow problems, such as analysis of power systems with distributed generation. HENRIQUE P. CORRÊA received the bachelor's degree in electrical engineering from the Federal University of Goiás (UFG), Brazil, in 2018, where he is currently pursuing the master's degree in electrical and computer engineering. His main research topics of interest are power distribution networks, photovoltaic systems, and optimization applied to power systems.