A Three-Phase Power Flow Algorithm for Transmission Networks: A Hybrid Phase/Sequence Approach

In this paper, the three-phase generalization of a single-phase power flow (named PFPD) developed by the first author is presented. This three-phase formulation is chiefly conceived for HV/EHV transmission network applications, but it preserves a general validity for any power system. An iterative method for the solution achievement is throughout expounded. The algorithm quantitatively aims at investigating the impact of the asymmetrical transmission structures on power systems. This impact is evaluated in terms of voltage and current sequence components. Moreover, discussions on possible improvement actions to enhance the power quality are developed. The algorithm is implemented in Matlab environment and tested by several fictitious networks. Eventually, extensive comparisons in terms of execution time, number of iterations and solution accuracy with the software DIgSILENT PowerFactory are presented.


Y , Y N , Y SGL
Total bus admittance, network admittance, shunt admittance matrices (three-phase). Y G , Y L Generator and load admittance submatrices. Y GG , Y GL, Y LG , Y LL Admittance submatrices of Y (three-phase).

Y Geq , Z Geq
Admittance and impedance equivalent matrices as seen at the generator buses (three-phase

II. INTRODUCTION
The three-phase power flow problem is a widely investigated topic in the technical literature, e.g. [1]- [19]. The hypothesis of perfectly balanced three-phase network, in fact, cannot be achieved in real power systems. Nonetheless, the great number of publications is chiefly focused on the distribution networks (some considerable examples are [1]- [7]). This fact is mainly due to the unbalanced distribution load configurations, the typical low voltage four-wire distribution system structures, and the growing presence of distributed generation [1]- [4]. Moreover, the high r/x ratios of the distribution lines make the problem particularly challenging and worthy of research [5], [6].
For the transmission networks, instead, few contributions are presented in the last fifty years (after a careful review in the international technical literature, only the contributions [8]- [19] were found). In fact, it is often taken for granted that the transmission systems must be systematically operated in a balanced manner [1], [6], [8]. Therefore, power flows of transmission power systems can be computed by means of their equivalent single-circuit at the positive sequence.
Notwithstanding, an accurate knowledge of the voltage/ current unbalance factors for the HV and EHV networks is fundamental to make power quality evaluations [9], [20]. This topic is becoming more and more important for transmission networks and ought not to be underestimated. In fact, the network unbalance factors in the transmission power systems are going to increase [21], mainly for the following reasons: • Increase of the existing transmission line loading due to the growing electricity demand (also due to the recent difficulty of erecting new OHLs), • The possibility of incrementing the lengths of the transmission lines, • The transpositions of transmission lines are very rarely adopted [20], [22].
In this paper, a three-phase power flow algorithm (named as PFPD_3P) is presented. Beyond the evaluation of power quality, this algorithm can have a further practical/industrial fall out, by foreseeing the values of negative sequence currents in each section of the network. This is particularly important for synchronous generator/compensator negative sequence protections [21], [23]. Too high values of negative sequence currents, in fact, can bring to undesired generator and synchronous compensator protection trippings.
Twenty years ago, the first author had developed a matrix three-phase power flow algorithm (in the following named as PF [24]_3P) [24]. In the present paper, instead, the new three-phase power flow algorithm PFPD_3P is discussed and it takes its inspiration from a single-phase power flow algorithm (PFPD) recently developed [25]. In this way, a compact matrix approach can be exploited in order to achieve more rapidly the power flow solution. In fact, the interpretation of the (three-phase) slack generator inside an ''all-inclusive'' (three-phase) admittance matrix allows reducing the CPU time, decreasing the number of iterations, and achieving greater precisions compared to those of the other methods. This englobing is made possible by treating the slack generator as quasi-ideal current source at positive sequence. In PFPD_3P, it is demonstrated that such way of modelling the slack generator is also possible if the three-phase extension of PFPD is considered. Thus, this paper further demonstrates the possibility of treating the slack generator as a quasi-ideal current source, at positive sequence, in power flow problems.
With regard to the existing resolution approaches, the three-phase power flow algorithms are basically divided into two resolution categories [10]. The first ones only operate in the phase component frame of reference [5], [10]-15], so both symmetrical and asymmetrical components are modelled by means of their phase matrices. With the phase sequence approach, however, it is not possible to exploit the main advantage of the symmetrical component approach by means of Fortescue Transformations. The second ones only operate in the symmetrical component frame of reference so, in order to model the asymmetrical devices (like non-transposed lines), proper compensation techniques are adopted [9], [16]- [19]. These compensation techniques allow ''symmetrizing'' the asymmetrical components, so that the problem can be studied by using the three single-phase sequence circuits.
In the present algorithm, instead, each iterative cycle alternates the use of a sequence approach with that of phases by means of the Fortescue's transform and its inverse. As it is explained in the following, this is possible because some formulae involve the symmetry of some network devices (e.g., the formulae for computing the correcting currents into the synchronous generators). Whenever the Fortescue transform cannot be adopted (for instance, when non-transposed lines are involved), a phase approach is chosen. Therefore, differently from the existing three-phase power flow methods [8]- [19], PFPD_3P can be considered a hybrid one, by combining the merits of each frame of reference and without englobing the demerits deriving from the choice of a unique frame of reference. Thus, PFPD_3P is a novel three-phase matrix algorithm suitable for the study of transmission networks, characterized from a computational (hybrid) paradigm different from the existing ones.
As it is explained in the paper, the algorithm is easily self-implementable, since only five matrix iterated formulae are employed. Eventually, such three-phase evaluations can be carried out efficiently (as it is described in the section dedicated to the computational performances) without the need of using the classical numerical techniques (e.g., Newton-Raphson and derived).

III. THE THREE-PHASE POWER FLOW ALGORITHM BASED ON PFPD
Recently, the first author published an AC single-phase power flow algorithm (named PFPD) based on an ''all-inclusive'' bus admittance matrix. This means that all the power flow data (i.e., the network elements and the technical constraints) are embedded in a unique matrix, and an iterative procedure to achieve the solution is developed from it. Therefore, a concise, efficient, and rapid power flow has been presented [25].
The reasons that persuade the authors to investigate this topic is understanding if the three-phase extension of PFPD preserves its computational advantages. Obviously, the size of the three-phase power flow problem is greater than the corresponding single-phase one. Thus, any computational advantage is useful, even if the three-phase power flow is often intended as an off-line planning tool [12]. In light of this, a three-phase power flow formulation (PFPD_3P) inspired by PFPD is described in the following. In general, power flow solutions are determined by solving a set of equations formulated from some technical constraints (typically power and voltages) known a priori. In the present approach, all the three-phase network sections are divided into three sets: the slack, the generator, and the load sections. For these ones, the following technical constraints are set: 1. For the SLACK section: the three-phase positivesequence voltage u a,P = u a,P is constrained, and the phasor u a,P is assigned as the system angle reference; 2. For the Generator sections: the positive-sequence injected active power p b,P , . . . , p g,P and the positivesequence voltage magnitudes u b,P . . . u g,P are constrained; 3. For the Load sections: the constrained quantities are the complex power p h +jq h , . . . , p m +jq m absorbed when they are subjected to their positive sequence nominal voltage (even null, when transit sections are considered). It is worth noting that all these constraints are specified in the symmetrical component frame of reference. All these elements representing the constraints are included in the phase admittance matrix Y ABC SGL represented in Fig. 1. For the generic generator/load (i.e., shunt) element connected to the section i, in fact, the following relation can be written: where Y ABC i is the (3×3) admittance matrix linking the phase vector of the three entering currents in the shunt elements i ABC S,i (i.e., passive sign convention) with the phase vector of the three phase-to-ground voltages u ABC i . Obviously, for the active generators, (1) gives a set of negative currents (since they are injected from the generator terminals). By considering the total number m of the network   sections, a (3m)×(3m) square diagonal block matrix Y ABC SGL can be defined. This matrix allows writing the following phase relation including all the m elements connected to the m sections: where i ABC S = is the (3m)×1 block vector of the entering currents in all the section phases of the elements and u ABC = is the (3m)×1 block vector of all the phase-to-ground voltages. In a section i where no elements are connected (i.e., transit sections), the matrix Y ABC i is null. These situations, for instance, occur in the busbar sections connecting power transformers with electrical lines. Moreover, it must be reminded that also the slack generator is modelled and included inside Y ABC SGL . This inclusion is due to the formal possibility to consider the slack generator as a quasi-ideal current source, similarly to [25]. The three-phase admittance matrix Y ABC N , instead, links the phase currents entering the network with its phase-to-ground voltages: where the network system is of course supposed linear and i ABC N = is the (3m)×1 block vector of the phase currents entering the network and u ABC is the same of (2). It is worth reminding that the construction of Y ABC N is similar to the single-phase one, by using the Linear transformation techniques [12], [24].
By summing (2) and (3) member to member, the following relation is obtained: where Y ABC is the three-phase ''all-inclusive'' admittance matrix containing all the network information, u ABC is the set of the phase-to-ground voltages, whereas i ABC = is the vector of the net three-phase currents injected at the sections a ÷ m of Fig. 2.
It is worth noting that all the current sub-vectors are null except the first one (J ABC a ): it represents the three-phase external current injection due to the current source modelling the slack generator at positive sequence. As in [25], Y ABC allows obtaining a series of matrix relations describing the steady-state regime of the entire three-phase power system of Fig. 2. By introducing the matrix partitioning shown in Fig. 3, the following two sets of linear equations are obtained: and by applying the standard Gauss-Rutishauser matrix reduction techniques, (6) can be written as: and substituted in (5), it follows: The phase matrix Y ABC Geq models the Ward equivalent network [33] as seen at the generator sections. By observing the reversed structure of (8), i.e.,: it can be seen that the vector of the injected currents by the slack current generator can be computed as follows: where Y ABC Geq,aa is the inverse of the diagonal first block matrix of the impedance matrix Z ABC Geq (see Fig. 4): it represents the three-phase concept of admittance as seen at the slack bus section.
It is worth noting that the above-mentioned mathematical steps can be thought as the three-phase generalization of the steady-state formulation of PFPD [25].
Notwithstanding, differently from PFPD, the current vector J ABC a of the slack bus is not immediately known by means of (10). In fact, only the positive sequence slack voltage u a,P is scheduled, whereas u ABC a is not completely known a priori, as it depends on the negative and zero sequence voltages of the slack section, i.e., In turn, the negative and zero sequence slack generator voltages depend on the negative and zero current circulations in its sequence networks, as shown in Fig. 5.
However, the slack section sequence voltages are not analytically determinable. In fact, if the network as seen at the slack section is assumed to be unbalanced, the extent of this network unbalance varies from case to case. Thus, an iterative scheme to determine such sequence voltages allows finding u ABC a . This iterative scheme is described thoroughly in Sect. IV. Once the convergence is achieved, the knowledge of u a,P , u a,N ,u a,0 allows computing u ABC a by means of (11). As a result, J ABC a can be computed by means of (10), so the vector i ABC is immediately known. The subsequent application of (9) and (7) allows computing u ABC G and u ABC L respectively. Therefore, the three-phase steady state regime of the network is completely defined.

IV. THE PHASE COMPONENT MODELLING OF THE NETWORK COMPONENTS
In this section, it is explained how both the power flow constraints and the network elements are modelled and included inside the three-phase matrices Y ABC SGL and Y ABC N . The Y ABC SGL matrix holds all the generator and load models of the network.

1) SLACK GENERATOR
One of the main novelties of this paper deals with the slack generator treated as a quasi-ideal current generator. This choice comes from the possibility to make a source transformation by considering the impedance of the equivalent positive sequence network exciting the system as infinitesimal [25] (e.g., j10 −5 p.u.). Hence, the slack generator positive sequence network can be considered as an ideal current source in parallel with an infinite shunt admittance (see Fig. 5). However, it was demonstrated in [25] that from an engineering point of view, an infinite admittance is an admittance with a very arbitrarily large value (e.g., -j10 5 p.u., as it will be detailed in Sect. IV).
Thus, the sequence component admittance matrix Y 0PN a modelling the three-phase slack generator can be built (see Therefore, the three-phase slack generator modelled as an admittance matrix can be embedded inside the ''all-inclusive'' Y ABC SGL matrix.

2) GENERATORS
Synchronous machines are the most employed devices for power generation in the transmission networks. From a structural point of view, a synchronous machine can be considered a symmetrical device, so its steady-state regime can be studied by means of the sequence networks represented in Fig. 6: only the positive sequence network is the active one, since the electromechanical power conversion takes place in this sequence. The presence of the negative and zero sequence passive networks, however, is fundamental to consider possible voltage distortions, due to negative and zero-sequence current flows. For steady-state transmission networks, these current flows are caused by the structure asymmetries (i.e., mainly non-transposed lines), and by the unbalanced loads. By considering a generator connected to the section g, its equivalent positive sequence admittance can be computed as in the following relations: Eq. (13) is the three-phase generalization of the PV constraints modelled in PFPD [25] as negative (sign, not sequence!) passive admittances: this generator modelling allows their inclusion into the unique Y ABC . Therefore, the building of the diagonal sequence admittance matrix Y 0PN g , shown in Fig. 6, is immediate. The phase matrix Y ABC g can be immediately found by means of the Fortescue transformations applied to Y 0PN g : Eventually, the phase admittance matrix Y ABC g can be stored in the g-th position of Y ABC SGL (Fig. 6).

3) LOADS
A three-phase balanced load in section m can be easily treated by means of the sequence approach [24]. For such loads, the complex power absorbed by each section terminal are equal to: As shown in Fig. 7, these loads can be thought as composed of three different types: (1) the asynchronous load, (2) the static, and (3) the reactive power compensation (to achieve a specific power factor target) devices. Such kinds of loads can be considered in parallel with each other and no mutual coupling is considered.
By denoting with S t the complex power absorbed by each of the three types of load under the nominal positive voltage (1 p.u.), its positive sequence star admittance can be computed as: Obviously, the positive complex power S P,t absorbed by each load typology must be known a priori.
Eq. (14) is the generalization of PQ constraints modelled in PFPD [25] as passive admittances: this load modelling allows computing its corresponding phase frame of reference matrix to be included into the unique Y ABC total bus admittance matrix (see (15)).
For the static load and the reactive power compensation devices, the positive and negative sequence admittances are equal. In the case of the asynchronous loads, instead, the negative sequence admittance can be considered equal to ξ y P,1 e jψ , where ξ = 5÷7 and ψ = −60 • ÷ −75 • [24]. Moreover, since each type of load is considered as being symmetrical, its corresponding sequence admittance matrix is diagonal (Fig. 7).
Eventually, by exploiting the linearity of the system, these sequence matrices can be summed, and the phase admittance matrix Y ABC m of the balanced load in the section m can be computed by means of the Fortescue transformations: Similarly to the generator shown in Fig. 6, Y ABC m can be stored in the m-th position of Y ABC SGL . Moreover, the conciseness of this matrix approach allows modelling all the radial subtransmission/distribution networks supplied by the EHV/HV transmission sections as single equivalent loads (see Fig. 8).
For instance, the equivalent load formation of the radial network of Fig. 8 is described. This network consists of ten sections, so the associated three-phase network admittance matrix is (3·10)×(3·10), and it must be summed to the shunt admittance matrix (in which the loads are stored in the 8-th, 9-th and 10-th positions).
Consequently the (3·10)×(3·10) ''all-inclusive'' matrix Y ABC R can be computed. Therefore, the impedance matrix It can be immediately noted (see Fig. 9) that the sub-vector of the absorbed currents at the section 1 can be computed as follows: where admittance matrix modelling the entire network as an ''equivalent load'' as seen at the section 1.
holds all the other elements of the network but the generators, and the loads. The steady-state three-phase modelling of the network devices and their insertion inside Y ABC N is a well-investigated topic [12], [24], and different modelling approaches can be adopted.
The modelling used in this paper is briefly presented in the following. Transmission OHLs, insulated cables, and GILs are firstly modelled by considering the presence of all their active and passive conductors (e.g., ground wires, metallic screens/armours, and enclosures). For this purpose, a multiconductor approach with n conductors in parallel with each other and with the ground is exploited [27]- [29]. Thus, single and double-circuit lines, and other complex configurations can be indifferently considered, without neglecting the real asymmetries of each line. The method is known as Multiconductor Cell Analysis (MCA) developed in 2009 by R. Benato. In order to provide some brief reminders, in Fig. 10, an elementary cell of a generic line is shown: it is modelled by means of three blocks: the longitudinal one (L), and the two transversal ones (T S and T R ). These blocks form a generalised multiconductor π-circuit, and are modelled through the admittance matrices Y L , Y TS , and Y TR . For electrical lines, it is almost always Y TR = Y TS = Y T . It is worth reminding that Y L , Y T are computed by means of the Carson/Schelkunoff-Pollaczek/Carson-Clem formulations [30]- [32], and the classical formulations for the capacitive and conductive couplings between all the line conductors respectively.
By combining Y L , Y T , the matrix formulation linking the entering currents at the sending and receiving ends of the cell with their voltages is represented in (18).
Then, the matrix cascade of each cell for computing the matrix of the entire line is applied. Fig. 11 shows the analytic formulation of the equivalent matrix of two cells: this operation must be repeated depending on the number of cells (if c is the number of cells, c-1 computations must be done). Thus, a unique (n × n) admittance matrix of the entire line as seen at its terminals can be computed.
It is important to highlight that this MCA allows representing each element inside each cell. All these elements must be modelled by means of suitable circuits: their admittance matrix representations are performed by considering electrotechnical evaluations [27], [28], [33], [34].
Once the (n × n) admittance matrix of the electrical line is computed, the Kron's matrix reduction [35] is applied in order to consider the behaviour of the system as seen only at the phase conductors (see Fig. 12).
This reduction is fundamental and unavoidable since this power flow involves only the three phase active conductors.
This reduction can be achieved under the assumption that either the voltages or alternatively the currents of the passive  . Matrix equivalent to the cascade of two multiconductor cells [31].
conductors at the ends of the electric line are null. Hence, the dimensions of the line admittance matrices decrease (e.g., (6 × 6)-matrices for single-circuit lines). Thus, the effects of the passive conductors on the phase conductors is not neglected, even if this assumption of null voltages of the passive conductors is not completely true. For instance, Fig. 12 shows this matrix reduction for a five-conductor OHL line (three phases and two ground wires) for two different ground wire arrangements: 1) the two ground wires are earthed at their ends, 2) the two ground wires are unearthed or insulated at their ends. In both cases, the null voltage/current components allow erasing the columns of the Z line /Y line matrices associated to the passive conductors. Similarly, the rows of the Z line /Y line matrices associated to the passive conductors can be erased also, since both the voltage and current quantities of the passive conductors must not be considered. Therefore, the (6×6) admittance matrices Y line,eq can be computed for both the cases. In Sect. VI, the influence on power losses of the two different ground wire practices are thoroughly evaluated. In this paper, this procedure is named as row/column elimination technique. It is worth noting that for the formation of Y line,eq in the first case (earthed ground wires), no matrix inversion is needed. This fact is relevant, since the great majority of the OHL ground wires in real power systems are earthed (at the substation earthing grids and also at each tower). It can be proved that the computation of Y line,eq can be obtained also by means of the standard Gauss-Rutishauser matrix reduction techniques. However, this procedure is more computational expensive than the row/column elimination one. This matrix reduction is also applicable to insulated cables, and to GILs. However, the following considerations must be kept in mind when the power flow solutions are assessed: • The Kron's reduction on OHL is generally based on a ''light'' assumption. In fact, the passive conductors of such lines are earthed at each tower and at the substation earthing grids, so the assumption of null voltage is quite accurate.
• The Kron's reduction on insulated cables is generally based on a ''heavy'' assumption. In general, in fact, the induced voltages in the passive conductors could not be negligible. In particular, they depend on the different screen bonding/earthing techniques (e.g., cross bonding arrangement for HV/EHV cables [36], single-point bonding for very short HV/EHV cables, and solid bonding arrangements for MV/LV cables). Differently for GILs, the solid-bonding arrangement (or even the multiple point-bonding 37], [38]) allows considering the enclosure voltages practically null.
In any case, these simplifying hypotheses are considered licit for a good three-phase power flow study of a transmission network, as shown in Sect.V. A very precise steady-state regime evaluations on the passive conductors can be achieved by means of MCA [27], [28], [33], [39] which is a powerful circuital tool, but not a power flow one. Further researches are ongoing to combine MCA with this power flow so to obtain a general multiconductor power flow. It is worth reiterating that the steady state regimes of passive conductors cannot be intrinsically known with the present power flow.

2) POWER TRANSFORMERS
As the synchronous machines, the two/three-winding transformers can be considered symmetrical devices from a structural point of view. Thus, their steady-state regime can be licitly studied by means of their sequence networks. For a two-winding transformer, its positive-sequence two-port network is immediately inferable from its nameplate data. Its negative-sequence two-port network is the same of the positive one, except for the phase shift angle which is always the opposite of the positive one. The zero-sequence two-port network, instead, depends on its particular winding earthing and configuration.
For these networks, their (2 × 2) admittance matrices can be computed, and opportunely put together in a unique (6×6) admittance sequence matrix Y 0PN 2w−tr . Then, the phase admittance matrix Y ABC 2w−tr can be computed by means of the Fortescue generalised transformation [24]: where F is the generalized (6 × 6) Fortescue matrix. The aforementioned procedure is also applicable to the three-winding transformers. The three-port sequence networks are described by means of (3 × 3) admittance matrices. These matrices are put together in a unique (9×9) admittance sequence matrix Y 0PN 3w−tr , and then the phase matrix Y ABC 3w−tr is computed.

3) SHUNT ELEMENTS
The following power system shunt elements are considered: • capacitive/inductive shunt compensation in electrical substations; • shunt compensation of long EHV/HV cable systems. Their admittance matrices can be easily computed by means of the inspection method or by inversion of VOLUME 9, 2021 their impedance matrices [28] and must be overlapped in the suitable positions of the three-phase admittance matrix Y ABC N [36].

V. THE HYBRID ITERATIVE PROCEDURE
The complete knowledge of Y ABC N and Y ABC SGL (and thus of the ''all-inclusive'' matrix Y ABC ) allows computing the three-phase power flow solution by means of the procedure described in Sect. II. However, the elements of Y ABC SGL are not completely known since they depend on the power-flow solution.
Thus, an iterative method to compute the three-phase power flow problem is proposed. The iterative method does not need any numerical analysis technique, as in [25].
It is well known that a good initial guess is necessary for any iterative algorithm. In particular, the initial guess of Y ABC SGL is built by considering the following values: For the Slack Section: y a,P = −j · 10 5 p.u.
As in [25], the slack generator is conceived as a quasiideal current source: hence, the slack section constraint can be stored inside Y ABC SGL . Since the scheduled slack voltage is the positive one, the quasi-ideal current generator is set in the positive sequence network, as shown in Fig. 5.
The discussion about the choice of the initial value of y a,P is the same of the one described in [25].
For the negative and zero sequence networks, instead, the typical synchronous machine admittance values are considered, according to the machine typology chosen for the slack one.
For the Load Sections: The load matrices Y ABC m of all the load typologies described in sect. IV are initially computed by considering their positive sequence nominal voltage, i.e., 1 p.u.
The three-phase voltages in the load sections, in fact, are actually the unknowns of the problem. However, if the network interconnections do not introduce high voltage drops, the choice to keep the load voltages to their nominal positive sequence value (1 p.u.) can be considered as a valid initial guess.
For the Generator Sections: Eq. (13) allows computing the first guess for the positive sequence admittance which models the generator in section g: Eq. (20) is the application of (13) for the initial guess of the iterative procedure.
For the positive sequence, it is worth noting that the only unknown value is the reactive power injected q g,P in the section g, since it depends on the power flow solution. Thus, all the generator reactive power initial guesses q b ÷q g (for k = 0) must be estimated in detail. The discussion on the choice of the initial reactive power of the generators is the same of [25]: the network is firstly considered as ideal (the matrix Y ABC N is imaginary), and (8) is applied. Then, after the Fortescue transformation, the vector of the positive sequence reactive power injected by the generators can be computed. Fig. 13 summarizes the reactive power estimation procedure.
In this paper, the three-phase power flow region of attraction of the solution can be also assessed [40]. In fact, also in this algorithm, the only initial guesses are the reactive power q b ÷q g, of the generator sections. Therefore, evaluations on how ITER. changes with the initial generator reactive power guess can be made (see Sect. V).

A. THE LOAD/GENERATOR CORRECTING CURRENT METHOD
After the initial guess, the correcting-vector iterative procedure [25] is achieved by means of the formulae (21)÷ (25) and depicted in the flow-chart of Fig. 14.
This iterative scheme is basically the generalization of the procedures in [25].
The formulae (21)÷(25) are represented for the k-th iteration: Eq. (25) represents the positive sequence quadrature component current vector that springs out the same injection of positive sequence reactive power of (24). These generator current injections computed in (25) are more suitable to correct the positive sequence reactive power which are the real unknowns of the three-phase power flow (the positive sequence active power are constrained, as explained in Sect. II). This is the three-phase generalization of the concept already presented in PFPD [25].
Eq. (25) must not be applied for the slack section since both the active and reactive power are unknown quantities: For each iteration, (21) and (22) allow computing the phase voltages of the generator and load sections, and can be obtained considering the generalization of (5) and (6).
Eq. (23), (24), and (25) allow computing the positive sequence correcting currents to be applied iteratively in each load and generation section. These correcting currents are computed starting from the technical constraints given to the sequence component frame of reference. Afterwards, these currents must be converted in the phase frame of reference to be embedded into (21) and (22) of the subsequent (k+1)-th iteration.
Differently from the other methods, these formulations are defined in both the sequence component and phase frame of reference. Therefore, the impact of the system unbalance due to the presence of negative and zero sequence currents/ voltages can be assessed. For this reason, this three-phase approach is referenced as a hybrid one. Fig. 14 schematically represents the three-phase power flow algorithm flow-chart. It is worth noting the alternations of the phase formulations with the sequence component ones through the F and F −1 transformations for each cycle (once again the hybrid approach).
In Fig. 14, it can be seen that the correcting currents for k = 0 are all set to zero except for the slack bus sub-vector. After using (21), the application of F allows computing the generator voltages in sequence component frame of reference. In this way, the T x transformation allows correcting the positive sequence voltage of all the generators: a new corrected vector u G1,c can be built by changing the calculated positive sequence voltage magnitudes with the constrained ones, but by keeping δ b . . . δ g angles unchanged.
Afterwards, the application of F −1 allows passing to the phase frame of reference and applying (22) to compute the phase voltage vector of loads u ABC L,c . Thus, all the phase components of the voltage vectors are computed. These block vectors are necessary to compute the correcting current vectors by means of F and F −1 transformations, and (23), (24), (25). The procedure is iterated until convergence, i.e., until any mismatch of positive sequence generator voltage is within the tolerance.

B. THE GENERATOR CORRECTING CURRENT METHOD
So far, the algorithm is based on the calculation of two correcting current sets (one for the generator sections i ABC G , and one for the load sections i ABC L ). In order to reduce the computational cost of each iterative cycle, this section investigates the possibility of using an alternative iterative scheme, by exploiting only the generator correcting current set.
Thus, by posing i ABC L = 0 the above-mentioned formulation becomes: i G,c_q,P = −j Im u G,c,P ⊗ i * G,c,P u * G,c,P (28) Eq. (28) (which is the same of (25)) must not be applied to the slack section since the active and reactive power are unknown quantities: i G,c_q,P (1, 1) = i G,c,P (1, 1).
As it will described in the following section, the series of (26), (27), and (28) allows computing the power flow, bringing to a CPU-Time and ITER reduction.  The authors think that this result is remarkable since a three-phase power flow solution is achieved by applying only three iterated formulae. Fig. 15 shows the flow-chart of the algorithm without correcting the load sections.
No matrix inversion is necessary for each iteration, differently from [24], and this advantage is fundamental when considering three-phase power flows.

VI. PFPD_3P PERFORMANCE ASSESSMENTS A. PRELIMINARY COMPARISON BETWEEN THE NEW PFPD_3P AND THE OLD PF[24]_3P
In order to check the improvements and the enrichments of PFPD_3P versus the old algorithm [24] (which is called as PF [24]_3P in the following), a performance comparison between their two implementations is shown. These two algorithms are implemented in Matlab environment and their performance are assessed under the same conditions. Several fictitious transmission networks are tested, which differ in topologies and load/generation scenarios. Table 1 shows ITER. and CPU-time comparisons between PFPD_3P and PF [24]_3P for a voltage tolerance equal to 2 mV, by testing the 18-section network (see Fig. 16 for the representation and App. II for the description) and other test networks. The convergence tolerance used is very small to produce a very accurate power flow solution. In these networks, all the loads are considered as ''equivalent loads'' (see Sect. III A3), hence only the generator correcting current method (see Sect. IV B) is employed. Furthermore, it is worth underlining that the old PF [24]_3P allowed only the generator correcting current method and not also the generator/load correcting current one. The tests are performed in a PC using an Intel (R) Core (TM) i7-7700K CPU @ 4.2 GHz processor, with a 32 GB RAM. By observing Table 1, it is possible to note a CPU-time improvement for PFPD_3P, and this fact is especially observed for the 300 section case: a considerable CPU time reduction is achieved. In fact, the novel iterative cycle is composed of only three formulae, i.e., (26)÷ (28), and this synthetic formulation is due to the three-phase slack generator inclusion inside the ''all-inclusive'' admittance matrix. Hence, the single-phase improvements of the algorithm presented in [25] are still valid in the three-phase generalisation presented in this paper. Moreover, in order to check the PFPD_3P solution consistency, extensive comparisons with the solutions of PF [24]_3P show negligible maximum solution differences (i.e., 10 −7 p.u. for the phase magnitudes, and 10 −6 deg. for the phase angles).

B. COMPARISON BETWEEN THE NEW PFPD_3P AND DGS
In order to validate PFPD_3P with a reliable software benchmark, all the analysed networks are also tested in the commercial software DGS. The DGS software implements the Newton-Raphson method to compute the three-phase power flow solution. As in [25] and in [41], a self-made interface procedure to automatically pass all the network data from Matlab environment to DGS software is exploited. Firstly, all the power flow solution differences between PFPD_3P and DGS are assessed, and maximum differences of the orders of magnitude equal to 10 −3 p.u. for the phase magnitudes, and 10 −2 deg. for the phase angles are found. The minimum differences are instead of the orders of magnitude of 10 −6 p.u. for the phase magnitudes, and 10 −4 deg. for the phase angles. Such differences confirm the very good agreement between PFPD_3P and DGS solutions. Table 2 compares ITER and CPU-time of PFPD_3P and DGS by considering both the generator and generator/load correcting current methods (see Sect. IV A, and Sect. IV B, respectively), thus the 18-section network loads are considered both as lumped and equivalent. When the network loads are considered as equivalent ones, it is not possible to run the simulations with the generator/load correcting current method (grey area in Table 2). Table 3 compares ITER and CPU-time of PFPD_3P and DGS by considering both the iterative approaches for four tested networks (18, 41, 60, 300 section networks). The 300-section network is a three-phase version of the classical    single-phase IEEE-300 bus network, and it is chosen to demonstrate the efficiency of PFPD_3P when dealing with large three-phase networks. Eventually, a 600 section network is considered for further testing the algorithm. Eventually, to use the same convergence criterion adopted by DGS, all the iterative cycles of Table 2 and 3 are stopped for a maximum acceptable power flow mismatch equal to 1 kVA. Table 2 and 3 allow inferring the following considerations: • The CPU times of PFPD_3P are always lower than the ones of DGS for the first four networks. However, the last case study shows a slightly trend reversal: this is due to the fact that DGS is an optimized software compared with a self-implemented Matlab-based algorithm. This means that each PFPD_3P iterative cycle (i.e., (21)÷(25) or (26)÷(28)) is faster than each DGS one.
• The generator correcting current method allows reducing the CPU-Time and ITER compared to the load/generator correcting current method. In fact, in the generator correcting current method, only three-iterated formulae are used (i.e., (26)÷(28)). About the two different above-mentioned iterative approaches, also their regions of attraction are different. This concept is graphically represented, for the 18-section network, in Fig 17. It is worth noting that the region of attraction for the generator correcting current method is almost always beneath the region of attraction of the other iterative method. But most importantly, the region of attraction of the generator correcting current method is more variable/sensitive. This is because, in this approach, the iterative scheme is based on the computation of the correcting currents in the generator section only. Thus, the initial generator reactive power estimation has a greater impact in the correcting current calculation, and consequently in the convergence achievement. VOLUME 9, 2021

VII. POWER QUALITY EVALUATIONS BY MEANS OF PDPD_3P
In this section, some power quality evaluations on the 18-section fictitious network are presented. These evaluations are increasingly important for the challenging future networks. Moreover, the remarks/comments for the 18-section network are still valid/confirmed by the other tested networks. Therefore, PFPD_3P makes possible the following technical evaluations:

A. POWER QUALITY VARIATION WITH THE LENGTH OF THE LINES
The great majority of HV/EHV electrical lines are not transposed [20,22]. This fact contributes to the power system unbalance increase, and this effect is stronger the longer the transmission lines are. Fig. 18 shows, for instance, the section voltage unbalance factors for different line lengths: all the considered line lengths are the 50%, or 100%, or 150% of the 18-section network ones. Fig. 18 allows confirming how much the voltage unbalance factors grow in all the sections (both generation and load sections) as the line lengths increase.

B. POWER QUALITY VARIATION WITH THE POWER SYSTEM LOADING
In order to have a safe power system operation, the load sections should not be characterised by excessively high power   absorptions. One of the reasons is due to the power system power quality. In fact, the more power the system absorbs, the higher the unbalance factors are (this is due to the higher current circulation in the power lines). Fig. 19 shows the section voltage unbalance factors for different values of the apparent power absorbed by each section (the apparent power values are the 50%, 75%, 100%, 125% of the 18-section network ones: the scheduled active power of all the generators are also changed proportionally to guarantee the balance between generation and loading). It can be noted that the voltage unbalance factors tend to grow in the network sections as the network loading increases.

C. POWER QUALITY VARIATION WITH THE PERCENTAGE OF ASYNCHRONOUS SHARE
The load composition affects the unbalance factors of the system. In fact, the presence of a remarkable percentage of electrical machines into a load section have positive effects on the network power quality.
It is known, in fact, that the rotating electrical machinery can play a key role in the lowering of the negative sequence voltage levels, due to their high negative sequence admittance values.
For instance, Fig. 20 shows a voltage unbalance factor comparison (considering the 18-section network) by considering all the loads with the asynchronous shares equal to 0% (static load) and 60% respectively, at the same conditions (same power flow technical constraints and same network data). It is worth noting that, in the case of an asynchronous share percentage of 0%, the voltage unbalance factors are higher. Physically, this is due to the higher negative sequence admittances of the asynchronous loads, which drain the negative sequence currents. Fig. 21 shows this concept: in the load sections the absorbed current unbalance factors grow (red line) when the percentage of the asynchronous share is turned up to 60%. The voltage symmetrisation, in fact, is due to the higher negative sequence current drain in the load sections (the transit sections are excluded from the representation of Fig. 21, since no current drain can occur in these network sections).

D. POWER QUALITY IMPROVEMENT BY MEANS OF SYNCHRONOUS COMPENSATORS
Synchronous compensators could play an important role in the future transmission networks, although the massive penetration of static devices which regulate the reactive power. In fact, the synchronous machines, in addition to regulating the reactive power flows, can lower the network voltage unbalance factors. Fig. 22 shows, for instance, how the voltage unbalance factor decreases in all the network sections when a unique synchronous compensator is applied at section 9 (see Fig. 17).
The presence of negative sequence currents in the network may have detrimental impacts on rotating equipment, i.e., synchronous machines or induction motors. In fact, negative sequence currents cause rotating field moving in opposite direction to the rotor rotation. These fields produce pulsating torques/vibrations, induced currents, and machinery temperature rise which reduce both the machine life and efficiency. Table 4 shows the negative sequence current flows in rotating equipment, by considering asynchronous shares equal to 0% and 60%. The case where a synchronous compensator is inserted at section 9 of the 18-section network is also shown. In addition to the considerations on power quality, an important industrial fall out can be derived from the application of PFPD_3P, e.g., on the operation of the protections. In fact, suitable protections are installed to preserve the generators from continuous negative sequence current, which are typically set in a range 8÷15% of the machine rated current [20]. Therefore, the ability to compute the negative sequence current values is a basic knowledge to predict protection behaviours.

E. POWER LOSSES COMPARISONS BY CONSIDERING DIFFERENT GROUND WIRE EARTHING METHODS
The passive conductors of the AC electrical lines are subjected to induced voltages. Therefore, induced currents can circulate in such conductors, and their entities depend on the passive conductor earthing methods. These currents globally affect the power losses of the power system and their evaluations can be assessed by the present PFPF_3P. The 18-section network is considered, and all the electrical lines are treated as OHLs (see Fig. 23). For this network, two different ground    wire earthing arrangements are considered: the former with all the ground wires unearthed and insulated from the towers and from the substation earthing grids, and the latter with all the ground wires earthed.
These two configurations are modelled by considering the two different procedures expounded in Sect. III B,1 (Fig. 12),  and by considering suitable resistance values to model the earth wire links along the lines by means of MCA.
In particular, the single-circuit line (see Fig. 23a) has two different ground wires (i.e., a steel and an OPGW ground wires: r = 2.66 /km and r = 0.28 /km, respectively), whereas the double-circuit line (see Fig. 23b) has only one steel ground wire (r = 2.66 /km). Moreover, the tower footing resistance (r tower = 15 ), and the soil resistivity (r soil = 100 m) are supposed to be equal for the two cases. The main difference between the two configurations is the ground wire/tower contact resistance which is r 1cr = 1 G for the insulated ground-wire towers, whereas r so = 1 m for the earthed ground wires. In Table 5 the power loss comparison between these two ground wire configurations is carried out: the network power losses, by considering all the ground wires earthed, are 3.17% greater than the case with the insulated ground wires. This percentage increase is the 0.03% of the global active power generated in the 18-section network (2.3722 GW). These simulations are performed by using the load/generation correcting current method with a convergence tolerance of 10 −7 p.u. However, a thorough power loss evaluation can be directly carried out by knowing exactly the steady state-regime of the passive conductors (that can be assessed with a passive conductor power flow tool, which it has never been developed). Therefore, the global losses can be estimated thoroughly. Moreover, the network losses due only to the ground wire presence can be estimated. A further research could combine PFPD_3P with the MCA method in order to achieve this goal.

VIII. OPEN QUESTIONS
From the knowledge of the phase voltages in all the sections of the analysed network, all the section phase currents are determinable by means of the nodal (three-phase) admittance matrix. This is valid only for the active conductors (the phases), but obviously not for the passive ones (for instance the ground wires of OHLs or metallic screens of insulated cables). Nevertheless, in this paper, the effects due to those passive conductors can be included in the admittance matrix modelling such lines, by the Kron's reduction matrix techniques [35]. However, further researches are developing a multiconductor power flow, by combining the present technique with the MCA [27]- [29], [42], [43].
Thus, that future multiconductor power flow would allow directly considering both the currents and voltages in each section of all the active and passive conductors (metallic screens, ground wires, or enclosures of GILs [37], [38]) of the transmission lines. That future result could be conceived as the final generalization of this research.

IX. CONCLUSION
A novel hybrid (with the simultaneous use of phase and sequence component frames of reference) three-phase power flow algorithm based on the slack generator treatment as a quasi-ideal current source, at positive sequence, is presented. This open algorithm combines the merits of each frame of reference without englobing the demerits deriving from the application of a unique frame of reference. Two different iterative procedures are developed (i.e., the generator and the load/generator correcting current methods), and their validation is carried out by means of extensive solution comparisons with the old power flow algorithm PF [24]_3P, developed in 2000, and the commercial software DIgSILENT PowerFactory. Moreover, the paper shows a strong CPU-time improvement compared with the other methods. This is due to the low number of formulae for each iteration of PFPD_3P. By using the generator correcting current method only, it is possible to include in the model also large radial subtrasmission/distribution systems by considering them as equivalent lumped load sections, so to reduce the dimensions of the system. By this new three-phase power flow, a detailed power quality assessment of the HV/EHV network can be performed in terms of unbalance factors, and power losses also considering the passive conductors (e.g., ground wires of OHLs). A possible improvement of power quality is also presented by means of a synchronous compensator insertion into the EHV grid: the forecast of negative sequence currents into the synchronous generator and compensator sections plays a key role in investigating possible tripping of negative sequence protections. PFPD_3P, namely the three-phase generalization of the single-phase PFPD, keeps all the merits of his ''parent'': • It is an open and fast algorithm, a tool for researchers and a powerful ally to Transmission System Operators; • The convergence procedure is based on the circuit theory and not on numerical analysis techniques; • All the load, generator and slack bus sections are englobed into a unique bus admittance matrix: this is possible also for slack generator, since it is considered as a quasi-ideal current source at positive sequence; • Convergence is guaranteed by the suitable injections of two correcting current vectors into all the generators (including slack one) and all the loads: alternatively, only correcting current vectors into all the generators (including slack one) can be used in order to further speed up the algorithm; • This open algorithm allows obtaining strong performances in terms of greater precision of the solution (tolerance up to 10 −15 p.u.), shorter execution times.
In addition to all the above-mentioned advantages, the algorithm could be further enhanced by considering all the passive conductor models inside the ''all-inclusive'' three-phase admittance matrix. These models would improve the algorithm, even if all the considered (''hard'' and ''soft'' as defined in Sect. IV) hypotheses of PFPD_3P are proven to be licit for a good three-phase power flow study of any transmission network.
Further researches, in fact, are on-going towards the combination of MCA with this power flow algorithm in order to evaluate any steady-state regime of both active and passive conductors without any simplifying hypothesis. This future result could be conceived as the final generalization of this research.

APPENDIX A
In the load/generator current iterative procedure, a clarification about (23) i.e., i 0PN L,c = −Y L,P 1 − u L,c,P . 2 u L,c,P * .
is necessary. In (23), the symbol '' '' takes on the meaning of ''positive sequence multiplication''. This notation is introduced to summarize the procedure giving the sequence correcting currents due to the positive sequence voltages u L,c,P only. Therefore, u L,c,P and the elements connected to the positive sequence voltages of the load admittance submatrix Y L must be considered in (23). All the elements connected to the positive sequence voltages of Y L are stored in the matrix indicated as Y L,P . In order to explain how to compute Y L,P , an example involving only the (3 × 3) block load admittance Y m modelling the m-th load is considered. By making a simple consideration on the matrix/vector multiplication, it is clear that the impact of the positive sequence voltage u m,P on the sequence correcting currents i 0PN m,c is due to the second column of the matrix Y m only (see Fig. 24). Therefore, each element of i 0PN m,c is simply the product of the second column of Y m for the scalar u m,P . However, (23) indicates that the above-mentioned procedure must be applied for each load section: thus, the complete vector i 0PN L,c representative of all the load sections can be computed, by means of (23) uniquely.
Eventually, it is worth clarifying that the squaring present in (23) is a point-wise squaring (each element of the vector u L,c,P must be squared).

APPENDIX B
The 18-section network (see Fig. 16) consists of three generator sections, fifteen load sections (eight of which are transit sections) connected each other by means of fifteen non-transposed electrical lines. The coexistence of both single and double (the line between sections 7 and 8) circuits is considered. All the loads are three-phase symmetrical ones. In section 9, it is possible to insert a synchronous compensator. Moreover, three-phase power flow studies with load sections or equivalent load sections can be assessed and the simulations can be carried out for both cases. In particular, when the equivalent loads are considered, the subtransmission/distribution network underlying their sections are characterized by the presence of seven feeders and seven transformers supplying seven loads. About that, Fig. 16 shows that the equivalent load in section 18 allows reducing all the fifteen underlying sections (from 19 to 33) in a unique lumped one. Therefore, by considering all the loads (transit ones excluded) of the 18-section network as equivalent, one hundred and five sections are globally reduced into seven sections only.