An Iterative Scheme for the Power-Flow Analysis of Distribution Networks based on Decoupled Circuit Equivalents in the Phasor Domain

: This paper presents an alternative solution for the power-ﬂow analysis of power systems with distributed generation provided by heterogeneous sources. The proposed simulation approach relies on a suitable interpretation of the power network in terms of a nonlinear circuit in the phasor domain. The above circuit interpretation can be solved directly in the frequency-domain via the combination of a standard tool for circuit analysis with an iterative numerical scheme, providing directly the steady-state solution of the power-ﬂow of a generic distribution network. At each iteration, the resulting circuit turns out to be composed by two decoupled subnetworks, a large linear part and a set of smaller nonlinear pieces accounting for the load characteristics, with evident beneﬁts in terms of the computational time. The feasibility and strength of the proposed simulation scheme have been veriﬁed on a large benchmark consisting of the IEEE 8500-node test feeder. Then it is applied to the statistical simulation of a power network accounting for the variability effects of renewable generators. According to the results, the proposed tool provides an effective alternative to the state-of-the-art approaches for power-ﬂow analysis further highlighting the beneﬁts of the application of well-established tools for circuit analysis to power-ﬂow problems.


Introduction
From a very early time, the electrical power industry has faced many issues including generation-demand gap, reliability, planning, and operational back-offs and blackouts caused by weather and other external factors. In this framework, specific attention has been paid to the availability of modeling and simulation methods for both the transient and the steady-state assessment of the network state, also including a possibly large number of different heterogeneous distributed generator sources (DGs), nowadays spread into the power distribution networks. For the latter, several different analyses have been proposed in the literature including, power-flow, three phase-power flow, and harmonic analyses which have been proven to be mature tools successfully adopted in recent application problems [1][2][3][4]. Specifically, a probabilistic load flow with distributed wind generators is discussed in [1], a hybrid load flow involving multiple energy carriers is presented in [2], and the stochastic effects of the behavior of electric vehicles on the distribution networks is considered in [3].
Without loss of generality, efficient and reliable power-flow solution techniques have been developed and widely used for the assessment of power system operation [5][6][7][8]. Low resistive/reactive This paper is organized as follows. Section 2 describes the proposed numerical simulation of distribution networks, including additional details on the definition of loads. Section 3 collects the validation of the method and its application to standard test cases. The three-phase IEEE 8500-node test feeder is considered as the initial state-of-the art benchmark example. A second validation involving a smaller test case with renewable generators is then included and solved by both our custom MATLAB implementation of the equations and SPICE. For the latter example, the statistical analysis of the distribution network aimed at assessing the inclusion of the variability effects of distributed generators concludes the application results. Final remarks and conclusions are given in Section 4. Table 1 lists all the acronyms used in the text.

Power-Flow Analysis and Numerical Simulation
This section presents the simulation technique proposed in this paper for the power-flow analysis of a generic distribution network with renewables. The proposed solution relies on a circuit interpretation of the network along with its distributed generators in terms of a nonlinear circuit in the phasor domain, which is then solved via the standard tools for the circuit analysis.

Circuital Equivalent of a Distributed Network in the Phasor Domain
For the sake of illustration, the following discussion focuses on the example network shown in Figure 1a. It consists of a simple radial power distribution network fed by an ideal sinusoidal voltage source with phasorÊ(ω), in which the transmission lines are described by lumped RL equivalents represented by the gray blocks, whereas the circuital elements connected to the nodes 3 and 4 correspond to either an absorbing load or a renewable source providing a complex power S 3 = P 3 + jQ 3 and S 4 = P 4 + jQ 4 , respectively. Further details on possible alternative characterization of the network loads are provided at the end of this section.
The above distribution network can be interpreted via the equivalent circuit in the phasor domain depicted in Figure 1b. It consists of a nonlinear circuit in frequency-domain which provides a nonlinear map in the complex domain between a given set of isofrequential excitations at frequency ω and the corresponding phasors of the voltage and current responses at each node of the network.
To this aim, the loads of the network connected at node 3 and 4, characterized in terms of their adsorbed or produced power, have been readily converted into the corresponding nonlinear voltage controlled current sources f (V k ) defined through the complex power S k =V kÎ * k , yielding: where P k and Q k are the active and reactive power for a load at node k (e.g., k = 3 and 4). It is important to remark that the circuit schematic in Figure 1b can be easily described via the MNA [23,25], which represents one of the most established formulations for the solution of linear and nonlinear circuits, since it is currently used in a number of commercial simulation software such as SPICE. Specifically, according to the procedure in [25] the MNA allows analyzing the circuit in Figure 1b from topological inspection only, leading to the following matrix form: whereŴ ∈ C 5 is a complex vector collecting the nodal voltages and the currents flowing through current controlled elements, i.e.,Ŵ = [V 1 ,V 2 ,V 3 ,V 4 ,Î e ] T ,Â 0 ∈ C 5 is a complex vector collecting the independent sources of the circuit, i.e.,Â 0 = [0, 0, 0, 0,Ê] T and M is a complex matrix, readily derived from the topology of the circuit in Figure 1b by inspection, which writes: In the above equation Y i = 1/(R i + jωL i ) for i = 1, 2 and 3, where the dependency from the angular frequency ω is omitted for compactness. Furthermore, the complex vectorÂ(Ŵ) accounts for the nonlinear I-V characteristic of the network load and renewable elements:

Fixed Point Iteration Via Circuital Interpretation
Unfortunately, a trivial closed-form solution of the transcendental equation given by the MNA in (2) is not available, therefore such an equation must be solved numerically. Without loss of generality, the numerical analysis proposed hereafter focuses on a fixed-point iteration scheme.
The discussion starts considering the nonlinear circuit in Figure 1c, in which the linear and nonlinear portions of the network have been split into two sub-networks coupled through current and voltage controlled sources, respectively. It is important to remark that the above circuit turns out to be equivalent to our original network in Figure 1b. However, the proposed interpretation is the first step toward the so-called relaxation of the nonlinear problem.
The new interpretation of the nonlinear circuit in the phasor domain can be relaxed by replacing the voltage controlled current source in the linear network with an independent current sources of Figure 1d. This means that we are decoupling the linear network form the nonlinear one. According to this, the nonlinear MNA equation in (2) can be approximated by the following iteration scheme [26]:

, as shown in
, 0] T is a vector collecting the amplitude at the i-th iteration of the independent current sources connected to the linear portion of the proposed network interpretation (see the left part of the schematic in Figure 1d) andÂ (i+1) is a vector collecting the phasors of the currents flowing through the voltage controlled sources in the non-linear portion of the circuit (see the right part of the schematic in Figure 1d) defined by the voltages inŴ (i) and (4) as follows: Iterating the above equation leads to the following update rule: The above equation is iterated until a given convergence criterion is reached. A widely used idea is to use the L ∞ -norm error between the value of the unknown vectorŴ at current i-th iteration and the new prediction (i + 1)-th is less than a given threshold ε: The iterative rule in (7) derived from the relaxed circuit interpretation of the nonlinear MNA equation, turns out to be an alternative interpretation of the fixed-point iteration in the complex domain (see Appendix A for additional details).
The main advantages of the proposed formulation is that the linear portion of the system, which is usually large, needs to be solved only once during the iteration (i.e., the inverse matrixM −1 must be computed once since all the entries are known and do not vary). Also, different from the Newton-Rapshon iterative scheme usually used for the analysis of non-linear circuits in time-domain, the proposed technique does not require to compute the Jacobian matrix of the system, thus simplifying its implementation. The above simplification is motivated by the fact that, since we are working in the frequency-domain, the MNA equation turns out to be a static nonlinear equation in the complex domain (i.e., the equation does not contain any differential operator), thus leading to a simpler solution compared to more common dynamical cases. However, if it is needed, the iteration rules of the proposed approach based on the fixed-point scheme can be easily rearranged in order to cope with the more computationally expensive Netwon-Raphson method, being the latter a different interpretation of the proposed fixed-point technique. Moreover, the adopted circuital interpretation of a distribution network readily allows the circuit implementation in any SPICE solver (e.g., HSPICE, PSPICE, LTspice) with the solution in AC at a single frequency point. Additional details on the above feature will be provided in Section 3.2.

Source and Load Definition
Distributed generation sources (DGs) of the network can be described by small generators connected at various nodes in the distribution system to reduce the generation demand gap and improving both the reliability and voltage profiles across the network. They are connected in different configurations like power factor-controlled, voltage-controlled, or current-controlled modes. DGs can be traditional generators or renewable energy devices (photovoltaic, wind, fuel cells, etc.). (e.g., the blocks P 3 , Q 3 and P 4 , Q 4 in the schematic of Figure 1a). At any node, there are four parameters, real power (P), reactive power (Q), voltage magnitude (V) and voltage angle (θ). During the analysis, any two variables are known and two unknowns, thus DGs are classified into the type of known variables as PQ, PQ(V), PV (or PI). Since this work focuses only on renewables, the DGs term used hereafter refers to renewable energy sources, only.
Based on the information provided for the definition of the various types of DGs, each generator can be readily converted into a current source as briefly done for the simple PQ case in Section 2. For the PQ(V) case, current injection is still done via (1) because DG's real power and reactive power are specified, with the exception that the reactive power turns out to be defined as a function of the voltage magnitude. For the PV type, instead, the specified values of DG are the real power output and bus voltage magnitude. To calculate the equivalent current injections for this type, an iterative process is required. The mentioned procedure, further explained in [28,29], calculates the reactive power of the DG necessary to keep the bus voltage magnitude on the specified value. After the required power has been obtained, the equivalent current injection is computed as in (1) and integrated into (7) to calculate power-flow results.

Numerical Results
In this section, the effectiveness of the proposed approach for the power-flow analysis is investigated by considering two test cases. The first one consists of the benchmark IEEE 8500-node test feeder which represents a realistic full-size distribution system [30,31]. This first distribution network is aimed at validating the method and discussing its performance in terms of accuracy, efficiency, and convergence. Then, a second test case consisting of a single-phase network with 90 nodes is solved via the proposed approach for the power-flow analysis by considering both the MATLAB and its implementation in a commercial SPICE-based solver (i.e., LTspice from analog devices is used in this comparison). The obtained results are then validated via OpenDSS. The same network is then used for a practical application in which the DGs are spread along the network. For this second example, a statistical simulation for the systematic assessment of the possible effects of the variability of power rating of the DGs is carried out. For all the validations and applications, the tolerance ε in (8) is set to 10 −4 .

Validation: IEEE 8500-Node Test Feeder
The IEEE 8500 benchmark is considered for the performance assessment of the proposed simulation scheme. It is a large three-phase unbalanced radial network having both MV (medium voltage, primary distribution) and LV (low voltage, secondary distribution) levels and includes most of the components available in a typical North American medium voltage distribution feeder. The placement of the regulators and capacitors on the test feeder are shown in Figure 2. The three-phase network is analyzed through the proposed method detailed in Section 2, yielding to an unavoidable increase in the size of the involved matrices. Both the magnitude and the angle of the nodal voltages associated with the three lines, hereafter labeled as phases a,b and c, are computed. It is important to point out that transformers and regulators are modeled as in [20,32]. In addition, fixed regulator taps are considered as published in [24]. In order to provide a fair validation, the obtained results are compared with the ones obtained by means of the Z Bus method [20,33] and of the open source distribution system simulator (OpenDSS) [34], being the latter a reference tool released by the Electric Power Research Institute (EPRI). The voltage magnitudes and angles computed by the three methods are shown in Figure 3, highlighting identical results. In addition, Table 2 proposes a compact summary of the performance of the differed methods in terms of simulation time and number of iterations. The above numbers are generated running the scripts on a PC equipped with a CPU Intel Core i7 (seventh generation) and 32 GB ram. In addition, for the proposed approach, Figure 4 shows the maximum voltage difference computed at each iteration, highlighting the typical decreasing staircase behavior of fixed-point convergence where the horizontal line shows the convergence tolerance set to 10 −4 . The above figures and table highlight that the proposed approach has been proven to provide accurate results without any simulation overhead in terms of simulation time. It offers the unavoidable benefit of being general and flexible, thus not requiring any custom solutions for any arbitrary network.

Validation: 90-Node Single Phase Network
The test case considered in this section is a 10 kV single-phase radial distribution system with no transformers and regulators, hereafter termed as Test90. The topological structure of this test network is shown in Figure 5 and the parameters and data can be found in [35]. The network has been implemented within the simulation framework described in Section 2. For the validation, no DGs are connected and all the loads are PQ type and their corresponding current injections are calculated using Equation (1). Figure 6 shows the voltage magnitude curve solved through the propsosed MNA method validated with OpenDSS program. The same network is also solved through LTspice by replacing the MNA solution with a SPICE-based solution, leading to overlapped responses. The above comparison also stresses the flexibility of the proposed approach in allowing either a custom MATLAB implementation of the equations or the direct use of any classical software for circuit analysis. The circuital simulation is carried out in AC at the single operating frequency. The netlist is automatically generated based on the topology of the network and on the load definitions at each step of the iterative procedure and is feeding the LTspice solver which produces the phasor responses of the nodal voltages.  As an additional analysis, the considered Test90 power distribution is suitably modified as a weakly meshed network by closing the connections between the three dashed lines that represent tie branches (see Figure 5). The results of this test are shown in Figure 7. It is worth noting that the convergence for radial and weakly meshed networks are on the same order of magnitude of the first example, leading in this case to few dozens of milliseconds of CPU time (e.g., 16.9 ms for the radial configuration and 18.2 ms for the weakly meshed one). What is more important, the branches with a loop are treated in the same way, stressing again the generality of the proposed approach. Voltage profile of the Test90 network; differences between radial and weakly meshed configuration.

Application: DGs and Stochastic Analysis
In this section, a practical application is carried out for Test90 network where DGs are spread along the network. A statistical simulation aimed at assessing the effects of DGs (either in terms of location but specifically in terms of power generation). Table 3 provides additional details on the different configurations considered for this alternative application in which DGs are not included in the system (case I) or placed in different locations (cases II and III). The latter two cases include PQ and PV types of DGs, for which reactive power is corrected at each iteration according to [28,29] as stated in Section 2.3. Figure 8 collects the results of the simulation of the network for the above-mentioned configurations, thus highlighting the deterministic impact of the DGs, along with their placements. As a second and more realistic application, the Test90 radial distribution network is analyzed for variability analysis of DGs power ratings, to simulate practical conditions when DGs are spread along the network. Specifically, the power of DGs is assumed not being constant but is varying throughout the day according to weather conditions, with the aim of observing the change in the voltage profile. Two tests are carried out, based on Monte Carlo simulations performed by varying DGs power in the range [−20, 20%] around their nominal rated capacity (see Figure 9), and in the range [0, 100%] (see Figure 10). The latter case corresponds to an extreme simulation stressing the wide range of effects arising from the possible unpredictable behavior of DGs. In addition, Figure 11 shows the probability density function of the voltage of node 7 (i.e., the node with the highest load). From Figure 10 the band in the voltage profile is observed and the behavior of the distribution network (which relies largely on DGs) can be predicted. This can help utility and stakeholders to understand the effects of renewables' variability and unhealthiness of the network. It is to state that Monte Carlo simulation is used to observe the stochastic behavior of the network because we keep only two or three parameters of DGs as varying so simple Monte Carlo method holds, however, if the number of parameters is large, other methods are preferable for analyzing the uncertain behavior [36][37][38][39].

Summary and Conclusions
This paper addressed the compact modeling of a power distribution system using an alternative approach relying on a circuital interpretation of the network in terms of a nonlinear circuit in the phasor domain. A standard MNA formulation for circuit analysis is adopted and the resulting matrix equation is solved via a fixed-point iterative scheme, offering a simple solution in which the original network is split into a large linear part combined with and a sparse nonlinear portion. The validation and applications performed, which focused on two test cases, the so called Test90 network and the widely adopted benchmark IEEE 8500-node test feeder, proved that the approach offers good performance in terms of both simulation time and accuracy and turns out to be an excellent alternative candidate to be effectively used for large networks with multiple dispersed generators. A final application involving the statistical simulation of the Test90 network accounting for the uncertain effects of renewable generators is carried out.
It is important to point out that the fixed-point numerical scheme is chosen for demonstrating that even for realistic large networks, this choice offers a viable and efficient solution that does not require the computation of auxiliary matrices (such as the Jacobian matrix) at each step of the simulation procedure. The proposed method is, however, independent form the numerical scheme adopted and if needed, for problematic cases, other iteration schemes such as Newton-Rapson can be effectively adopted. In addition, the interpretation of the original network as the circuit of Figured possibly allows replacing the MATLAB MNA formulation of the network with a more intuitive circuital formulation which can be suitably implemented within any commercial SPICE-like solver. The above feature, which introduces additional user's flexibility and also allowed cross-validating the proposed technique, has been applied to the Test90 case. As a final summary, the key highlights and outcomes of this paper are:

1.
Verify the feasibility of using the waveform relaxation technique [26] to the power-flow analysis of power networks. This allows decoupling the inherent nonlinear network into the combination of a large linear subpart and a nonlinear portion, unavoidably enabling any tool for speeding up the simulation of the linear part; 2.
Demonstrate that the load flow analysis can be effectively carried out via standard tools for circuit analysis such as the MNA [25] (as done in other papers in literature as [23]), but also possibly replacing the simulation of the decoupled network via any SPICE-based commercial tool; 3.
Offer a flexible iterative procedure, which is combined with the simplest possible numerical scheme such as the fixed-point iteration. The mentioned procedure can be readily modified to incorporate other schemes such as Newton-Rapson; 4.
Verify that without working on the fine optimization of the routines, the implementation of the MNA-based formulae in MATLAB still produces the same accuracy and computational efficiency of reference commercial tools (e.g., OpenDSS) and of other custom implementations of state-of-the art load flow analyses (e.g., [20]).
formulation by recasting the MNA equation (2) in term of a fixed-point equation. Specifically, let us start reinterpreting the MNA equation in terms of the following nonlinear equation, whereF is a continuous vector-valued map which writes: Solving the above equation is equivalent to solving the following fixed-point equation, where:Φ (Ŵ) =Ŵ −F(Ŵ).
The vectorŴ * which solves the above equation is called fixed-point. The fixed-point of the transcendental equation can be estimated via a fixed-point iteration scheme starting from a first guess of the approximated solutionŴ (0) by iterating the following sequence, W (i+1) =Φ(Ŵ (i) ). (A5) The above iteration scheme converges after few iterations, if we can assume that the starting valuê W (0) is sufficiently closed to the actual fixed-pointŴ * . Substituting (A2) and (A4) in (A5) it can be easily proved that the above equation is equivalent to the one provided by the circuital-based iterative approach in (7).