A Sensitivity-Based Three-Phase Weather-Dependent Power Flow Algorithm for Networks with Local Controllers—Part I: Algorithm Development

—Local voltage controllers (LVCs) are important components of a modern distribution system 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. This Part I presents the theoretical development of the proposed algorithm, which 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, which accelerates the convergence speed of the algorithm, c) it considers the precise switching sequence of LVCs based on their reaction time delays, and d) the nonlinear influence of weather variations in the power flow is also taken into consideration. Simulations and validation results presented in Part II 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.


I. INTRODUCTION
OCAL OLTAGE controllers (LVCs) are important components in modern distribution networks. An LVC is a physical apparatus such as on-load tap changer transformer (OLTC), step voltage regulator (SVR), switched capacitor, distributed generator (DG), which regulates the voltage of a specific bus within permissible limits [1] [2]. For safe, reliable, and efficient operation, distribution management systems (DMS) regularly carry out optimization functions such as Volt/Var control (VVC) and optimal feeder reconfiguration (OFR), which require precise power flow solutions with accurate modeling of LVCs [3]. Therefore, accurate modeling of LVCs in power flow studies is essential to understand their interactions in the distribution network and have appropriate understanding of the power system states.
This research is co-financed by Greece and the European Union (European Social Fund-ESF) through the Operational Programme «Human Resources Development, Education and Lifelong Learning» in the context of the project "Strengthening Human Resources Research Potential via Doctorate Research -2nd Cycle" (MIS-5000432), implemented by the State Scholarships Foundation (ΙΚΥ).

A. Challenges of Incorporating LVCs in Power Flow
Two common approaches for incorporating LVCs into the power flow are the error-feedback adjustment method and the automatic adjustment method [1] [2] [4]. 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 is the existence of multiple power flow solutions resulting from the different ways that a system state is reached. More specifically, LVCs do not react instantaneously but they are assigned to react with different reaction delays [1]. For instance, an SVR placed near the substation may react faster after a voltage violation occurs than an SVR placed in a farther distance. Moreover, LVCs are assigned different voltage bandwidths, inside which the LVC state remains unchanged. Depending on the reaction sequence of LVCs, the power flow could yield multiple solutions, with all of them satisfying the voltage constraints. Therefore, in order to find the most probable power flow solution in networks with LVCs, the actual reaction sequence of the LVCs should be considered based on their time delays, voltage settings and bandwidths [1].
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 [5]- [7]. 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 [8]. 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.

B. Literature Review
Two NR-based automatic adjustment methods are proposed in [9] and [10] by incorporating the LVCs into the Jacobian matrix. Another automatic adjustment method is proposed in [11], 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 [4], 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 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 [12]. 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] and [2] produce the most probable power flow solution by accurately simulating the actual switching sequence of LVCs. It is pointed out that all methods described above (except [10] and [12]) neglect the network unbalances and the variety of DG operational modes investigating only single-phase networks.

C. 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. [5], [7], [13]- [20]. However, the field of power flow in distribution networks with LVCs has not been extensively studied, and a lot of work is still needed to improve the accuracy and computation time of power flow. This two-part paper aims to contribute 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 (presented in Part II) that the proposed power flow method produces more accurate results than the methods of [1] and [2] which are considered the state of the art in steady-state modeling of LVCs [2].
• 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. It is shown through simulation results in Part II that fewer iterations are needed for the proposed algorithm to converge, thus significantly reducing the computation time. This property is very important since power flow with LVCs is widely applied in real-time DMS applications, in which the time is a crucial factor.
• Unlike the existing approaches, the proposed approach considers the weather and magnetic effects, which have a great influence on the impedance of overhead conductors [21] [22], thus significantly affecting the power flow results and the switching state of LVCs. The matrix decomposition technique [23] is applied to avoid continuous re-factorization of YBUS matrix, due to the impedance variation, thereby reducing the computational burden of the power flow.

D. Paper Structure
The structure of the paper is as follows: Section II presents a theoretical analysis of the applied weather-dependent power flow algorithm. Section III presents the simulation models of the LVCs, while Section IV proposes a generic and accurate, but computationally expensive, algorithm for calculating the power flow in networks containing LVCs. Section V introduces the three-phase sensitivity parameters, based on which the novel sensitivity-based power flow algorithm is introduced in section VI. Finally, Section VII concludes this Part I.

II. WEATHER AND MAGNETIC DEPENDENT POWER FLOW
In this section, the weather and magnetic dependent power flow solver applied in this paper is described. The algorithm utilizes the nonlinear thermal equations of IEEE Std 738-2012 [24] for calculating the conductor temperature and the nonlinear impedance-temperature-current relation proposed in [21] for calculating the conductor impedance. The power flow is solved using the implicit ZBUS method [21] and matrix decomposition [23].

A. Thermal Modeling of Overhead Conductors [24]
The temperature of an overhead conductor can be calculated from the nonlinear heat-balance equation of IEEE Std 738-2012 [24]. The nonlinear heat-balance equation is as follows: 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. IEEE Std 738-2012 provides several equations for calculating the convective heat loss [24]. The convective heat loss ( ) is caused either by forced convection or natural convection depending on the wind speed ( ). Two alternate equations are recommended by IEEE Std 738-2012 for the calculation of forced convection, as shown in Eq. (2) and (3), respectively.
The equation for natural convection ( ) is as follows: In Eqs.
(2)-(4), is the thermal conductivity of air, is the wind direction factor, 0 is the conductor diameter, is the air density, is the dynamic viscosity of air, is the conductor temperature, and is the ambient temperature.
In order to obtain the value of , Eq.
The solar heat gain ( ) is given in Eq. (6), and concerns the energy transferred from the sun to the conductor through solar radiation. In Eq. (6), is the solar absorptivity, and is the solar irradiance.
The Joule loss ( ) is given in Eq. (7). It refers to the active power losses generated in the conductor due to the current ( ) flow through the conductor resistance ( ).
The nonlinear heat-balance Eq. (1) can be solved to obtain the conductor temperature ( ) 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 ( ) is defined as follows: The heat-balance function in Eq. (8) is then solved using Newton's method to obtain the conductor temperature ( ).
After the conductor temperature 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 subsection.

B. Modeling of the Magnetic Effects [21]
IEEE Std 738-2012 [24] 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 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 [21].
The aforementioned effect is most pronounced in singlelayer 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. Since all ACSR conductors with a current capacity less than 400A are single-layer conductors [25], these are widely 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.
An algorithm to accurately calculate the impedance of overhead conductors considering the magnetic properties of the steel core is described in [21,Algorithm 1]. The algorithm is executed offline for each type of conductor and outputs two 2-dimensional matrices, which include the values of resistance and reactance for a wide range of conductor temperature and current.
Using this algorithm, the resistance and self-reactance of the Rabbit ACSR and Ant AAC conductors are depicted in Fig. 1 and Fig. 2, respectively, as a function of the current and conductor temperature [21]. Looking at the figures, the following observations are derived: • The resistance and reactance of ACSR varies nonlinearly with the temperature and current. • The resistance of AAC varies linearly with the temperature, while the self-reactance is constant. The power flow solver applied in this paper, takes into consideration the weather and magnetic effects on the line impedances, thus generating more precise results than the conventional power flow algorithms, e.g. [1]- [4], [9]- [12].  C. Power Flow Solver using Matrix Decomposition [23] The three-phase configuration of a distribution network with m+1 buses is shown in Fig. 3. Let us define the current and voltage vectors of bus , as follows: where and denote the complex voltage and current of bus i at phase r = {a, b, c}, respectively, as shown in Fig. 3.
The load vectors can be expressed as a function of the voltage vectors according to Eq. (11): where coresponds to the slack bus voltage. Subsequently, for the power flow solution, we remove the first three rows of Eq. (11) that correspond to the slack bus and Eq. (12) is obtained.
In order to avoid time-consuming repetitive factorizations of the ZBUS matrix of Eq. (13), as a result of the impedance variation, the matrix decomposition technique is applied, as follows: where the branch-path incidence matrix K is given in [23]. # for i={1 , ... , m} includes the 3x3 impedance matrix of each three-phase line section [23].

III. MODELING OF LOCAL VOLTAGE CONTROLLERS
The modeling and operational modes of LVCs are briefly described below.

A. 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 [26]. In reactive power mode, the capacitor is switched depending on the reactive power value and direction through the supervised line. In power flow formulation, capacitors are simply modeled as constant impedance loads with = ( ) −1 .

B. Step Voltage Regulators
Step voltage regulators consist of autotransformers connected in several configurations e.g., open-delta, closeddelta, wye [27]. A three-phase SVR model is proposed in [27] and it is depicted in Fig. 4.
The equations of current sources are provided in [27, Table  III], while the admittances , in [27, Table II]. The tap variables are included in the current sources, thus avoiding the continuous factorization of YBUS matrix after each change of taps. Furthermore, the inclusion of taps into the current sources enables the construction of sensitivity parameters relating the tap variables with the voltages, as explained in the next section.

C. On-Load Tap Changing Transformers
A three-phase OLTC transformer model is presented in [28] and is depicted in Fig. 5

D. 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. Inverter-based DGs can generate balanced phase-to-phase voltage, balanced phase-to-neutral voltage, or balanced current. On the contrary, the synchronous generator-based DGs present nonzero finite negative-and zero-sequence admittances, thus they usually generate unbalanced voltages and currents [8].
2. Power profile. DGs can operate in constant PQ, constant PV (conventional PV bus modeling), or constant P-Q(V) mode [6] [7]. The first case does not pose any challenge to the power flow since the DG is simply treated 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.

A. Time Delays of LVCs
LVCs in a network are configured to react with prespecified time delays. This is done in order to coordinate the switching actions of LVCs and to avoid unnecessary switchings, due to the temporary voltage variations. Switching capacitors, for example, are set intentionally with a time delay such that the capacitor switching is activated after a pre-specified 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-specified bandwidth, considering a pre-specified time delay. More specifically, as soon as a voltage violation is detected, the intentional time delay is counted before the first stepup/down switching action. 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 present a very fast reaction time compared to the intentional or mechanical delays of the other LVCs. Thus, it is assumed that they react instantaneously to any load/generation variation [5], [14], [16], [17].
Switching capacitors, SVRs, and OLTCs can operate in two different ways: a) their three-phases are independently controlled; thus, each device has three local controllers acting independently on each phase, or b) each device has one local controller, which controls the three phases simultaneously [29]- [31].

B. 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 of them control an SVR (one or more phases depending on whether the three phases are independently or simultaneously controlled), of them control an OLTC, while of them control a capacitor so that ⊆ , ⊆ and ⊆ as well as + + = . The proposed algorithm defines a priority indicator ( ) for each controller ∈ . The switching sequence of the local controllers is executed based on this indicator. More specifically, only one switching action is executed at each time from the local controller with the lowest indicator value.
The steps of the proposed algorithm are described below: INITIALIZATION: Initialize the taps of SVRs and OLTCs as well as the capacitor states, and set ( ) → ∞ ∀ ∈ . Specify the voltage bandwidth ( ), the time delays STEP 1: Run the power flow considering the DGs. In this step, the power flow is executed until convergence considering the operational mode of DGs (refer section III.D).

STEP 2:
Calculate ( ) for each controller ∈ as follows: The function min(x, y) in Eq. (15) returns the lowest value between x and y. and ( ) are the calculated and desired voltage magnitude of local controller , respectively.
If the voltage is outside the limits, the priority indicator ( ) is updated taking the lowest value between the intentional delay of the local controller and the current value of ( ) . If the voltage is inside limits, ( ) is set to infinite.

STEP 3:
Find the controller with the minimum priority indicator such that ( ) ≤ ( ) ∀ ∈ . Then, for each controller ∈ , set ( ) = ( ) − ( ) . As shown, the priority indicator with the lowest values is subtracted from the other priority indicators. Consequently, controller with the minimum ( ) 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.

STEP 4:
In this step, controller with the minimum ( ) value, undergoes switching action according to Switching Logic 1. Similarly, if one or more additional controllers have a zero ( ) value, they will be also subject to a switching action according to Switching Logic 1.

STEP 5:
In this step, ( ) of controller j is updated based on its mechanical delay, as follows: Similarly, if one or more additional controllers have also undertaken a switching action in the previous step, their priority indicator value will be updated according to (16).

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.

V. 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 (refer to Part II). To accelerate the convergence, a novel sensitivity-based algorithm is proposed here.
To derive the three-phase sensitivity parameters of LVCs, let us firstly assume the network of Fig. 6. It includes an OLTC connected between the buses h-w, three single-phase 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 Fig. 4 at buses p, m, s. Let us define the vector of LVC control variables ( ) and voltages as: where: In Eq. (20), the sensitivity matrix , is of dimension 14x14 corresponding to the schematic presented in Figure 6. 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 Eq. (20) is presented in a supplementary document in [32].

VI. 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 sensitivitybased 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 in Part II that the sensitivity-based algorithm presents identical results with the generic algorithm of Section IV, while the iterations and computational effort are significantly reduced. Furthermore, we also compare our proposed algorithm with dynamic simulation to validate its accuracy in Part II of this research.
The algorithm consists of the following steps: INITIALIZATION: Similar to the initialization in Section IV. STEP 1: Similar to STEP 1 of Section IV. STEP 2: Similar to STEP 2 of Section IV. STEP 3: Similar to STEP 3 of Section IV. STEP 4: In this step, the change in voltage of the LVCs is predicted based on the operation of the DG utilizing the sensitivity relationship. As the DGs in the network react instantaneously, we firstly predict the variation of the positive sequence reactive power ( ) as well as the negative-and zero-sequence currents ( ) of DGs based on Eq. (21). In Eq. (21), and are the deviation of and from their set reference voltages, respectively.
Due to the coupling of LVCs, the variations and calculated in (21) cause voltage variations at the other LVC buses. By utilizing the sensitivity parameters and the results of Eq. (21), the predicted voltage variations ( , , ) of the LVCs are calculated, as shown in Eq. (22).
The new predicted voltages , , of LVCs are then calculated in Eq. (23): 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 ( ) , as presented in Switching Logic 2.

Switching Logic 2
If ∈ or ∈ then { It should be noted that the Switching Logic 2 is similar to Switching Logic 1 of Section IV. However, Switching Logic 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, Switching Logic 2 uses the predicted voltage value of controller ( ), which is calculated from Eq. (23), 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 Eq. (24). Note that Eq. (24) Finally, the DG control variables are updated as shown in Eq. (25), using the result of (24).
STEP 7: Similar to STEP 5 in Section IV. STEP 8: Execute one iteration (only) of the power flow with the updated LVC state from Switching Logic 2 and the updated , from Eq. (25). Go to STEP 2 until no further switching actions take place.
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 IV. 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 Part II.

VII. CONCLUSION
This Part I presents the theoretical development of 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 faster convergence and low computational complexity and is therefore expected to serve as an important tool in offline and real-time applications of distribution systems. Simulation case studies presented in Part II highlight the superior performance of the proposed approach against the existing power flow methods. Εq. (24) is directly derived from (A1), by making some simple algebraic manipulations.