A Sensitivity ‐ Based Three ‐ Phase Weather ‐ Dependent Power Flow Algorithm for Networks with Local Voltage Controllers

: Local voltage controllers (LVCs) are important components of a modern distribution sys ‐ tem for regulating the voltage within permissible limits. This manuscript presents a sensitivity ‐ based three ‐ phase weather ‐ dependent power flow algorithm for distribution networks with LVCs. More specifically, the proposed algorithm has four distinct characteristics: (a) it considers the three ‐ phase unbalanced nature of distribution systems, (b) the operating state of LVCs is calculated using sensitivity parameters having accelerated convergence, (c) it considers the precise switching se ‐ quence of LVCs based on their reaction time delays, and (d) the nonlinear influence of weather var ‐ iations in the power flow is also taken into consideration. Simulations and validation results pre ‐ sented indicate that the proposed approach outperforms other existing algorithms with respect to the accuracy and speed of convergence, thus making it a promising power flow tool for accurate distribution system analysis. The proposed approach was compared against some state ‐ of ‐ the ‐ art power flow methods using a 7 ‐ bus, an 8 ‐ bus and the IEEE 8500 ‐ node network. The results indicate that the existing algorithms can produce inaccurate power flow results and LVC’s states, due to the imprecise consid ‐ eration of LVC’s switching sequence. Finally, it was deduced through simulations, that for a more precise estimation of LVC’s states and power flow results, the weather varia ‐ tions should be considered in the power flow. Author Contributions: Conceptualization, E.E.P.; methodology, E.E.P.; software, E.E.P.; validation, E.E.P.; formal analysis, E.E.P.; investigation, E.E.P.; resources, E.E.P.; data curation, E.E.P.; writing— original draft preparation, E.E.P. and A.A.; writing—review and editing, E.E.P. and A.A.; visualiza ‐ tion, E.E.P. and A.A.; supervision, M.C.A.; project administration, E.E.P. and M.C.A.; funding


Introduction
Local voltage controllers (LVCs) are important components in modern distribution networks. These are physical apparatus such as on-load tap changer transformer (OLTC), step voltage regulator (SVR), switched capacitor, distributed generator (DG), which regulate the voltage of a specific bus within permissible limits [1][2][3].
For safe, reliable, and efficient operation, distribution management systems (DMS) regularly carry out state estimation and optimization functions such as Volt/Var control (VVC) and optimal feeder reconfiguration (OFR). These functions require the online execution of power flow, by generating what-if scenarios. As soon as an assessment of the network's state is estimated via distribution system state estimation [4], DMS optimization functions calculate new switching states for the local controllers to optimize the network operation. The real-time power flow constitutes an integral part of VVC and OFR functions. More specifically, it is applied to estimate how the switching actions of LVCs affect the performance of the network, and accordingly, it is sent (or not) for execution via SCADA [2]. As a result, precise and fast computation of power flow in networks with LVCs is crucial for the effective execution of DMS optimization functions and state estimation of distribution networks.

Challenges of Incorporating LVCs in Power Flow Computation
Although several power flow methods have been presented in the literature for simulating the steady state of distribution networks, the overwhelming majority of them has not considered the exact operation of LVCs e.g., SVRs, OLTCs, switched capacitors [5][6][7][8][9].
For the incorporation of LVCs into the power flow, there are two common approaches: the error-feedback adjustment method and the automatic adjustment method [1,2,10]. In the error feedback methods, the control variables e.g., taps of SVR/OLTCs, states of capacitors and reactive power of DGs are updated after each iteration of power flow, based on the deviation of the regulated voltage from its target value. For instance, after the execution of each power flow iteration, the deviation of the voltage value of each LVC is computed and the respective tap or switching state is updated accordingly. The error feedback methods result in high oscillations of the power system states if the interactions between the LVCs are strong [2,3].
On the other hand, in automatic adjustment methods, the control variables are incorporated into the power flow equations and are solved simultaneously. This method is usually applied in the Newton-Raphson (NR) approach, in which the control variables are incorporated into the Jacobian matrix. Although the interactions between the control variables are considered in automatic adjustment methods and the oscillations are avoided, the discreteness of control variables is ignored, e.g., taps of transformers. The variables are discretized to the nearest discrete value after the power flow has converged, resulting in final voltage solutions, which may lie outside the permissible limits [2].
Another important challenge when incorporating LVCs into the power flow computation is the existence of multiple power flow solutions resulting from the different ways that a system state is reached. More specifically, LVCs are usually switched with different switching delays [1]. For example, an SVR placed near the substation may react faster in a voltage violation than an SVR placed at a farther distance. Furthermore, LVCs are set with different voltage bandwidths, inside which the state of LVC is unvaried. Due to the different reaction sequence of LVCs, the power flow could yield multiple steady state solutions, with all of them satisfying the voltage constraints. As a result, the actual switching sequence of LVCs should be considered according to their time delays and bandwidths, in order to compute the most probable power flow solution [1,2].
Another important challenge is the variety of operational modes of DGs. They may operate in constant voltage mode (modeled as PV nodes) [3] or in droop Q(V) mode according to the newly developed standard IEEE Std 1547-2018 from March 2018 [11][12][13][14][15][16][17][18]. In both cases, they interact with the other LVCs. Finally, the DGs may operate in different configurations e.g., 3-wire or 4-wire and may generate different voltage profiles e.g., balanced voltage or balanced current, which significantly affects the power flow results [14]. Therefore, precise modeling of the DGs is also required, which should simulate the actual operational characteristics of the DGs as well as the interactions with the other LVCs.

Literature Review
Two NR-based automatic adjustment methods are proposed in [15,16] by incorporating the LVCs into the Jacobian matrix. Another automatic adjustment method is proposed in [17], where the authors incorporate the LVC variables into the admittance matrix of the system forming a hybrid model, in which the LVC variables are solved simultaneously with the power flow.
An error-feedback method is proposed in [10], by adjusting the LVC variables after each iteration of power flow. However, a possible hunting oscillation may be observed between different adjustments due to the interaction of LVCs. To overcome the hunting, authors in [3] adopted sensitivity parameters that consider the interactions between the LVCs in single-phase formulation. In this way, the oscillations between consecutive power flow iterations are reduced, and the speed and robustness of the power flow solution are enhanced. All the aforementioned works ignore the time delays of LVCs, thus the actual switching sequence of the LVCs is not represented. As a result, the aforementioned power flow approaches do not yield the most probable power flow solution but rather a random one.
The time delays of LVCs were introduced for the first time into the power flow studies in [1]. The authors divided the LVCs into different delay groups. The LVCs with the fastest reaction belong to the first group, while the slowest ones belong to the last group. A similar division of LVCs into delay groups is applied in [18]. The method of [1] was improved in [2] by introducing sensitivity parameters to accelerate the convergence of the power flow. Authors in [2] claim that the power flow methods of [1,2] calculate precisely the power flow by estimating the most probable switching states of LVCs. However, it will be shown through dynamic simulations in Section 7 of this paper that none of the above-mentioned references calculate precisely the power flow since they ignore the actual switching sequence of LVCs. Furthermore, all methods described above (except [16,18]) neglect the network unbalances and the variety of DG operational modes investigating only single-phase networks.

Technical Contribution of This Paper
A lot of effort has been put so far by the researchers to study the optimal operation of LVCs in distribution networks e.g., [11,13,[19][20][21][22][23][24][25][26]. However, power flow computation in distribution networks with LVCs has not been extensively studied, and scope to improve remains. This manuscript contributes by proposing an error-feedback power flow method with the following distinct characteristics:  In contrast to the existing literature, the proposed method is a three-phase algorithm, which considers the network's unbalances. Moreover, the proposed power flow method accurately simulates various DG operational modes.  It calculates the actual switching sequence of LVCs based on their time delays. It is shown through dynamic simulation in Simulink that the proposed power flow method produces more accurate results than the methods of [1][2][3] which are considered the state of the art in the steady-state modeling of LVCs.  Three-phase self-and mutual-sensitivity parameters are derived relating the control variables (e.g., taps of OLTC and SVRs, capacitors, reactive power of DGs) with the regulated voltages of the network. The sensitivity parameters are then applied in the proposed algorithm to accelerate the convergence of the power flow, reducing significantly the computation time. This property is very important since power flow computation with LVCs is widely applied in real-time DMS applications, in which time is a crucial factor.  Unlike existing approaches, the proposed approach considers the weather and magnetic effects, which have a great influence on the impedance of overhead conductors [27,28], thus significantly affecting the power flow results and the switching state of LVCs.

Weather and Magnetic Dependent Power Flow
In this section, the weather and magnetic dependent power flow solver applied in this manuscript is briefly described. The algorithm utilizes the nonlinear thermal equations of IEEE Std 738-2012 for calculating the conductor temperature and the nonlinear impedance-temperature-current relation proposed in [27] for calculating the conductor impedance. The power flow is solved using the implicit ZBUS method [27] and matrix decomposition [29].

Thermal Modeling of Overhead Conductors
The temperature of an overhead conductor can be calculated from the nonlinear heatbalance equation of IEEE Std 738-2012 [30]. The nonlinear heat-balance equation is as follows: (1) where is the convective heat loss, is the radiative heat loss, is the heat gain due to solar radiation, and is the heat gain due to joule heating. The nonlinear heat-balance Equation (1) is solved with Newton's method to obtain the conductor temperature ( ), for any amount of current flow under any given weather condition experienced by the conductor (more details are given in [27,30]).

Modeling of the Magnetic Effects
IEEE Std 738-2012 [30] assumes a simplified linear resistance-temperature relationship that ignores the nonlinearities arising from the magnetic effects in the core of the conductors. Practically, both ACSR (aluminum conductor steel reinforced) and AAC (allaluminum conductor) 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, and thus, increasing the losses and impedance of conductor [27]. This effect is most pronounced in single-layer ACSR conductors (e.g., below 400 A [31]).
An algorithm to accurately calculate the impedance of overhead conductors based on the conductor temperature (calculated by Equation (1) in Section 2.1) and current (calculated by the power flow in Section 2.3) is described in [27].

Power Flow Solver Using Matrix Decomposition
The equation of implicit Z-Bus power flow is described by (2) below: where k is the iteration number, and are the vectors including the three-phase complex voltages and currents of bus i. is the self-or mutual-admittance matrix between bus i and j.
is the voltage of the slack bus. Equation (2) is iteratively solved until convergence. More details are provided in [27].
In order to avoid time-consuming repetitive factorizations of the ZBUS matrix of Equation (2), as a result of the impedance variation, the matrix decomposition technique is applied on the ZBUS matrix, as follows: where the branch-path incidence matrix is given in [29]. # for ∈ 1, … , includes the 3 × 3 impedance matrix of each three-phase line section. More details are provided in [29].

Modeling of Local Voltage Controllers
The modeling and operational modes of LVCs are briefly described below.

Switched Capacitors
Switched capacitors in distribution networks usually operate in voltage or reactive power control. In voltage control mode, the controller switches the capacitor ON when the measured voltage is less than a minimum value and OFF when it is more than a maximum value [32]. In reactive power mode, the state of capacitor depends on the reactive power flowing through the controlled line.

Step Voltage Regulators
A three-phase SVR model is proposed in [33], the equivalent of which is depicted in Figure 1. It is a 3-bus equivalent that simulates the taps as current sources, enabling the construction of sensitivity parameters relating the tap variables with the voltages, as explained in the next section.

On-Load Tap Changing Transformers
A three-phase OLTC transformer model is presented in [34] and its 2-bus equivalent circuit is shown in Figure 1. The taps are simulated as current sources in a similar sense as in SVRs.

Distributed Generators
The steady state modeling of DGs is usually categorized depending on the generated voltage/current profile as well as the power profile, as follows: 1. Voltage and current profile: Electronically coupled DGs can generate balanced phase-to-phase voltage, balanced phase-to-neutral voltage, or balanced current. On the opposite, synchronous generators present nonzero finite negative-and zero-sequence admittances, and therefore, they generate unbalanced voltages and currents [14]. 2. Power profile: DGs can operate in constant PQ, constant PV (conventional PV bus modeling), or constant P-Q(V) mode [12,13]. The first case does not pose any challenge to the power flow computation since the DG is simply modeled as a negative constant power load. However, DGs with controllable reactive power may interact with the other LVCs increasing the number of power flow iterations or even leading to divergence.

Time Delays of LVCs
LVCs in a network are regulated to react with pre-defined time delays. This occurs in order to coordinate the operation of LVCs and to avoid unnecessary switchings, due to temporary voltage violations. Switching capacitors, for example, are set intentionally with a time delay such that the capacitor switching is activated after a pre-defined time delay from the instant of voltage (or reactive power) violation. In the following, the delay of capacitor is denoted as Τ . Similarly, SVR and OLTC step the voltage up or down when a voltage violation occurs beyond a pre-defined bandwidth, considering a pre-specified time delay. For instance, once a voltage violation at the supervised bus is noticed, the counter of LVC starts counting. The first step-up/down switching action is executed as soon as the counter reaches the intentional time delay. If the voltage remains outside the bandwidth, additional step-up/down actions are executed with a mechanical delay. The intentional and mechanical delays of SVR and OLTC are denoted here as Τ DGs are assumed to react instantaneously to any load/generation variation [11,20,22,23].
Switching capacitors, SVRs, and OLTCs can operate in two different ways: (a) their three phases are independently controlled, namely every phase has its own LVC, or (b) The three phases are controlled simultaneously by the same LVC [35][36][37].

Presentation of the Proposed Algorithm
Let us firstly assume a network with local controllers, which control the switching capacitors, SVRs and OLTCs. Let us further assume that N of them control an SVR (one or more phases depending on whether the three phases are independently or simultaneously controlled), N of them control an OLTC, while N of them control a capacitor so that N ⊆ N, N ⊆ N and N ⊆ N as well as N N N N.

INITIALIZATION:
Initialize the taps of SVRs and OLTCs as well as the capacitor states, and set T → ∞ ∀i ∈ N. Specify the voltage bandwidth ( STEP 1: Compute the power flow considering the DGs. In this step, the power flow is executed until convergence considering the operational mode of DGs (refer to Section 3.4). STEP 2: Calculate T for each controller ∈ N as follows: The function , in Equation (4) returns the lowest value between x and y. and are the calculated (from step 1) and desired voltage magnitude of local controller , respectively.
If the voltage is outside the limits, the priority indicator T is updated taking the lowest value between the intentional delay of the local controller and the current value of T . If the voltage is inside limits, T is set to infinite.

STEP 3:
Find the controller with the minimum priority indicator such that T T ∀ ∈ . Then, for each controller ∈ , set T T T . As shown, the priority indicator with the lowest value is subtracted from the other priority indicators. Consequently, controller with the minimum T becomes 0 and it will undergo switching action. It should be noted that if one or more controllers have the same priority indicator value with controller , their priority indicator also becomes 0 in this step, and hence, they will also undergo switching actions simultaneously.
This step prioritizes the switching actions based on their time settings. A zero T value means that the counter of LVC j has reached its time delay, and therefore, a switching action should be undertaken from j th LVC.

STEP 4:
In this step, controller with the minimum T value, undergoes switching action according to Algorithm 1.

Algorithm 1. Switching Logic 1
In Algorithm 1, for instance, if the LVC j is an SVR or OLTC and its regulated voltage is below bandwidth ( ) as well as its tap setting is below the maximum value ( ), then the tap is increased by one ( 1) in order for the voltage to move toward the bandwidth. If the LVC j is a capacitor, and its regulated voltage is below bandwidth as well as the capacitor is OFF ( ), then the capacitor is turned ON to increase the voltage of the regulated bus toward the bandwidth.

STEP 5:
In this step, T of controller j is updated based on its mechanical delay, as follows: In this step, the time setting of LVC j is set equal to its mechanical delay. This step represents the real operation of LVCs, based on which, after the intentional delay of the first switching action, the LVC is subsequently switched with a mechanical delay.

STEP 6:
If at least one switching action took place in STEP 4, go to STEP 1 and run the power flow. If no switching action was executed, go to STEP 2. The algorithm terminates when no further switching actions are required and error tolerance is met.

Three-Phase Sensitivity Parameters of LVCs
The algorithm presented in the previous section is accurate but requires a very large number of power flow iterations since it executes the power flow after each individual switching action of capacitors, SVRs, and OLTCs. To accelerate the convergence, a novel sensitivity-based algorithm is proposed here.
To derive the three-phase sensitivity parameters of LVCs, let us first assume the network of Figure 1. It includes an OLTC connected between the buses h-w [34], three singlephase capacitors connected to the three phases of bus q, and a DG connected to bus g. Furthermore, an SVR is also connected through the 3-bus equivalent circuit of [33]. Figure 1 describes the positive sequence current that is injected by the DG at the three phases. Similarly, 0 and 2 describe the zero-and negative-sequence currents injected into the three phases. It is pointed out that the OLTC model of [34], represents an OLTC transformer that connects a 3-wire MV and a 4-wire multi-grounded LV network. Nevertheless, in this paper, we investigate distribution MV networks, which are usually 3-wire. Therefore, in Figure 1, the model of [34] is simplified and the current sources of neutral and grounding conductors are ignored.
Let us define the vector of LVC control variables ( ) and voltages as: where:  includes the tap variables of each phase of the OLTC transformer. The vector directly regulates the vector | | | | | | that includes the magnitude of the three-phase voltages of bus w. The variation of the vector is related to the variation of the through a sensitivity matrix ( , ), as shown in Equation (8).
In Equation (9), the sensitivity matrix , is of dimension 14 × 14 corresponding to the schematic presented in Figure 1. The sensitivity matrix , relates the variations of the control variables to the voltage's variations. This relationship is utilized in the proposed novel power flow algorithm to accelerate the convergence, as presented in the following section. Due to space constraint, a detailed derivation of the sensitivity parameters in Equation (9) is presented in a supplementary document in [38]. It is important to clarify that the elements of , matrix are derived, assuming that the control variables , , are continuous. However, this assumption does not affect the accuracy of the results since the distinct nature of the , , is accurately considered in the sensitivity-based algorithm of Section 6 and the discrete variables are updated in a discrete manner.

A Novel Sensitivity-Based Algorithm with Accelerated Convergence
The sensitivity parameters, which are presented in the previous section, are utilized here to derive a novel sensitivity-based three-phase weather-dependent power flow algorithm for distribution networks with LVCs. The algorithm has accelerated convergence due to the utilization of the sensitivity parameters. It is shown through simulation results that the sensitivity-based algorithm presents identical results with the generic algorithm of Section 4, while the iterations and computational effort are significantly reduced.
The algorithm consists of the following steps: INITIALIZATION: Similar to the initialization in Section 4. Equation (10). In Equation (10), and are the deviation of and from their set reference voltages, respectively. The notation denotes an estimate value.
Due to the coupling of LVCs, the variations and calculated in (10) cause voltage variations at the other LVC buses. By utilizing the sensitivity parameters and the results of Equation (10), the predicted voltage variations ( , , ) of the LVCs are calculated, as shown in Equation (11) The new predicted voltages , , of LVCs are then calculated in Equation (12): (12) where , , denote the voltages of OLTC, SVRs and capacitors, respectively, as calculated in the previous iteration.

STEP 5:
In this step, based on the predicted voltages (from STEP 4), we perform switching action at the LVC with the minimum T , as presented in Algorithm 2.

Algorithm 2. Switching Logic 2
If ∈ or ∈ then It should be noted that the Algorithm 2 is similar to Algorithm 1 of Section 4. However, Algorithm 1 uses the voltage value of controller ( ), as calculated from the power flow. Therefore, a large number of iterations are required since the power flow is executed after each individual switching action. On the other hand, Algorithm 2 uses the predicted voltage value of controller ( ), which is calculated from Equation (12), without requiring the execution of a power flow.

STEP 6:
If a switching action took place in STEP 5, it would make the DG voltages deviate from their reference values. Therefore, further DG action is required to maintain DG reference voltages. This is achieved by correcting again the DG control variables and , as follows in Equation (13). Note that Equation (13) is derived from Equation (A1) presented in the Appendix A.
Finally, the DG control variables are updated as shown in Equation (14), using the result of (13).
In case the positive-sequence reactive power ( ) of DG exceeds its upper or lower limit, it is set equal to the respective limit, so that the limits are not exceeded. The negative-and zero-sequence currents ( ) are not restricted since they are a small fraction of the positive sequence current, and therefore, the current of DG is mainly determined by the positive-sequence component.  (14). Go to Step 2 until no further switching actions take place. A flowchart of the proposed method is depicted in Figure 2.
With the proposed sensitivity-based approach, a full power flow is required to be executed only at the initialization step and not after each individual switching action as in Section 4. Therefore, the total number of iterations in the proposed sensitivity-based algorithm has been significantly reduced. The distinct features of the proposed algorithm are further highlighted via simulation and validation in the following sections.
Finally, it is pointed out that in this paper, we have considered a temperature-dependent power flow algorithm for distribution networks with bare overhead conductors, using the IEEE Std 738-2012 [30]. However, the proposed approach can be also applied in distribution networks with underground cables. In this case, the model that was recently proposed in [39] can be applied, instead of the IEEE Std 738-2012, for the estimation of conductor temperatures.

Validation and Performance of the Proposed Sensitivity-Based Algorithm
First and foremost, the proposed algorithm is validated against dynamic simulation in MATLAB ® Simulink using an 8-bus balanced and a 7-bus unbalanced network. In addition, the proposed approach is compared against the power flow methods of [1] and [3], with respect to result accuracy and convergence speed. It is noted that the method proposed in [2] presents almost identical results with the method of [1]. Therefore, the method of [2] is not simulated in this paper due to its strong similarity with [1] with respect to the produced results.
All algorithms were coded and implemented in MATLAB ® . It should be noted that for the purpose of validation against the dynamic simulation, the weather-dependent impacts were neglected in this section as traditional algorithms and Simulink do not consider weather-dependent modeling.
where , , , are the positive sequence reactive power, the reference voltage, the droop gain, and the positive sequence voltage of DG i, respectively. The topology of the network is similar to the one investigated in [21]. Data about the network, the controllers of SVRs, and the DG are provided in Table 1. The delays of the controllers are set based on their distance from the substation. The slack bus in Figure 3 is assumed to be a substation. The SVR near the substation has the fastest reaction time and as the distance from the substation increases, the reaction time increases [23,40]. All SVRs and OLTCs of the paper are considered to have, in total, 33 distinct taps (e.g., −16, …0, …+16) and can boost or reduce the voltage by ±10%. So, each tap variation changes the voltage by ±0.0625%. All SVRs are in wye configuration and are initialized to the 0 position for the purposes of simulations. Each phase of an SVR is modeled with its own local controller, which is independently controlled. The DG at bus 8 operates in mode and generates balanced phase-to-neutral voltages [14]. The network supplies four balanced three-phase loads. All the loads are modeled as constant impedance loads, as shown in Table 2.
The output phase-to-neutral voltages of the three SVRs as well as the reactive power of DG are depicted in Figure 4. It should be noted that since all the loads are balanced, the three phases of the SVRs have equal voltages and undergo similar tap changes. Initially, all the voltages are below their bandwidths, and therefore, the SVRs undertake step-up actions with the intentional and mechanical delay. Initially, due to the low voltage at bus 8, the DG generates reactive power. As long as the SVRs boost their voltage, the reactive power of DG is reduced based on the droop curve (Equation (15)). The voltage of all SVRs is stabilized inside the bandwidth after 30 s. After the stabilization of SVRs, an overvoltage occurs in bus 8 due to the active power generation of DG, and thus, the DG consumes reactive power according to Equation (15).
The tap positions and switching sequence for all SVRs of the 8-bus network obtained via the dynamic simulation of Simulink are presented in Figure 5. SVR 1 executes the first switching action after an intentional time delay (10 s), followed by subsequent switching actions with a mechanical time delay (2 s), until the voltage lies within the bandwidth. Similarly, for SVR 2 and 3. Since SVR 1 has the lowest intentional and mechanical delay, it undergoes the highest number of tap changes. The final tap positions obtained via dynamic simulation for SVR 1, SVR 2 and SVR 3 is 11, 3, and 1, respectively.
In Figures 6-9, the results of the SVR tap change versus iteration number are presented for the proposed algorithm as well as for the algorithms of [1,3]. Figure 6 presents the tap change profile of all the SVRs versus the iteration number for the proposed algorithm without consideration of sensitivity parameters (refer to Section 4). In Figure 7, the tap change profile versus iteration number is presented for the proposed algorithm with the consideration of the sensitivity parameters (refer to Section 6). The proposed algorithm, for both with and without sensitivity parameter, considers the actual reaction delays of LVCs. However, as observed in Figure 7, the consideration of sensitivity parameters yields accelerated convergence. It is reminded that the proposed algorithm, both with and without sensitivity parameters, considers the same switching algorithm, thus, it presents identical results in both cases. Their main difference is that the usage of sensitivity parameters makes the execution of a complete power flow after each switching action unnecessary, reducing the total iteration number required for the final power flow solution.   It is observed that the proposed power flow algorithm yields correct final tap positions when compared to the dynamic simulation, both with and without sensitivities. This is due to the consideration of the actual switching sequence of LVCs, as presented in Section 4. On the other hand, the LVC's states estimated by the power flow algorithms of [1] and [3] do not conform to the dynamic simulation, as shown in Figures 8 and 9.
Algorithm [1] is only able to correctly estimate the final tap position of SVR 2 as shown in Figure 8. It is reminded that the authors in [1] divide the LVCs into delay groups based on their reaction delays. In this example, SVR 1, SVR 2 and SVR 3 belong to the first, second and third delay group, respectively. Initially, the SVR of the first delay group reacts, by varying its taps until its voltage lies inside the bandwidth. Subsequently, the SVR of the second delay group undertakes switching actions, and so on until all SVR voltages lie inside their bandwidths. However, this process does not outline exactly the actual sequence of switching of LVCs, because in real distribution networks, the LVCs are not separated into delay groups and their reaction is not executed in the way that the algorithm of [1] describes. For instance, in the examined 8-bus network, the actual switching sequence is depicted in Figure 5 using the dynamic simulation of Simulink, which is quite different than the switching sequence predicted in [1].
The algorithm of [3] reaches final tap position in only a few iterations (Figure 9), because the algorithm updates all the LVC states simultaneously, and thus, reduces the total power flow iterations. It is reminded that algorithm [3] neglects completely the switching delays of LVCs, and updates the control variables of LVCs (e.g., taps of SVRs/OLTCs and reactive power of DGs) simultaneously, using sensitivity parameters. Due to the discrete nature of taps, after the convergence of the power flow, algorithm [1] rounds the taps to the nearest tap position.
The final voltages of the network for the investigated methods and Simulink are provided in Table 3 at the end of the paper. The proposed approach with and without sensitivities present almost identical results with those of Simulink, in contrast to the methods of [1,3] that show significant deviations.
In Figure 10, the total three-phase reactive power of the DG connected at bus 8 is presented. It operates in inductive mode according to droop Equation (15), to mitigate the voltage rise caused by the high amount of generated active power [14]. It is observed that the proposed algorithm indicates a consumption of 542 kVar, which is exactly the same as in the dynamic simulation. References [1,3] yield different results. The deviation of the reactive power in algorithms [1,3] is caused due to the inaccurate state estimation of SVRs, which inevitably leads to an imprecise reactive power calculation.

(B) 7-Bus Unbalanced Network
An unbalanced 7-bus network consisting of an OLTC, a voltage-controlled capacitor, an SVR, and a DG is considered for further simulation and validation, as shown in Figure  11. This network was selected for simulation since it includes all kinds of LVCs (OLTC, SVR, Capacitor, DG), it is unbalanced and also its simple topology facilitates the comparison and interpretation of the simulation results. Moreover, its small size allows the execution of simulations in the time domain environment of MATLAB ® Simulink. The DG operates in droop control mode (see Equation (15)). Data about the network, the controller of OLTC, SVR, capacitors and DG are presented in Table 4. The SVR is connected in wye configuration, while OLTC in Yg-Yg connection [41]. Each phase of OLTC, SVR, and capacitor has its own local controller and is independently controlled. The DG generates balanced phase-to-neutral voltages [14]. The network supplies three balanced and one unbalanced three-phase constant impedance load, as shown in Table 5.   Table 6 summarizes the calculated LVC states of each phase in the 7-Bus network using dynamic simulation, the proposed algorithm without and with sensitivity parameters, as well as the methods of [1,3]. It is observed that the proposed algorithm presents identical results with those of Simulink, confirming its accuracy. On the other hand, the calculated states of the approaches in [1,3] deviate from those of Simulink. Moreover, Table 7 at the end of the manuscript depicts indicatively the voltages of each phase for the last three buses of the network for the investigated approaches. As shown, the proposed method with and without sensitivity parameter yields near identical results with those of Simulink, while the other investigated power flow methods present significant deviations.

Algorithm [3]
(11, 11, 10) (OFF, ON, OFF) (3, 3, 2) 625 kvar (inductive) Table 7. Voltage profile (of the last three buses) of the investigated methods for the 7-Bus unbalanced network. Finally, the number of iterations required by the algorithms to converge with an accuracy of 10 −4 pu are presented in Tables 8 and 9, for the 8-Bus and 7-Bus networks, respectively. As shown, the algorithm of [3] presents the fastest convergence since it calculates all LVC's state simultaneously ignoring their reaction delays. The proposed algorithm without the sensitivity parameters has the slowest convergence due to the successive power flow executions after each switching action. The proposed method with sensitivity parameters combines the high accuracy with the fast convergence, as no other power flow method so far.

Influence of Weather on the LVC States and Power Flow
We conduct a case study on the large IEEE 8500-Node [42] network to compare the proposed weather-dependent algorithms (with and without sensitivities) against the conventional methods of [1,3]. Moreover, the accuracy of the proposed sensitivity-based algorithm in a large-scale network is investigated. Figure 12 depicts the IEEE 8500-Node Network, which has been modified by including 4 DGs operating in several modes. The network originally includes 4 SVRs and 4 capacitors [42], as shown in Figure 12.

(A) Network description
The 4 additional DGs are connected to the buses 100, 350, 835, 1600. Data about the power profile and voltage/current profile of DGs are provided in Table 10, while data about the reference voltages, the droop gains, the active powers, etc. are given in Table 11. For DG 1, which is a synchronous generator (SG), the negative-and zero-sequence impedance are 0.2 0.2. DG 1 and DG 4 operate in constant voltage mode (treated as PV buses) and generate constant active power and positive-sequence voltage. DG 2 and DG 3 operate in droop control mode generating balanced currents (refer to Section 3.4 for more details about the operational modes of DGs).
Data of the LVCs are provided in Table 11. Each phase of the LVCs is independently controlled. The capacitors are voltage controlled and the SVRs are in wye configuration. The time delays of SVRs and capacitors are set based on their distance from the substation [23,40], while DGs react instantaneously. LVCs near substation have faster reaction times as shown in Table 11. Each controller regulates the voltage of the connection point of LVC and no line drop compensator or any other remote voltage control is applied.
Finally, all the lines of the network were replaced with the Penguin ACSR [43]. It is a single-layer conductor with a cross-sectional area 125.1 mm composed of a 6 Aluminum and 1 steel wire. Penguin is the largest single-layer ACSR conductor and can successfully withstand the full load of the network in both investigated environmental conditions, without thermal violation. This modification was necessary in order to simulate the influence of weather and magnetic effects into the power flow results. With the original lines of IEEE 8500-node network, this would not be possible since all lines consist of constant impedances with unknown conductor-specific details.
In Figures 13 and 14, the resistance and self-reactance of Penguin ACSR conductor are presented, respectively, as a function of conductor temperature and current. The methodology to calculate the resistance and self-reactance of Figures 13 and 14 is presented in [27]. It is observed that the resistance of Penguin ACSR is strongly related to both current and conductor temperature (conductor temperature is a function of both current and weather conditions [27]). Nevertheless, the impact of weather is significant as demonstrated in the following subsection.
The case study involves the simulation of two scenarios with two different weather conditions (summer and winter), as presented in Table 12. The first one represents a typical summer day in Thessaloniki, Greece, while the second represents a typical winter night in the same region [44]. The directions of the lines were considered exactly as they are given in Figure 12, while the wind direction was considered parallel to the x-axis of Figure 12. The original network loads were considered for all buses [45].

(B) Comparison between the investigated algorithms
Simulations were executed in a generic PC with an Intel Core i7, 3.4 GHz CPU and 16 GB DDR3 RAM. As mentioned before, all investigated algorithms were coded in MATLAB ® . More specifically, the following algorithms are compared:  Proposed weather-dependent power flow algorithm without sensitivities.  Proposed weather-dependent power flow algorithm with sensitivities.  Algorithm [1] with fixed impedances.  Algorithm [3] with fixed impedances.
Usually, utilities apply traditional power flow algorithms assuming line impedances of fixed values, which are seasonally updated e.g., summer, winter to account for the weather variations. In this study, we assume fixed-line impedances for algorithm [1] and [3], with different values in summer and winter period, as follows:  In Summer: Due to the large size of the network, modelling and dynamic simulations in Simulink are not possible here. In all approaches, tap position 2 was assumed as initial for all SVRs, while all the capacitors are initially considered OFF.  The tap positions and switching sequence for SVR 3 in the modified IEEE 8500-node network are presented in Figure 15. Figure 15 presents indicatively the tap positions and switching sequence for SVR 3 (all phases) obtained via the proposed algorithm, with and without the sensitivity parameters, for the summer scenario. As observed, the final tap positions are in agreement for the proposed algorithm with and without the use of sensitivity parameters. This is because both cases use similar algorithms to estimate the exact switching sequence of SVRs and capacitors (see Section 4 and 6). It is also observed that the use of sensitivity parameters yields a much faster solution.
A summary of the final tap positions of LVCs, for all examined algorithms, in winter and summer conditions, is presented in Table 13, while the reactive power of DGs is quoted in Table 14. As observed, for the same loading conditions, algorithms [1] and [3] derive different power flow results compared with the proposed approach, due to the neglect of weather and magnetic impacts as well as the actual switching sequence of LVCs. It should be noted that the proposed algorithm produces the same power flow results for the two cases, with and without sensitivity parameters. The difference observed between the two cases was in the convergence characteristics and the evolution of the solution states. The numbers of iterations and the total computation time required by the algorithms to converge with an accuracy of 10 −4 pu are presented in Table 15, for the modified IEEE 8500-node network. As shown, the proposed algorithm with the usage of sensitivity parameters appears to be the best power flow tool for networks with LCVs, since it combines high accuracy with fast convergence and low computation time. Although algorithm [3] is by far the fastest, the results are not accurate for the aforementioned reasons.

(C) Accuracy of the Sensitivity-Based Algorithm
In this sub-section, we investigate how accurate the proposed sensitivity-based algorithm is, compared with the proposed algorithm without sensitivities, when applied in the large IEEE 8500-node network. Six scenarios were simulated for this purpose in the IEEE 8500-node network, explained below. In all scenarios, the winter condition was assumed for calculating the line impedances. The results of LVC's state are presented in Table 16 for the proposed method without and with the use of sensitivity. As shown, the algorithms produce identical results in the first four and 6th scenarios. In most of the simulations presented in this manuscript, there has been full agreement between the proposed algorithm without and with sensitivities. However, a slight difference is observed in the state estimation of SVR 4 of scenario 5 (highlighted in bold in Table 16), where the two algorithms output (3, 2, 0) and (3, 1, 1), respectively. This difference is attributed to the fact that the usage of sensitivities involves predictor and corrector steps, which are approximations. Nevertheless, this difference is negligible in comparison to the major reduction in the number of iterations required to obtain the power flow solution, when sensitivities are adopted, as observed in the last column of Table 16.

(D) Weather Data Acquisition
The real-time measurement of local weather parameters is feasible today since the number of automatic weather stations (AWSs) has been considerably increased worldwide [44]. For instance, in Greece, National Observatory of Athens (NOA) has installed, over the last 15 years, more than 430 automatic weather stations over Greece (see Figure  16). All weather stations record (every 10 min) the 10-min average values of wind speed, wind direction, ambient temperature and solar radiation, namely all the weather parameters required by the proposed approach. The data are provided in a real-time manner to potential users, through web-based platforms [44], and therefore, they are currently very easily accessible. Figure 16. Geographical distribution of the network of automatic weather stations over Greece [44]. The red points indicate the position of weather stations.

Conclusions
This paper presents a sensitivity-based three-phase weather-dependent power flow algorithm for simulating distribution networks with LVCs. The proposed algorithm has four distinct characteristics: (a) it considers the three-phase unbalanced nature of distribution systems, (b) the operating state of LVCs is calculated using sensitivity parameters accelerating the convergence speed of the algorithm, (c) it considers the exact switching sequence of LVCs based on their reaction time delays, and (d) the influence of weather variations on the power flow is also taken into consideration. The proposed algorithm has high accuracy and low computational complexity and is therefore can constitute an important tool in offline and real-time applications of distribution systems. The proposed approach was compared against some state-of-the-art power flow methods using a 7-bus, an 8-bus and the IEEE 8500-node network. The results indicate that the existing algorithms can produce inaccurate power flow results and LVC's states, due to the imprecise consideration of LVC's switching sequence. Finally, it was deduced through simulations, that for a more precise estimation of LVC's states and power flow results, the weather variations should be considered in the power flow.

Conflicts of Interest:
The authors declare no conflict of interest.