A Three-Phase Weather-Dependent Power Flow Approach for 4-Wire Multi-Grounded Unbalanced Microgrids With Bare Overhead Conductors

Conventional power flow algorithms assume that the network resistances and reactances remain constant regardless of the weather and loading conditions. Although the impact of the weather in power flow analysis has been recently investigated via weather-dependent power flow (WDPF) approaches, the magnetic effects in the core of Aluminum Conductor Steel Reinforced (ACSR) conductors have not been explicitly considered. ACSR conductors are widely used in distribution networks. Therefore, this manuscript proposes a three-phase weather-dependent power flow algorithm for 4-wire multi-grounded unbalanced microgrids (MGs), which takes into consideration the impact of weather as well as the magnetic effects in the core of ACSR conductors. It is shown that the magnetic effects in the core can significantly influence the power flow results, especially for networks composed of single-layer ACSR conductors. Furthermore, the proposed algorithm explicitly considers the multi-grounded neutral conductor, thus it can precisely simulate unbalanced low voltage (LV) and medium voltage (MV) networks. In addition, the proposed approach is generic and can be applied in both grid-connected and islanded networks. Simulations conducted in a 25-Bus unbalanced LV microgrid (MG) highlight the accuracy and benefit of the proposed approach, while its computation performance is tested in the IEEE 8500-Node network.


I. INTRODUCTION
P OWER flow is a fundamental tool for several power system studies such as planning, design and operation, state estimation, contingency and stability analysis, etc. Therefore, the accuracy of the power flow is very important for a better Manuscript received April 19, 2020; revised August 4, 2020; accepted October 1, 2020. Date of publication October 6, 2020; date of current version April 19, 2021. This research is co-financed by Greece  understanding of power systems and improving their analysis. Conventional power flow (PF) is based on the assumption of fixed line impedance (the line impedances have pre-specified and invariable values) without considering the influence of weather and magnetic effects in the core of the conductors. As such, inaccuracies resulting from the weather and magnetic effects are present in PF. The importance of incorporating weather conditions into the power flow analysis has been demonstrated in various literature studies [1]- [7]. These algorithms present the advantage of more accurate evaluation of power system states, power losses, and power flows in the network. In addition, they exhibit potential benefits in power system planning, analysis, and operation such as network upgrade deferment, reduction of total generation cost, maximization of useable capacity of overhead lines, improvement of transient stability analysis, etc. [1]- [9].
However, the algorithms presented in the existing literature do not fully consider the nonlinearities of line impedance as a result of the magnetic effects in the core of bare overhead conductors. For instance, the recently proposed weather-dependent power flow (WDPF) [4] algorithm considers explicitly the impact of the weather parameters, but assumes a linear resistance-temperature relationship. Furthermore, the existing WDPF algorithms [1]- [7] are for single-phase analysis of balanced transmission networks. As such, they are not suitable for three-phase analysis of LV and MV networks. Similarly, they are also not suitable for the analysis of islanded microgrids, due to the particular characteristics the islanding presents in power flow analysis, such as the absence of a slack bus and the droop-control of distributed generators (DGs) [10], [11].
In this paper, we highlight the nonlinear relationship that exists between the impedance, temperature, and current of overhead conductors. A large portion of the electric network is composed of bare overhead conductors due to their lower cost in comparison to cable installations [12]. These conductors consist of several aluminum (or aluminum alloy) strands wrapped helically around a central core. Depending on the material of the core, they are sub-divided into two major categories: a) All Aluminum Conductors (AAC) with an aluminum core of low relative permeability (μ Al = 1) and b) Aluminum Conductor Steel Reinforced (ACSR) with a steel core of high relative permeability (μ s >> 1).
As a result of the spiraling structure of the aluminum strands around the core in both types of conductors, an alternating 0885-8950 © 2020 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission.
longitudinal magnetic field is created in the core [13]. The extent of this longitudinal field depends on the lay length and the number of layers of the helical aluminum strands [14]. With an even number of layers, there is a partial cancellation of the magnetic field due to the opposite spiraling directions of the aluminum layers [13]. On the other hand, with an odd number of layers, there is a significant resultant magnetic field in the core that creates hysteresis and eddy current losses increasing the conductor resistance and reactance nonlinearly [13], [14]. While the magnetic effect can be neglected in AAC conductors due to the low relative permeability of the core, it is significant in ACSR conductors with an odd number of layers, especially the single-layer [13], [14]. Single-layer conductors are mainly used in distribution networks due to their lower current capacity, between 100A and 400A [15]. As a result, the magnetic effects are more pronounced in LV and MV networks. Therefore, modeling and analysis considering the nonlinear relation between the conductor impedance, conductor temperature, and loading as a result of the magnetic effects in the core is crucial for an accurate distribution system analysis.
To date, a three-phase weather-dependent power flow to study 4-wire multi-grounded unbalanced MGs has not been proposed. This manuscript proposes a more realistic modeling of unbalanced microgrids incorporating DGs and the effect of weather conditions as well as the nonlinear magnetic effects in the core of ACSR conductors. Furthermore, the 4-wire multi-grounded nature of distribution networks is realistically represented, thus the algorithm can precisely simulate LV and MV networks. It is worth mentioning that the proposed power flow approach is generic and applicable in all kind of configurations (meshed and radial) and operational modes (grid-connected and islanded). A series of power flow simulations are executed to highlight the novel findings in power flow analysis via this approach.
The rest of the paper is structured as follows: Section II briefly describes the power flow approach proposed in [10], which is the power flow solver utilized. Section III outlines the nonlinear heat-balance equation of bare overhead conductors, while section IV models, more accurately, the nonlinear relationship between the impedance, conductor temperature, and conductor current by taking into consideration the magnetic effects in the core of the conductor. Section V summarizes the steps of the proposed power flow approach. Section VI presents simulation results to highlight the contribution of the proposed power flow approach. Section VII investigates the computational performance of the proposed approach in a large network, while Section VIII concludes the paper.

II. CONVENTIONAL POWER FLOW ALGORITHM
This section includes a short description of the power flow solver applied in this paper. More details about the solver are provided in [10]. The algorithm is comprehensive and can be applied in both grid-connected and islanded networks, regardless of the network configuration.

A. Theoretical Background of Islanded MGs
MG is a group of interconnected loads and distributed energy resources within clearly defined electrical boundaries that acts as a single controllable entity with respect to the grid [16]. It can operate either in grid-connected or in islanded mode. In the first case, it is connected to the main grid, while in the second case it is intentionally or unintentionally disconnected from the main grid forming an autonomous network [10].
In Islanded mode, the DGs operate in droop control to equally share the active and reactive power determined by the P-f and Q-V droop curves, as follows [10]: where f, f 0i , K pi , P Gi , V Gi , V 0i , K Qi and Q Gi are the output frequency (Hz), reference frequency (Hz), frequency droop gain (Hz/W), active power output (W), positive sequence voltage magnitude (V), reference voltage (V), voltage droop gain (V/VAR) and reactive power output (VAR) of DG i, respectively. As a result of the droop operation, the consideration of a conventional slack bus for power flow analysis is invalid in islanded MGs.

B. The Concept of Virtual Slack Node [10]
In the absence of a real slack node in islanded MGs, a virtual slack node is assumed to be connected in the network so that the formulation of power flow equations is the same as in gridconnected mode. The only difference is that in islanded mode, the voltage of the virtual slack node is not constant but equalized in each iteration k+1 with the voltage of its adjacent node as it was calculated in the previous iteration k. More specifically, assuming the islanded network of Fig. 1, the voltages of the virtual slack node (e.g V a(vir) , V b(vir) , V c(vir) , V n(vir) , V g(vir) ) are equalized with their adjacent voltages according to Equation (3).
Moreover, in islanded operation the active power that flows through the virtual slack node is calculated in each iteration and it is assigned to DGs based on their droop equations. As a result, after the algorithm has converged, the virtual slack node acquires the same voltage as its adjacent node and can therefore be discarded (no physical significance). More details are provided in [10], Section III].
It is important to note that the virtual slack node concept is very convenient in networks with frequent transition between grid-connected and islanded mode, since it allows the same power flow equations to be applied in both modes. The only difference is that in grid connected mode, the slack node has a constant reference voltage (based on conventional power flow assumptions), while in islanded mode it is updated in every iteration until convergence.

C. Power Flow Equations
The full configuration of an unbalanced LV network is presented in Fig. 1, including the neutral conductor and the grounding resistances. Let us define the current and voltage vectors of node i as follows: where I Lir and V iy denote the load current and voltage (in complex form) of node i at phase r = {a, b, c} and conductor y = {a, b, c, n, g}, respectively, as shown in Fig. 1. In Fig. 1, V y(vir) and Z y(vir) denote the voltage and impedance of conductor y = {a, b, c, n, g} of the virtual slack node, respectively. For a network with m nodes, the load vectors can be expressed as a function of the voltage vectors according to Equation (6): where V 0 coresponds to the virtual slack node. Subsequently, for the power flow solution, we remove the first five rows of Equation (6) that correspond to the virtual slack node and Equation (7) is obtained.
Then, by transferring all the voltage variables of the left-hand side of Equation (7) to the right-hand side, Equation (8) is derived, where I new and Y new are the modified current and admittance matrices.
As a last step, we define the final matrices Y fin and Y fin . The first one consists of the first five columns of Y new , while the second one consists of the remaining columns so that (9), we finally derive Equation (10), which is iteratively solved until a certain preset tolerance is reached. In Equation (10), k denotes the iteration number, and the vector V fin contains the voltages of all nodes except the virtual one.
The vector V virt is fixed in grid-connected mode. For example, in grid-connected LV networks V virt = [ 230, α · 230, α 2 · 230, 0, 0], where α is the phasor rotation operator α = e j 2 3 ·π . In islanded mode, V virt is equalized in each iteration with its adjacent node, according to Equation (3). Further details about this power flow algorithm are provided in [10], while its applicability is extended in [17]- [19] to hybrid AC/DC networks and to networks containing voltage regulating devices e.g. step voltage regulators (SVR) or on-load tap changer transformers (OLTC).

III. THERMAL MODELING OF OVERHEAD CONDUCTORS
The thermal modeling of conductors is crucial for incorporating the effect of weather in power flow analysis. The temperature of an overhead conductor can be calculated from the nonlinear heat balance equation of IEEE Std 738-2012 [4], [20]. The nonlinear thermal equation is as follows: where Q c is the convective heat loss, Q r is the radiative heat loss, Q s is the heat gain due to solar radiation, and Q j is the heat gain due to joule heating. IEEE Std 738-2012 provides several equations for calculating the convective heat loss [4], [20]. The convective heat loss (Q c ) is caused either by forced convection or natural convection depending on the wind speed (V w ). Two alternate equations are recommended by IEEE Std 738-2012 for the calculation of forced convection, as shown in Equation (12) and Equation (13), respectively.
The equation for natural convection (Q cn ) is as follows: In order to obtain the value of Q c , Equations (12), (13), and (14) are solved and the maximum value is considered i.e. max {Q c1 , Q c2 , Q cn }, as recommended by IEEE Std 738-2012 [20].
The radiative loss (Q r ) is given by Equation (15), and includes the losses transmitted from the conductor to the environment by radiation.
The solar heat gain (Q s ) is given by Equation (16), and concerns the energy transferred from the sun to the conductor through solar radiation.
The Joule loss (Q j ) is given by Equation (17), and refers to the active power losses generated in the conductor due to the current flow.
The nonlinear heat-balance Equation (11) can be solved to obtain the conductor temperature (T c ) for any amount of current flow under any given weather condition experienced by the conductor. To solve the heat-balance equation, a heat-balance function H(T c ) is defined as follows: The heat-balance function in Equation (18) is then solved using Newton's method to obtain the conductor temperature (T c ) as follows: In Equation (19), k is the iteration number, H (T c ) k is the derivative of the heat-balance Equation (18) evaluated at the k th iteration. After the conductor temperature (T c ) is obtained, accurate calculation of the conductor impedance is performed by considering the nonlinear relationship between conductor impedance, conductor temperature, and current via modeling the magnetic effects in the core of the conductor, as presented in the following section.

A. Theoretical Analysis of the Applied Model
IEEE Std 738-2012 [20] assumes a simplified linear resistance-temperature relationship that ignores the nonlinearities arising from the magnetic effects in the core of the conductors. In reality, both ACSR and AAC overhead conductors consist of many aluminum wires stranded helically around a central core. Thus, the current flowing through these wires follows a helical path resulting in an additional longitudinal magnetic field in the core of the conductor. If the core consists of a ferrous material (e.g. steel), which is the case in ACSR conductors, additional losses are created in the core.
The aforementioned effect is most pronounced in single-layer ACSR conductors. This is because in multi-layer conductors a partial cancellation occurs due to the alternating directions of lay of the different layers of aluminum wires [13]. Since all ACSR conductors with a current capacity less than 400A are single-layer conductors [15], these are commonly used in LV and MV overhead distribution networks. Therefore, for a more precise analysis of distribution networks, the accurate estimation of conductor impedance through the inclusion of magnetic effects and weather-dependent impacts is crucial.
A precise model [14], to accurately calculate the impedance of a three-layer ACSR conductor considering the magnetic properties of the steel core, is depicted in Fig. 2. The authors in [14] have experimentally validated the model. In Fig. 2, R S is the resistance of the conductor core, while R 1 , R 2 , and R 3 , are the DC resistances of layer 1, 2, and 3, respectively. I S is the current in the core, while I 1 , I 2 , and I 3 , are the currents of layer 1, 2, and 3, respectively. The inductances of the conductor due to the longitudinal magnetic field are shown in the horizontally parallel branches, while the inductances due to the circular magnetic field are shown in the vertical branches of Fig. 2.
The longitudinal reactances are calculated as [14]: In equation (20), ω is the angular frequency, α s the crosssection area of the steel core, α L is the cross-section area enclosed by layer p, λ p the lay length of layer p, λ q the lay length of layer q [14]. The longitudinal reactance includes the reactive and active power (e.g. hysteresis and eddy current losses) losses in the core of the conductor.
The reactive and active power losses are separately represented in Equation (20) by means of a complex relative permeability (μ s ) instead of a real one. The complex relative permeability of the steel (μ s ), in Equation (20), consists of a real and an imaginary part [14] as follows: The real-part (μ s ) is associated with the reactive losses in the steel core, while the imaginary-part (μ s ) with the active losses. According to [21], (μ s ) and (μ s ) are dependent on the core temperature and the longitudinal magnetic field in the steel core of the ACSR conductor. For a typical steel core, both (μ s ) and (μ s ) are experimentally calculated in [21] and are  depicted in Fig. 3, as a function of the longitudinal magnetic field (H) and temperature. The longitudinal magnetic field (H) in the core of a threelayer conductor, as the one in Fig. 2, is calculated as follows [14]: In Equation (22), λ 1 , λ 2 , and λ 3 , are the lay lengths of layer 1, 2, and 3, respectively. The negative sign in the second term of the right-hand side of Equation (22) is due to the reverse spiral direction of the second layer, resulting in a partial cancelation of the longitudinal magnetic field. However, in single-layer overhead conductors, only the first term of the right-hand side of Equation (22) exists, thus a small current causes a significant longitudinal magnetic field in the core of a single-layer conductor.
In Fig. 3, (μ s ) and (μ s ) of the steel core are observed to be highly nonlinear as the longitudinal magnetic field in the core increases at different conductor temperatures. Thus, the active and reactive losses (therefore the resistance and inductance) vary nonlinearly with the increasing current (due to the longitudinal magnetic field) at different conductor temperatures.

B. Equations of Resistance and Reactance
Proceeding from the previous section, the calculation of the nonlinear conductor impedance is undertaken here, by considering also the magnetic properties of the core of ACSR and AAC conductors.
Assuming a single-layer ACSR conductor, only the first two layers of Fig. 2 are considered, namely the core layer (the layer of current I s ) and the first aluminum layer (the layer of current I 1 ). Equating the voltages of these two layers, Equation (23) is obtained.
The conductor current (I CON ) is subdivided to the steel and aluminum layer according to Equation (24).
Equations (23) and (24) can be written in a matrix form as follows: Equation (25) calculates the current distribution between the layers (I S , I 1 ) for a given conductor current (I CON ). Equation (25) is nonlinear since the parameter X 11 (see Equation (20)) includes the term μ s , which is dependent on the magnetic field inside the core (see Fig. 3), which in turn is dependent on the current I 1 (see Equation (22)). In other words, the current I 1 affects nonlinearly the parameter X 11 through the relative
Step 3: Solve the nonlinear Equation (25) for I 1 , I S .
Step 4: Calculate voltage (V CON ) from Equation (26) Step 5: Calculate R CON and X INT from Equations (27)- (28) and save them in the matrices R CON (I CON , T c ) and X INT (I CON , T c ).
Step 6: If I CON < I max increase I CON by 1 and go to step 3, otherwise go to step 7.
Step 7: If T C < T max increase T C by 1 and go to step 2, otherwise go to step 8.
Step 8: Take the completed matrices R CON (I CON , T c ) and X INT (I CON , T c ) for the desired current and temperature range.
The conductor voltage is then calculated from Equation (26) based on the circuit of Fig. 2 (25)).

and the layer currents (Equation
The resistance and internal reactance of the conductor are calculated from Equation (27) and Equation (28) [13], respectively, which include the longitudinal and circular magnetic fields, the magnetic core losses, and the DC resistance of the strands.
An algorithm, to calculate the R CON and X INT for a wide range of conductor currents (from 0 to I max ) and temperatures (from 0 to T max ), is presented in Algorithm 1. Algorithm 1 is executed offline and outputs two 2-dimensional matrices i.e. R CON (I CON , T c ) and X INT (I CON , T c ), which include the values of resistance and internal reactance, respectively, for every selected combination of conductor current (I CON ) and temperature (T c ).
Finally, the one-foot equivalent reactance (or self-reactance) of the conductor is calculated as: where f is the frequency, and D 1 is the diameter of the conductor. The matrices X CON (I CON , T c ) and R CON (I CON , T c ) can be stored and utilized in a power flow analysis for more precise results.

C. Impact on Resistance and Reactance
Using Algorithm 1, we present the resistance and reactance of two types of conductors (ACSR and AAC), considering the  different magnetic properties of their cores. More specifically, the resistance and self-reactance of the Rabbit ACSR [22] and Poppy AAC [23] (also referred as Ant in [22]) are shown in Fig. 4 and Fig. 5, respectively, as a function of the current for various conductor temperatures. Both Rabbit ACSR and Poppy AAC are single-layer bare overhead conductors with similar geometrical characteristics. The Rabbit ACSR has 6 aluminum and 1 steel strand, while Poppy AAC has 6+1 aluminum strands.
As observed in Fig. 4 and Fig. 5, both the resistance and self-reactance vary nonlinearly with the current and temperature for Rabbit ACSR conductor. Fig. 4 depicts that the linear relationship between resistance and temperature adopted in [1]- [9] and [20] produces inaccurate estimation of resistance since the nonlinear influence of current at any given conductor temperature is neglected. Moreover, as shown in Fig. 5, the consideration of a constant reactance adopted in [1]- [9], is also not suitable for networks with single-layer ACSR conductors, since the reactance varies nonlinearly with the current at any given conductor temperature.
In contrast to Rabbit ACSR, the resistance of Poppy AAC, as shown in Fig. 4, is dependent only on the conductor temperature and varies linearly with it. The self-reactance of Poppy AAC is not dependent on the conductor temperature, and remains constant over the whole current range, as shown in Fig. 5. This is due to the low relative permeability (μ Al = 1) of the aluminum core in AAC. Thus, for AAC, the proposed method presents identical results to [1]- [9] and [20].

V. PROPOSED WEATHER-DEPENDENT POWER FLOW
The proposed three-phase weather-dependent power flow approach is outlined as follows: Step 1: Firstly, the matrices R CON (I CON , T c ) and X CON (I CON , T c ) are calculated and saved offline, according to Algorithm 1, for a wide range of expected conductor currents and temperatures. Thus, the accurate impedance value for every current-temperature combination is available.
Step 2: Then, the iterative process begins by calculating the three-phase power flow using the algorithm described in Section II.
Step 3: From the currents obtained via the power flow solution in Step 2, we calculate the conductor temperatures by solving the heat-balance model described in Section III.
Step 4: Using the conductor currents obtained in Step 2 and conductor temperatures in Step 3, we extract the corresponding value of impedance of each branch conductor from the R CON (I CON , T c ) and X CON (I CON , T c ) matrices. This iterative process is repeated from Step 2 until convergence.

VI. SIMULATION STUDIES & RESULTS ANALYSIS
In this section, various simulation results are presented, to contrast the proposed power flow approach against the PF [10] and the WDPF algorithms, using a 25-Bus unbalanced LV MG. Moreover, the difference in the thermal rating and conductor sag, between the proposed approach and IEEE Std 738-2012 [20], is also presented. Finally, the computational times of the algorithms are compared.

A. Network Description
The topology of the 25-Bus unbalanced LV MG is depicted in Fig. 6 [10]. It operates either in grid-connected or islanded mode. Every node supplies an unbalanced three-phase load with a power equal to (P a , P b , P c ) = (1.3 kW, 1.5 kW, 1.7 kW). All the loads operate with an inductive power factor 0.8.
Three DGs are connected to the nodes 13, 19, and 25, respectively. The DGs operate in constant PQ mode with a unity power factor during the grid-connected operation and in droop-control during the islanded operation. The maximum generated power of the DGs of nodes 13, 19, 25 is 60 kW, 50 kW and 40 kW, respectively. All DGs generate balanced phase-to-neutral voltage both in grid-connected and in islanded mode to mitigate the voltage unbalances [10]. In islanded mode, each DG operates in droop control, with the following P-f and Q-V droop parameters: (f 0i , K P i ) = (50 Hz, 0.5·10 −6 Hz/W) and (V 0i , K Qi ) = (250 V, 2·10 −4 V/Var).
The Rabbit ACSR conductor is chosen for all the lines of the network. It is a single-layer 6/1 ACSR conductor with a cross-sectional area of 61.7 mm 2 and a lay-length of 120.6 mm. The length of the lines 2-6, 3-18 and 4-23 is 400 m, while all the other lines are 150 m.

B. Weather Conditions
Two base environmental conditions were considered for the simulation study: a) T a = 36.3°C, q s = 847 W/m 2 , V w = 1.33 m/s, K angle = 0.855, b) T a = 7.3°C, q s = 0 W/m 2 , V w = 4.03 m/s, K angle = 0.855. Here, T a is the ambient temperature, q s the solar irradiance, V w the wind speed, K angle the wind direction factor [20]. The first environmental condition corresponds to a typical summer day as measured at a weather station in Sindos (a suburb of Thessaloniki, Greece, with zero elevation) on 23/07/2018 at 14:00, while the second one corresponds to a typical winter night measured from the same station on 18/01/2018 at 21:00 [24]. The terms "summer" and "winter" in the rest of the paper denote the above-mentioned weather conditions.

C. Power Flow Simulation Results
In order to highlight the importance and accuracy of the proposed power flow approach, the PF [10] and the WDPF algorithms were coded and compared with the proposed approach in MATLAB. More specifically, the voltage profile, power losses, and voltage stability were compared.
In PF, the resistance and reactance of the conductors are considered constant without any influence of the weather on the power flow. We assume a thermal limit of 80°C for the ACSR conductors based on the recommendations of several manufacturers [22], [23]. Hence, in order for the power flow results to be on the conservative side, the corresponding resistance and reactance for PF were selected from Figs. 4 and 5 at the maximum conductor temperature of 80°C i.e. R = 0.77 Ω/km, X = 0.3298 Ω/km. In WDPF, the three-phase power flow solver of Section II was used, while the weather dependency is considered through the thermal equations of Section III. The linear resistance equation suggested by the IEEE Std 738-2012 is applied, as shown in Equation (30).
where R(T high ) = 0.88 Ω/km, R(T low ) = 0.55 Ω/km, T high = 120°C, and T low = 0°C for the Rabbit ACSR conductor. Similar to the PF, the WDPF algorithm assumes a constant reactance of X = 0.3298 Ω/km.
As mentioned earlier, the proposed power flow approach utilizes offline calculated values of resistance and reactance. Thus, the resistance and reactance are modeled exactly as shown in Fig. 4 and Fig. 5.

1) Voltage Profile of the Network:
The voltage profiles (phase c) of the network in both grid-connected ( Fig. 7(a)) and islanded ( Fig. 7(b)) mode in the summer time, for the three investigated approaches, are depicted in Fig. 7. In grid-connected mode, DGs inject their full power into the network, thus causing a voltage rise in some nodes [25]. In islanded mode, the whole energy demand is supplied by the DGs, thus causing higher voltage drops in remote nodes.
As shown in Fig. 7(a), for the grid-connected mode, deviations up to 4V are observed in some nodes between the proposed and the WDPF, and up to 6.2 V between the proposed and the PF. In islanded mode shown in Fig. 7(b), significant deviations up to 7V are observed between the proposed and the WDPF, and up to 15.5V between the proposed and the PF. These deviations cannot be neglected since they can completely mislead engineers with regards to distribution system analysis. This highlights the importance of correctly modeling the weather and nonlinear magnetic core effects of conductors in power flow analysis.

2) Power Losses of the Network:
The total losses of the network using the three investigated power flow approaches in the summer time are shown in Fig. 8. It is noted that in grid-connected mode, the active losses are increased compared to the islanded mode for two reasons: First, in grid connected mode, the DGs inject their full power into the grid increasing the line currents and therefore the losses. Second, in grid-connected mode, the DGs operate at a unity power factor, which is usually the case in real LV networks. Thus, the reactive power of the loads is transferred from the main grid increasing the losses. On the contrary, the DGs in islanded operation generate reactive power according to their droop settings, which is locally consumed by the loads.
As shown in the figure, important deviations exist in active and reactive power losses, between the three approaches. Indicatively, the active losses in islanded mode for the PF, WDPF and proposed approach are 20.44 kW, 13.82 kW and 11.04 kW, respectively. The PF deviates by 85.1% and the WDPF by 25.2%, when compared with the proposed approach.
3) Voltage Stability Analysis: The PV-curve is a very important indicator for studying the stability of a system. It indicates the maximum load limit, beyond which a voltage collapse occurs [26]. The maximum loadability of node 15 (the weakest node of the network) is investigated. The PV curve of phase c of node 15 is presented in Fig. 9(a) for the three power flow approaches, while the resistance of phase c of the line between the nodes 14-15 is presented in Fig. 9b. It should be noted that the power in the horizontal axis of Fig. 9 denotes an additional three-phase balanced load (with a cosϕ = 0.8) connected to node 15.
Based on Fig. 9(a), the three investigated approaches reach to completely different conclusions regarding the loadability of the network. More specifically for node 15, the maximum loadability estimated by the PF and WDPF is 9.69 kW and 18.86 kW, respectively. However, the proposed approach indicates a maximum loadability of 25.98 kW. This significant deviation resulted from the different estimation of the line impedances by the three investigated approaches, as shown in Fig. 9(b).

D. Comparison of Conductor Thermal Rating
In this study, a comparison of the thermal line rating (for the Rabbit ACSR conductor) obtained via IEEE Std 738-2012 and our proposed approach is performed. IEEE Std 738 utilizes a linear resistance-temperature relationship presented in Equation (30) to calculate the conductor resistance. The same linear relationship is used in the WDPF algorithm [4]. To investigate the influence of weather parameters on the thermal rating of the line, each parameter is varied stepwise, while the remaining parameters remain constant. The winter condition was considered. The results are depicted in Fig. 10. It should be noted that the maximum conductor temperature for Rabbit ACSR is 80°C [22], [23].
As shown, the wind speed is the most important parameter for the current rating. This is very important in networks with high wind DG penetration since the maximum current rating of the conductors usually coincides with the maximum DG power. On the other hand, the least influence is caused by the solar radiation, which only slightly affects the current rating of the conductor.
Generally, the moderate differences in Fig. 10 show that as far as thermal rating calculations are concerned, both IEEE Std 738-2012 and the proposed approach can be utilized with sufficient accuracy. However, the use of worst-case thermal rating, as performed in PF, is very conservative and may not yield realistic studies. For example, the manufactures in [22], [23] recommend a current rating of Rabbit conductor 269 and 190A, respectively, which is a very conservative estimation.

E. Comparison of Conductor Temperature and Sag
In this Section, we investigate the conductor temperature and the corresponding sag of the conductor calculated using IEEE Std 738-2012 and our proposed approach. Fig. 11 depicts the conductor temperature and sag of the Rabbit ACSR conductor in the winter condition, for a current up to 480 A. A maximum temperature difference between the two models of around 10°C is observed closer to 480 A of current.
For the calculation of sag, we assumed that the distance between the spans (S) is 75 m and the length of the conductor at 0°C (L 0 ) is 75.005 m. The coefficient of thermal expansion of Rabbit conductor is 19.1×10 −6°C−1 [22]. The length of the conductor is calculated by (31): where L and L 0 are the conductor lengths at T c and 0°C, respectively. The sag is calculated by Equation (32), based on the result of Equation (31), as follows [27]: As shown in Fig. 11b, IEEE Std 738-2012 overestimates the sag of the conductor by around 10 cm for currents near 480 A.

F. Comparison of Computation Time
The computation times of the three power flow approaches are presented in Table I, for an accuracy of 10 −5 pu. All simulations were executed in a PC with a 64-bit Intel Core i7, 3.4G Hz CPU and 16 GB RAM. In grid-connected mode, the PF is six times faster since the admittance matrix is formed and inverted once to be reused in every iteration. On the contrary, the WDPF and the proposed approach require an update and inversion of the admittance matrix in every iteration due to the variation of impedances. In islanded mode, the difference between the PF and the other two methods is smaller because the admittance matrix necessarily must be updated and inverted in all methods due to the variation of system frequency and therefore the inductance. It is also observed that the WDPF and the proposed approach present identical computation time, thus validating that the consideration of the magnetic effects through the proposed approach does not increase the computation time at all.

VII. APPLICABILITY IN LARGE NETWORKS
To investigate the applicability of the proposed power flow algorithm in large networks, it is applied in the IEEE 8500-Node network operated in grid-connected mode [30]. The network was modified so that all the conductors of the network are Penguin ACSR [15]. Moreover, the SVRs and capacitor banks were considered inactive. Table II depicts  r Step 2 is the most time-consuming step of the algorithm due to the iterative recalculation of the impedance matrix (namely Y −1 fin in equation (10)) of the network. The matrix decomposition technique was applied in this paper to speed up the calculation of impedance matrix [28], [29]. r Step 4 is instantaneously executed since it involves only the extraction of two values from two matrices. As shown in Table II, the proposed approach presents the same computation time as the WDPF, while it is 3 times slower than conventional PF. All algorithms converge in 9 iterations with an accuracy of 10 −5 pu, thus the total computation time of the PF, WDPF and proposed approach is 2.16, 6.39 and 6.39 seconds respectively.
An important thing to note is that the matrix decomposition is applicable only in radial or weakly meshed networks [28], [29]. In the case of highly meshed networks, the impedance matrix is computed using LU factorization. Several software packages exist, which can efficiently solve sparse matrices of large systems e.g KLU solver [32], [33]. Moreover, with the rapid progress of parallel computing, GPU-based parallel LU factorization solvers can complete the LU factorization of matrices with many thousands (even billion) of elements in a few milliseconds [31], [32].

VIII. CONCLUSION
This paper proposes a three-phase weather dependent power flow algorithm for unbalanced microgrids, which considers the impact of both the weather and the magnetic effects in the core of ACSR conductors. In this type of conductor, the resistance and reactance present a nonlinear relation with the temperature and current, which needs to be considered in the power flow for a more precise analysis.
Simulations were performed for a rural 25-bus unbalanced LV MG operated in either grid-connected or islanded mode. Simulations show that the proposed approach presents improvements in power flow results compared with the other existing approaches. The thermal rating and sag results show acceptable differences. Although, the computation time of the proposed method is increased compared to PF, it remains reasonable in light of result accuracy.