Determination of the Thevenin equivalent in power grids using real records based on short circuit power

In recent years, methods for loadability determination based on real-time measurements have been proposed for long-term voltage stability scenarios. One of these methods corresponds to the reduction of the complete system into a basic Thevenin equivalent at a load bus. This paper proposes a method to determine the Thevenin equivalent at a load bus in power systems based on short circuit power ( S SC ) and its application to real records. The proposed Thevenin equivalent evaluation is performed by processing two consecutive sets of voltage magnitude and power measurements, both recorded at the same bus of a power system. This method establishes a Thevenin equivalent based on the determination of system’s equivalent voltage ( V S ) which is obtained from S SC . This condition properly accounts for variable power factor load conditions and constant power factor alike in contrast with existing methodologies which assume constant power factor loads. The performance of this method is evaluated with real records from the Ecuadorian electric power system and contrasted by simulation.


INTRODUCTION
Voltage instability is produced as a result of the inability of the transmission and generation systems to provide power to the loads. For a long-term voltage stability scenario, the loadability determination is an important topic continuously analysed for power grids [1]. In [2] the authors proposed a classification of the methods for voltage stability detection; these methods are classified by means of two measurement strategies: local and wide areas.
In both strategies, usually the fundamental criterion for loadability determination in real-time corresponds to the reduction of the complete system into a basic Thevenin equivalent (TE) at a load bus. Then, the power load is maximised when the load impedance (Z L ) matches the TE impedance (Z TH ).
Thus, Thevenin equivalent determination (TED) becomes an essential part to apply this impedance criterion. Several approaches have been proposed in this field. For instance, the authors in [3] proposed the application of least squares estimation to track the TE, while the authors in [4] applied the Tellegen theorem at two consecutive synchronised measurements. Algorithms to determine a TE valid for extra high voltage (EHV) systems are proposed in [5,6], where E TH is continuously updated based on the sign of the load variation. In [7] the authors obtain a TE using three local measurements to obtain an analytic solution for the Z TH . In [8,9] Kalman filter is used to improve the accuracy of least squares techniques for TE. A comparative study of some of these TE techniques is presented in [10]. The Thevenin impedance match criterion is still used as an extension for multiport systems [11], as applied for flowgates [12] or for the visualization of voltage stability margins as exposed in [13,14]. Similarly, some voltage stability indexes have been proposed [4,[15][16][17][18] based on the same basic Thevenin reduction concept.
It is important to indicate that all of these techniques consider load changes at constant power factor. Moreover, in [19] the authors demonstrated that the calculation of TED is possible only if there exist continuous load changes. Additionally, TE criterion must be adapted to give a coherent response under normal power system conditions which includes the effects of active power dispatch and changes on load power factor.
In [14,20,21], the authors indicated that the active power dispatch can produce an incorrect determination of TED because of the non-linear response in the sharing of power supplied between generators following load changes. Similarly, in [22] the authors exposed that changes on load power factor can produce incorrect TED results because the maximum loadability is not reached when the Z L and the determined Thevenin impedance (Z THd ) are equal. Furthermore, the authors explained in detail the advantage of determining loading from the non-linear power domain solution perspective with better accuracy than the existing methodologies.
On the other hand, a few practical implementations have been proposed to apply TED on real measurements. For instance, in [23,24] the authors presented a methodology where E TH is assumed to be almost constant. Then, by assigning different values of X TH , the E TH is evaluated until the constant condition is accomplished. This methodology is only applicable to EHV because the authors have neglected the effects of R TH . In [25] the authors have proposed a method for TED by means of PMU measurements based on the assumption that the power factor is nearly constant. The authors in [26,27] show a TED by means of using a modified TE impedance based on coupled single-port circuit which considers constant power factor load model.
In this context, this paper greatly improves the analysis of loadability proposed in [22]. In addition, we introduce a novel scheme to TED which starts from the determination of the system's equivalent voltage (V S ) obtained directly from S SC . To prove the suitability of this scheme, a practical implementation is applied to real-world measurements obtained in the electrical Ecuadorian power system.
It is important to indicate that our objective in testing the proposed method by means of simulation was not to find which bus would exhibit first a voltage instability problem; instead, we focused our interest on determining a possible range for TE and giving a reference to compare with results obtained from real records. This paper is organised as follows. Section 2 expands the analysis of loadability based on the S SC and gives the definition of V S . Both are equally important to construct the methodology of the TED. Section 3 details the proposed methodology to determine the V S and TE. The TED is evaluated by simulation and with real records from the Ecuadorian grid. These results are detailed in Section 4. Finally, the main conclusions of this work are summarised in Section 5.

SYSTEM'S EQUIVALENT VOLTAGE (V S ), LOADABILITY AND THEVENIN EQUIVALENT IMPEDANCE
This section greatly expands the analysis previously presented in [22] related to the methodology for loading estimation based on short circuit power (S SC ). The basic power system FIGURE 1 Two-bus transmission system depicted in Figure 1 was considered in [22] to unravel all the analysis.
Here, the basic voltage equation can be defined as in (1): Let m be the load-system coupling factor, which is defined as the cosine of the difference between the TE system angle (φ) and the load angle (θ) as shown in (2). Then it is possible to rewrite the basic voltage Equation (1) yielding (3).
The main idea, expressed in [22], corresponds to treating the basic power equations into a parametric space valid for any power system. The basic parameter corresponds to the short circuit power (S SC ) at the load bus defined as in (4): The basic voltage equation can be written as in (5): The variables V S and S S correspond to (6): The voltage solution for the quadratic Equation (5) is found to be (7). The variable V S is called system's equivalent voltage, and it is defined if the variables m and S S are known, i.e. if the S SC can be determined and because the load power S is known, then V S would be also defined and allowing to determine the TE: From (7), it is possible to define the power flow solution domain as shown in (8). Considering typical values for system FIGURE 2 Profiles for power system in S S space: (a) power flow domain and loadability limit and (b) load trajectories for exponential and ZIP load models from [22] angles 80 • ± 7 • [15] and reasonable load angles of 10 • ± 20 • [28], the factor m could take values between the range: −0.12 ≤ m ≤ 0.73, then the power flow domain could be visualised in detail in Figure 2(a): In the limit of (8), there is a unique power flow voltage solution. Consequently, this point represents the loadability limit, and it corresponds to the point of maximum power flow S Smax transferred between two points on the transmission grid in per unit of the S SC as indicated in (9): Under this domain, it is possible to draw the S S trajectories for load increments starting from an initial point until the loadability limit for any load model. However, there are some differences in these trajectories for load models with constant and variable power factor. Figure 2(b) shows the S S trajectories for ZIP and exponential load models analysed in [22] and based on the test model proposed in [29].
From Figure 2(b) for the ZIP load model (blue trajectory), corresponding to a constant power factor condition, the S S trajectory is linear until it reaches the loadability limit. This condition considers that active and reactive power varies at the same proportion to produce a fixed power factor.
Meanwhile, for the exponential load model (green trajectory), corresponding to a variable power factor model, the S S trajectory is non-linear. A small load change from an initial loading point depicted in Figure 2(b), leads to a change of m, S S and loadability limit. This means that the power grid behind the load has also changed.
The variable power factor condition is not considered in the existing methodologies for TED [22]. In consequence, for small load changes, we propose a ratio between Z L and the Thevenin impedance determined for those methods (Z THd ) as shown in (10), where m 0 corresponds the evaluation of m at the initial point: The basic criterion where the power load is maximised when Z L matches Z TH is always accomplished. For constant power factor (m 0 = m), the Z THd is determined correctly. Then, the Z THd must be adapted as shows (10) for variable load power factor.
Moreover, the loadability could be determined in the S S domain. For this, consider a loading fraction x as a per unit respect to the maximum power limit given by (9), i.e. x could have a value between: 0 ≤ x ≤ 1. This defines a locus for a per unit loading called S S0 as in (11) and depicted in Figure 3: From Figure 3, the loading locus is a function of x. If this value is close to the maximum loadability S Smax (e.g. x = 0.8), this locus is like a parallel line to the maximum loadability. On the other hand, if x is close to the origin (e.g. x = 0.3), the loading locus shows a non-linear behaviour due to m. This adaptive feature includes the effects of variable load power factor and determines the loadability and TE from the S SC .

Transformation and maximization of voltage-power stability indicator (VPSI S )
The system voltage-power stability indicator (VPSI s ), defined in (12) is established as a means to bring the system loading into the SS domain as described in [22]. This domain includes all the effects of system and load changes as given in (6): For the sake of simplicity, let A be A maximum VPSI S occurs when the derivative of VPSI S with respect to S S is equal to zero. Thus, by solving this condition with respect to n yields: The voltage exponent n can be used to generate a maximum VPSI S in a predefined loading condition by considering S S equal to S S0 in (13) as indicated in [22]. Furthermore, when considering real-measured voltage V and apparent power S values, (12) can be rewritten as (15): where the VPSI is Note that VPSI can be determined with real voltage magnitude and power measurements, both recorded at the same bus of the power system. Now, let us consider a logarithmic transformation of VPSI S (17): (17) that is being proposed here to improve the error criteria encountered in [22] for low loading values.
Equation (18) describes the trajectory of ln(VPSI S ) for small displacements around an operating point: where ΔS, ΔZ TH and ΔE TH are defined in (19). These equations consider a small displacement from a starting point assigned as 1 to a final point assigned as 2: Considering the partial derivatives of (18) with respect to S, Z TH and E TH , it is possible to obtain (20) from (18): This equation defines the behaviour of VPSI S for any changes in the power system, and a TE method based on VPSI S maximization as will be described next.

V S evaluation and TE determination based on VPSI S maximization
The voltage exponent n has a distinctive feature to generate a maximum VPSI S in a predefined loading S S0 as indicated in (11). Equation (20) can be solved by continuously increasing the loading fraction x, until (20) is approximately equal to zero. This ALGORITHM 1 To determine TE based on VPSI and S SC Consider φ 0 according to the power system nature, for example for EHV φ 0 = 90 • Step 1: Take magnitude values of V i , V i+1 , P i , P i+1 , Q i and Q i+1 from two measurements. Calculate apparent power S i , S i+1 and load angle θ i , θ i+1 Step 2: Starting from an initial loading fraction x = x 0 until x < 1, calculate the following terms: Step 2.1: m i , m i+1 , S S0 i , S S0 i+1 , A i , A i+1 , n i and n i+1 applying (2), (11), (13) and (14), respectively.
Step 3: Search the S S0 for minimum error and calculate S SC from (22).
Step 4: Determine V S as a function of S S from (7).
Step 5: Determine E TH as a function of V S from (23).
Step 6: Determine Z TH as a function of E TH and S SC from (24).
condition is used to determine the S SC and to make a TE based on VPSI S . Consider a system where the load is changing and the power system remains almost unchanged, so (20) can be rewritten as (21) dVPS .
By monitoring the trajectory of VPSI by means of (21), it is possible to determine the value of S when a VPSI maximum occurs and consequently a maximum VPSI S is reached. From this point the S SC is directly calculated using (22): Once S SC is known, it is possible to find directly V S using (7). Then, the TE can be constructed. At first E TH is obtained using (23) and finally Z TH is determined in (24): Based on the adaptive feature of maximising VPSI under system changes, the S SC and the TE are possible to obtain by applying the following process:

THEVENIN EQUIVALENT CALCULATION ON THE ECUADORIAN GRID USING REAL RECORDS
This section studies the application of the proposed method to the Ecuadorian power system. For this study the Ecuadorian grid is modelled on DIgSILENT PowerFactory based on A radial sub-grid was selected for the application of the methodology explained in Section 3. This corresponds to the Southeastern Zone, which includes the following substations: Gualaceo, Limón, Méndez and Macas, as shown in Figure 4(b). For this sub-grid, real measurements were obtained from the historical SCADA system by CELEC EP -Transelectric for a randomly selected day (January 22, 2019). The case was conditioned to adjust the initial power flow as close as possible to the real data provided on the SCADA system before steady-state simulation using continuation power flow. Then, the algorithm detailed in Section 3 was applied to the continuation power flow results to determine possible TE ranges. Next, this algorithm was applied to real records obtained from the SCADA system for the selected day. Finally, the results with real measurements were compared against the results obtained from simulations.

Steady-state simulation
Our objective in testing the proposed method was not to find which bus would exhibit a voltage instability problem first; instead, our interest was to determine a possible range for TE in the bus under study and find a reference point to compare with actual records. To demonstrate the validity of our approach, two buses were selected based on data availability and quality, and they correspond to: Gualaceo 23 kV Bus and Méndez 13.8 kV Bus. The initial power flow was adjusted as close as possible to the real SCADA measurements for the day under study and for two specific scenarios at 03:00 and 20:00 as detailed in Table 1.
Taking into account this adjusted grid, continuation power flows were performed considering that the load increases in two stages: Stage 1: Increase only the load terminal under study. Stage 2: Increase all loads of the system simultaneously.
There could be a lot of more possible stages, each one with different results. Because the distance to the maximum power transfer depends on how the load is increased, some authors recommended increasing the load in a predefined manner, representing the most probable stressing scenario based on historical and forecast data [30]. Real-world records always contain a certain degree of uncertainty and will affect the definition of a load increase scenario completely consistent with reality. In this sense, if the differences between simulations and real system are small, they can be used to modify and enhance the model, within the processes of system validation [31].
Accordingly, the simultaneous load increments are established to calculate a range for the Thevenin equivalent of the bus under study, which is considered sufficient for comparison purposes. Figures 5 and 6 show the simulation results for Gualaceo 22 kV Bus, for stages 1 and stage 2, respectively. The power and voltage bases correspond to 100 MVA and 22 kV, respectively. Figure 5(a) shows that the maximum loadability is reached when the load impedance and the Thevenin impedance match for the two scenarios. It is clear that in the presented methodology the Thevenin impedance and the loadability are correctly determined. Concurrently, Figure 5(b) displays detailed profiles for Z TH and E TH , whose calculated values present low differences between the two hours considered, these parameters are around 0.535 and 0.95 pu, respectively.
Finally, Figure 5(c), depicts the S SC and the loading fraction value x as a function of lambda. The S SC is around 1.68 pu and the loading increase in a non-linear profile from less than 0.1 until 1 pu, where the maximum power transfer is reached. Because most of the generation is away from the Gualaceo 22 kV Bus, the value of Z TH and S SC determined from the present methodology correspond to the values of Thevenin impedance and three phase short circuit power determined directly by means of any short circuit program, as it was detailed in [32].
As it was expected, for Stage 2 the results give different values. Figure 6(a) shows that Gualaceo 22 kV Bus does not reach its maximum loadability because other buses in the grid are closer to the maximum loadability. In consequence, the power flow has a solution until this global lambda factor is reached.  Méndez 13.8 kV Bus is inside the Southeastern Zone, and its results will be compared to the real results similarly as Gualaceo 22 kV Bus. Table 3 summarises the principal range results for Pascuales 69 kV Bus is one of the most loaded buses of the Ecuadorian grid and for this reason its loading is close to the Based on the simulation results, the loading depends on how the load changes. For the two load increment stages presented here, the results show very different solutions. In the first stage, the load is increased in an unusual scenario. However, it validates the solution found since the load power is maximised when Thevenin and the load impedances match. It should be noted that the simulation was performed considering the modelling of the entire Ecuadorian and Colombian system, power limits of the generators, voltage dependency of loads and, as indicated above, loads with a variable power factor model.
For the second stage, the criterion of maximising the load power based on impedances is also met, but in some of the buses of the system as a whole. However, this condition can be observed in some of the buses with greater system load.
Considering these results, the loading determined in each bus is a local value, but it involves the complete response of the system by means of the S SC and will depend on how the real load increases. In this sense, better approximate loading increase scenarios are possible as shown by contrasting the simulation results with real records.

Test with real records
Available records for voltage, active and reactive power profiles were extracted from the SCADA system for a randomly selected day (22 January 2019). These data correspond to the Gualaceo 22 kV and Méndez 13.8 kV Buses. The records were obtained from bay control units (BCUs) located on these buses as shown in Figure 4(b). The BCUs are connected to analogue current and voltage transformers with a measurement precision of 0.1%, and they reported valid values sampled at 1 minute; thus, there are 1440 data points per day. These records are processed as displayed in Figure 8. This figure shows two ways to report the results, one obtained by direct application of the methodology on the records, and the other based on the average determination by means of a moving window at 10 minutes. The 10-min processing window is established to determine a continuous profile for all calculated variables that will be used to compare with the simulation results.
The power and voltage records, for Gualaceo 22 kV and Méndez 13.8 kV Buses, are detailed in Figures 9(a) and 10(a),  Tables 2-4, the values indicated as real records correspond to the data processing for mobile average sliding window. Table 2 shows the values determined in Gualaceo 22 kV Bus for real records at 03:00 and 20:00. These results show that for this bus the loading increments proposed in the simulation Stage 2 are closer at 03:00 than 20:00 with respect to real load changes. Table 3 shows the values determined in Méndez 13.8 kV Bus for real records at 03:00 and 20:00. These results show that for this bus the loading increments proposed in the simulation Stage 2 are closer at 20:00 than 03:00 with respect to real load changes.
It is important to indicate that, in the definition of the loading changes scenarios, not only their direction (increments or decrements) should be considered but also their magnitudes which are a function of the power at each analysed load, as it was observed in the results for Gualaceo and Méndez Buses. Table 4 shows the total ranges for the variables for Gualaceo 22 kV and Méndez 13.8 kV Buses determined from real records. The maximum and minimum values from the day-long data records could be considered as a reference range for the TE in these buses for this complete studied day.

CONCLUSIONS
This paper proposes a novel method to determine the TE, based on V S and S SC . The S SC parameter space expands the basic matching impedance condition to include the effects of active power dispatch and changes on load power factor. This method has been applied to simulated as well as real-measured data. The loading determined by the proposed method corresponds to a local value evaluated in each bus, but it involves the complete response of the system by means of the S SC and will depend on how the load changes. In this sense, contrasting the simulation results with real records constitute a better approach to get more approximate loading increase simulation scenarios. Finally, it is important to mention that future work will include the global analysis of TED and the possible determination of fault conditions from power flow real records.