A Comprehensive Three-Phase Step Voltage Regulator Model with Efficient Implementation in the Z-Bus Power Flow

— This paper presents a comprehensive three-bus equivalent circuit model of three-phase step voltage regulators. The proposed model can be efficiently integrated in the Z-bus power flow method and can accurately simulate any configuration of step voltage regulators. In contrast to the conventional step voltage regulator models that include the tap variables inside the Y BUS matrix of the network, the proposed model simulates them in the form of current sources, outside the Y BUS matrix. As a result, the refactorization of the Y BUS matrix is avoided after every tap change reducing significantly the computational burden of the power flow. Furthermore, possible convergence issues caused by the low impedance of step voltage regulators are addressed by introducing fictitious impedances, without, however, affecting the accuracy of the model. The results of the proposed step voltage regulator model are compared against well-known commercial softwares such as Simulink and OpenDSS using the IEEE 4-Bus and an 8-Bus network. According to the simulations, the proposed model outputs almost identical results with Simulink and OpenDSS confirming its high accuracy. Furthermore, the proposed 3-bus equivalent model is compared against a recently published conventional step voltage regulator model in the IEEE 8500-Node test feeder. Simulation results indicate that the proposed step voltage regulator model produces as accurate results as the conventional one, while its computation time is significantly lower. More specifically, in the large IEEE 8500-node network consisting of four SVRs, the proposed model can reduce the computation time of power flow around one minute for every tap variation. Therefore, the proposed step voltage regulator model can constitute an efficient simulation tool in applications where subsequent tap variations are required.


Introduction
TEP voltage regulators (SVRs) are widely applied to low-(LV) and medium-voltage (MV) distribution networks to maintain voltages within permissible limits. They consist of autotransformers with adjustable turn ratios connected in several configurations e.g open-delta, closedelta, wye [1].

Literature Review on the steady-state SVR modeling
Only a few models have been proposed in the literature so far for the mathematical integration of SVRs into the power flow [1]- [4]. More specifically, the authors in [2] and [3] propose a mathematical formulation relating the primary and secondary voltages and currents, which is solved using the backward/forward sweep (BFS) method. In [4], a mathematical model of SVR is proposed, which is directly integrated in the Jacobian matrix, while the power flow is solved using the Newton-Raphson (NR) method.
However, these approaches are not applicable to the Z-Bus power flow since the SVR equations of [2]- [4] are not compatible with the formation of YBUS matrix. The Z-Bus power flow approach is a fixed-point iterative method, which is widely applied in distribution network applications, due to its high robustness, fast convergence, ease of implementation and low computation time [5]- [9]. Therefore, accurate and computationally efficient modelling of SVRs is significant in order for the Z-Bus method to maintain its superior characteristics.
Considering the integration of SVR models in the threephase implicit Z-Bus power flow, only a limited number of solutions has been proposed in the literature. A comprehensive SVR model is proposed in [1], which is applicable in all SVR configurations, while an open-delta SVR model is proposed in [10]. Although the aforementioned models are applicable in the Z-Bus power flow approach, the tap variables of SVR are simulated in the YBUS matrix of the network. As a result, a re-factorization (or inversion) of the YBUS matrix is required after every tap variation, thus increasing significantly the overall computational burden.

Challenges of efficiently modeling SVRs in the Z-Bus power flow
Ideally, the tap variables of the SVRs should be simulated outside the YBUS matrix of the network so that the refactorization of the YBUS matrix is performed offline, only once, accelerating the power flow calculation. The S compensation technique is usually implemented in the ZBUS power flow to avoid the re-factorization of the YBUS matrix every time that one or more elements change [11]. According to this technique, the variable elements are modeled as fictitious current sources outside the YBUS matrix, thus the YBUS matrix remains always constant.
Τhe compensation technique is applied in [7] and [11] to simulate transformers equipped with on-load tap changer (OLTC) in balanced networks (using single-phase power flow) and by the authors in [12] to simulate three-phase OLTC transformers in Yg-Yg configuration. However, in case of the SVRs, the implementation of the compensation technique results in the divergence of the power flow, due to the small impedance of SVRs, which in turn increases unacceptably the fictitious current sources. Similar conclusions about the divergence issues associated with the modeling of SVRs are derived in [13].

Technical contribution of this paper
This paper proposes a novel model, which overcomes the aforementioned divergence issues. The proposed model consists of a 3-bus equivalent circuit, which presents the following distinct characteristics: • It simulates the taps of SVR as fictitious current sources outside the YBUS matrix, thus avoiding its continuous refactorization after every tap variation. As a result, several real-time distribution management system (DMS) applications that require sequential tap variations are accelerated.  [20] [21], optimal feeder reconfiguration (OFR) and optimal power flow (OPF) [7], which are executed following the load variations and the variation of the volatile distributed generator (DG) power [22] [23]. More specifically, as soon as an estimation of the system load and the volatile DG power is obtained, DMS optimization functions compute new switching states of the local voltage controllers (e.g SVRs, OLTCs, capacitor banks, distributed generators) to optimize the system performance. The real-time power flow serves as an internal tool in the optimization functions since it is employed to simulate the effect of control actions on the state of the system and determines, if a new switching state of local voltage controllers (LVCs) is to be implemented [11]. Therefore, accurate and computationally efficient integration of LVCs in the power flow is essential for the efficient implementation of DMS optimization functions. The proposed SVR model can constitute an important tool in DMS applications, due to its high accuracy and its ability to maintain the structure of YBUS matrix constant saving significant computation time.
Another important application of DMS is the state estimation and monitoring of the distribution system. It is realized through the execution of power flow considering the different time delays of the LVCs of the network. More specifically, LVCs do not react instantaneously but they are assigned to react with different reaction delays [18] [19]. 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, 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. This process requires the execution of power flow many times by subsequently varying the taps of LVCs based on their reaction time delays, until the final state is reached [11] [18] [19] [35]. With the existing SVR models, a high amount of computation time would be required, due to the continuous refactorization of YBUS matrix every time that a tap variation occurs. This limitation is completely addressed with the proposed SVR model.

Paper structure
The rest of the paper is structured as follows: Section 2 explains the difference between the impedance value of an autotransformer and a two-winding transformer. Section 3 presents the proposed SVR model, while Section 4 explains how it is integrated in the Z-Bus power flow. Section 5 validates the proposed SVR model against Simulink and OpenDSS. Section 6 presents an analytical case study, in which the performance of the proposed SVR model is highlighted in the IEEE 8500-network. Finally, section 7 concludes the paper.

Impedance of autotransformers
Before we proceed with the analysis of the proposed SVR model, it is important to clarify the difference between a two separate winding transformer and an autotransformer. Fig. 1 depicts an ideal two-winding transformer and an ideal autotransformer. The impedances of the transformers are denoted as Z2wT and Zauto and they include the total leakage impedance of the primary and secondary winding referred to the secondary winding. The primary voltage (V1), the secondary voltage (V2) and the number of turns of the primary windings (N1) are identical for both transformers. Voltages and turn numbers for the two types of transformers can be determined according to (1) and (2). (1) It can be observed that the number of turns of secondary windings (N2wT and Nauto) differs for the two transformer types. Furthermore, assuming that the two transformers have an equal magnetic reluctance, the corresponding impedance ratio for the two transformer types is given by (3).   Fig. 2 depicts the impedance ratio of the two transformer types with respect to the output/input voltage ratio. The figure is obtained using (4) for an output/input voltage ratio between 0.9 and 1.1, which is a typical voltage range for SVRs. It can be observed that the autotransformer presents significantly lower impedance than the two-winding transformer. As a result of the low impedance of the autotransformer, the ZBUS power flow diverges when the compensation technique is applied to the SVR modeling [13].

Proposed SVR model
where the matrices , , are defined in Tables 1 and 2  The square matrix in (5a) includes the matrix , which is not constant due to the tap variations. Therefore, it is moved outside the square matrix and (6a) is obtained, where the current sources , are given in (6b).
Finally, the voltages and currents between the nodes m-s are expressed in (7), where the matrix is given in Table 2.
Eqs. (6) and (7) can be comprehensively expressed in the form of an equivalent circuit as shown in Fig. 5. The parameters of the equivalent circuit for both types of SVR are given in Tables 1, 2 As shown in Table 3, the current sources depend on the matrices , . To ensure the convergence of the power flow, the values of these matrices should be sufficiently small so that the current sources do not take an extremely large value. That is achieved by selecting a sufficiently large value of . Essentially, the role of is to enhance the convergence by artificially increasing the value of .     Table 3 Current sources for the two types of SVR

Implementation of the proposed SVR model in the Zbus power flow: A step-by-step example of a 4-Bus network
In this section, the implementation of the proposed model is practically explained in the 4-Bus network of Fig. 6a. A closed delta SVR of type A is assumed to be connected between the nodes 2-3.
Step 1: Modification of the network to include the proposed SVR model. The conventional SVR configuration of Fig. 6a is modified, as shown in Fig. 6b, to represent the proposed three-bus equivalent SVR model. As shown in the figure, the proposed SVR model has an additional fictitious bus 5. Note that no load is connected to this fictitious bus. Buses 2, 5 and 3 correspond to buses p, m and s of the equivalent circuit of Step 2: Formulation of the YBUS matrix of the network. The network of Fig. 6b is described by (8), as follows: The square matrix of (8) is the YBUS matrix of the network and it includes the 3x3 admittance matrices of the lines between the buses 1-2 ( ) and 3-4 ( ). Moreover, it includes the 3x3 matrices and calculated according to Table 2, where and correspond to and , respectively. It is noted that in (8) Step 3: Formulation of the current vectors. The current vectors of (8) are given in (9). Step 4: Converting the YBUS matrix to a solvable form. Assuming that the bus 1 is the slack bus, the first three rows of (8) are removed and (10) is obtained, as follows: Subsequently, the product [ ] • is removed from both equation sides of (10) and (11) is derived. Please note that the column [ ] corresponds to the first column of the YBUS matrix of (10).
Step 5: Inversion of the YBUS matrix. By inverting the square YBUS matrix at the right-hand side of (11), equation (12) is obtained.
It is noted that the square YBUS matrix consists of constant elements since the tap variables of SVR are simulated as compensating currents in the current sources , , .
Therefore, YBUS is inverted only once and not after each tap change. Selecting a sufficiently large value of , the current sources , , are adequately small to ensure the convergence of the power flow.
Step 6: Iterative power flow solution. Eq. (12) is iteratively solved until a certain preset tolerance is reached. In (12), k denotes the iteration number. In case of a new tap variation, only the current sources , , are updated.

IEEE 4-Bus network with 1 SVR
The results of the proposed model are compared against Simulink (version 2018a) and OpenDSS (version 9.1.0.1), using the IEEE 4-Bus network. Simulink and OpenDSS were selected for the validation since they are two of the most widely-used and accurate commercial simulation software tools. The former is a powerful tool performing time-domain analysis, while the latter is a well-established software for the power flow analysis of grids. In Simulink, the SVR was formed by connecting two or three ideal transformer blocks depending on the SVR configuration e.g closed-delta, opendelta, wye. The transformers have pre-specified winding ratio to correspond with those of SVR. The impedance of SVR was modeled as an external impedance connected in series with the secondary winding, as shown in Fig. 6a. A similar rationale was followed when modeling SVRs in OpenDSS. Some clarifications about the applied parameters are provided below: • The unbalanced loads of [17] are used in this paper. In Open-Delta and Closed-Delta configuration the loads are phase-to-phase connected, while in wye configuration they are phase-to-neutral.
• The SVR is placed in the position of the transformer, between the buses 2-3 [17]. The impedance of the SVR ( ) is calculated based on the data of the single-phase two-winding transformer given in [17], which has the following parameters:  Table 4 for the SVR of type A, and in Table  5 for the SVR of type B. Indicatively, only the phase-to-phase (or phase-to-neutral) voltages of the 4 th bus are depicted here. The remaining buses present similar accuracy. As shown, the results of the proposed approach are in full agreement with those of Simulink and OpenDSS confirming the accuracy of the proposed model.

8-Bus network with 3 SVRs
The proposed SVR model is further validated in the 8-Bus network of Fig. 7, consisting of 3 SVRs and 4 unbalanced constant impedance loads. It is a 4-wire network with the following assumptions: 1) every bus is considered perfectly grounded, 2) the mutual impedance between the neutral and phase conductors is zero. The parameters of the network are quoted in Table 6, while the loads are given in Table 7. The taps of SVRs for the simulated configurations are presented in Table 8.  The proposed approach is compared against Simulink and OpenDSS and the results are presented in Table 9. Indicatively, only the voltage magnitude of the three phases of bus 8 are depicted, nevertheless, the remaining buses present similar accuracy. In Table 9, the results of the proposed approach are presented in regular font, the results of Simulink in bold, while the results of OpenDSS in green color. As shown, although the network consists of 3 SVRs, the proposed model produces almost identical steady-state power flow results with Simulink and OpenDSS confirming its accuracy. Finally, the currents of the lines of the 8-bus network are shown in Table 10 for the proposed approach (regular font), Simulink (bold) and OpenDSS (green). The SVRs are in closed-delta connection of type A and their tap settings are provided in Table 8. It is confirmed that the proposed approach is very accurate since all results are identical with those of commercial software tools.

Convergence properties of the proposed SVR model
The convergence characteristics of the proposed model are presented here, with respect to the values of the additional impedance Zadd introduced in Section 3. Fig 8 depicts indicatively the voltage of phase A at the last bus of the 8-bus network, throughout the iterative process of power flow, for several values of . More specifically, Fig. 8a depicts the voltage for low values of (from 0.1j to 0.3j), while Fig.  8b for high values (from 0.5j to 10000j). Although Zadd takes only inductive values in Fig. 8, similar results are derived with resistive values of similar magnitude. It is clarified that a flat start was considered for all the nodes of the network. Based on the figure, the following observations are derived: • For low values (see Fig. 8a) the Z-Bus power flow diverges. It occurs because for low values, the matrices , are increased (see Table 2). As a result, the current sources of Fig. 5 (see Table 3) are increased as well, resulting in the divergence of the power flow.
• The algorithm presents fast convergence for all ≥ 0.5 , as shown in Fig. 8b. Thus, can be arbitrarily selected over a very wide range of values, making the selection of a proper value a very simple process. • The algorithm converges to the same solution for all ≥ 0.5 . This is logical since an impedance − is connected in series with (see section 3), canceling completely the influence of . In fact, is added to reduce the values of the current sources of Table 3, and therefore, enhance the convergence of the power flow, without affecting the power flow results.

Description of the network
The performance of the proposed SVR model is tested in the IEEE 8500-Node network, which is a real large-scale distribution network. It includes originally 4 SVRs and 4 capacitor banks. SVR 1 is assumed to be connected in wye configuration, SVR 2 in open delta configuration, while SVRs 3 and 4 are connected in closed delta configuration. SVR 1, SVR 2, SVR 3 and SVR 4 are connected between the buses 1-2, 405-406, 1057-1058 and 1311-1312, respectively. The capacitors have been connected in buses 15, 146, 644 and 1298. The topology of IEEE 8500-network is shown in Fig.  9. In order to investigate the interactions of SVRs with the DGs, we connected 4 DGs to the buses 100, 350, 835 and 1600.
To get a sense about the position of each bus inside the network, their distance from the substation is presented in Fig. 10. More details about the data of the network e.g lines, loads are provided in [17]. All buses of the network are considered perfectly grounded.

Modeling of the local voltage controllers
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 [24]. 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 . In this study, all capacitor banks are assumed to operate in voltage control mode, and all have the same capacitances. The capacitance of each phase is given in Table 12.
The steady state modeling of DGs is usually categorized depending on the generated voltage/current profile as well as the power profile, as follows:  [26].
In this study, the DGs are assumed to generate different current and voltage profiles as shown in Table 11. The active power and the droop gains of DGs 2 and 3 are provided in Table 12.

Switching delays of local voltage controllers
LVCs are usually configured to react with pre-specified time delays. This is done in order to coordinate the switching actions of LVCs and to avoid unnecessary switchings, due to temporary voltage variations.
The switching capacitors are set with an intentional time delay such that the capacitor switching is activated, after a pre-specified time delay, from the instant that a voltage violation is sensed [18] [19] [27]. The intentional time delay, as well as the reference voltage and bandwidth of capacitors are quoted in Table 12.
Similarly, SVRs 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 step-up/down switching action. If the voltage remains outside the bandwidth, additional step-up/down actions are executed with a mechanical delay [18] [19]. The intentional and mechanical time delay, as well as the reference voltage and bandwidth of SVRs are quoted in Table 12.
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 [18] [19] [27]- [30].
It is clarified that the intentional time delays of capacitors and SVRs were selected based on their distance from substation. Near to substation LVCs have shorter intentional time delays [31, page 29], while delays are increased for longer distances.
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 [32]- [34]. In this study, we assume that each phase is independently controlled.

Power flow algorithms for distribution networks with local controllers
The inclusion of the LVCs into the power flow is known to create multiple power flow solutions, as a result of the different reaction speed of LVCs and the discrete switching states of SVRs, OLTCs and capacitor banks. The state of-theart for simulating local voltage controllers computes the most likely operational solution by considering the different precedence of the controlled devices, according to their reaction time delays [11].
The time delays of LVCs were introduced for the first time into the power flow studies in 2000 by a research group of Siemens [35] and it is considered the state-of-the-art so far in the steady state modelling of distribution networks with LVCs. More specifically, LVCs were divided into different groups based on their switching time delays. The LVCs with the fastest reaction belong to the first group, while the slowest ones belong to the last group. The method of [35] was improved later in [7] and [11] by introducing sensitivity parameters to accelerate the convergence of the power flow. In the aforementioned approaches, a single-phase power flow is applied neglecting the network unbalances.
Authors in [18] [19] developed a three-phase power flow algorithm, which precisely considers the actual switching sequence of LVCs based on their switching time delays. Moreover, it considers the various operational modes of DGs as well as the unbalances that usually exist in distribution networks. In this case study, we applied the method of [18] [19] for calculating the power flow of the examined 8500node network considering the LVCs and DG characteristics of Tables 11 and 12.

Power flow results
Initially, the power flow is executed considering all the LVCs at their initial switching state. In this example, the taps of all SVRs are initially assumed at 0 position, while the capacitors are OFF. After the convergence of the power flow (at 43 th iteration), the voltages of LVCs are checked and a switching action is undertaken from the LVCs with the shortest time delay, whose voltage bandwidth was violated. Then, the power flow is executed again, the voltages of LVCs are re-checked and a switching action is activated in order of priority, until all voltages lie inside their bandwidth and no other switching actions are required.
All switching actions of the examined network are summarized in Table 13 and explained below. Furthermore, the phase voltages of SVR 1, SVR 3 and capacitor 4 as well as the reactive power of DGs during the power flow process are indicatively depicted in . The voltages of the other LVCs are not depicted due to space limitation. The algorithm converges in 309 iterations. The proposed model as well as the three-phase model of [1] are applied for the modelling of SVRs and the results of both models are discussed in the next sub-section. Switching action 1: The first switching action is taken over by the SVR 1 since it violates its voltage bandwidth and has the shortest time delay. More specifically, as shown in Fig.  11, after the power flow has converged at 43 th iteration, the voltage of each phase of SVR 1 is near 7200 V, which is below the bandwidth of 7465V-7535V. Therefore, SVR 1 step up the voltage of all phases by one tap in order to increase its voltage toward its bandwidth. Switching actions 2-6: After the first switching action of SVR 1, its voltage remains below the bandwidth. Therefore, SVR1 steps up the voltage of all phases sequentially, with a mechanical time delay, until the voltages of all phases lie inside bandwidth, as shown in Fig. 11. Switching actions 7: After the switching actions of SVR 1, the voltage of SVR 2 lies inside the bandwidth and no switching is performed by SVR 2. On the opposite, the phase voltages of SVR 3 remain outside the bandwidth resulting in a step-up function by SVR 3, after the time delay of 60 sec (the time delay of SVR 3). 14) is reduced forcing the voltage of phase C of capacitor 2 below its bandwidth. Thus, phase C of capacitor 2 is switched ON 45s later at 180s.

Discussion about the power flow results
The voltage profile of the network after all LVC switching actions have been executed is depicted in Fig. 15, using the proposed model as well as the model of ref. [1] for the simulation of SVRs. The two models produce exactly the same results confirming the accuracy of the proposed model. However, a huge difference is observed in the computation time of the two models. As shown in Table 13, 19 switching actions are executed in the examined case study by the SVRs until their voltage is stabilized inside their bandwidths. The model of [1] includes the taps of SVRs in the YBUS matrix, and therefore, the re-factorization (or inversion) of the YBUS matrix is executed 19 times during the power flow process. On the contrary, the proposed model simulates the taps as fictitious current sources and the YBUS matrix remains constant throughout the whole power flow process. As a result, the power flow is executed in only 154.5s using the proposed SVR model (0.5s per iteration for 309 iterations) and in 1408.5s using the model of [1]. It is clarified that the inversion of the YBUS matrix of the examined network is performed in around 66s in Matlab, due to the large size of the network. A summary of the computation time of the power flow for the two SVR models is provided in Table 14. All simulations were conducted in a PC with an Intel Core i7, 3.4GHz CPU and 16GB DDR3 RAM.
Finally, another important thing to note in Figs 11-14 is the fast convergence of the power flow when the proposed SVR model is applied. More specifically, the power flow converges very fast to the new solution after each SVR tap variation. Although the proposed model uses fictitious current sources to simulate the taps of SVRs, the convergence of power flow is not sacrificed.

Conclusion
This paper presents a novel 3-Bus equivalent circuit for modeling three-phase SVRs in the ZBUS power flow. The model is applicable in all SVR configurations. Contrarily to the conventional SVR models that include the tap variables inside the YBUS matrix, the proposed model simulates them in the form of current sources, outside the YBUS matrix, avoiding its re-factorization or inversion after each tap change. Furthermore, possible convergence issues caused by the low impedance of step voltage regulators are addressed by introducing fictitious impedances, without, however, affecting the accuracy of the model.
The proposed model has been validated against wellestablished commercial softwares such as Simulink and OpenDSS confirming its high accuracy. Furthermore, it has been compared against a recently published conventional SVR model. Simulations results in the IEEE 8500-node network indicate that a time reduction around 1 minute for each tap variation is achieved from the proposed SVR model compared with the conventional one.
To summarize, the proposed SVR model combines three important characteristics: a) increased accuracy since it considers the precise three-phase SVR circuit without any simplification, b) generic applicability since it can be applied in all SVR types and configurations, c) reduced computation time. Therefore, it can be applied to accelerate the execution of several real-time DMS applications such as state estimation, OPF, VVC, OFR, voltage stability, heuristic optimization, which require the subsequent tap variation of SVRs.