Analyzing Impact of Distributed PV Generation on Integrated Transmission & Distribution System Voltage Stability—A Graph Trace Analysis Based Approach

The use of a Graph Trace Analysis (GTA)-based power flow for analyzing the voltage stability of integrated Transmission and Distribution (T&D) networks is discussed in the context of distributed Photovoltaic (PV) generation. The voltage stability of lines and the load carrying capability of buses is analyzed at various PV penetration levels. It is shown that as the PV generation levels increase, an increase in the steady state voltage stability of the system is observed. Moreover, within certain regions of stability margin changes, changes in voltage stability margins of transmission lines are shown to be linearly related to changes in the loading of the lines. Two case studies are presented, where one case study involves a model with eight voltage levels and 784,000 nodes. In one case study, a voltage-stability heat map is used to demonstrate the identification of weak lines and buses.


Introduction
Traditionally, power transmission, and distribution networks have been modeled and solved separately. Legacy transmission system power flows and tools model distribution systems as a single-phase bulk load on a transmission bus. Likewise, in distribution system power flow analysis, the transmission system is assumed as a slack bus. In the past, such modeling was justified using assumptions that the transmission system operates balanced, all generation lies in the transmission network, and the power flow is unidirectional from transmission to distribution. Most of the voltage stability work that matured in the 1980s and 1990s was developed for transmission networks under these assumptions [1][2][3][4][5]. However, the landscape of today's grid is rapidly changing with the advent of Distributed Energy Resources (DERs), and the assumptions made previously are not justified for smart grids.
Recent studies have demonstrated phase imbalance in transmission systems for infrastructures with increased penetration of DERs, and bi-directional power flows between transmission and distribution [1]. Specifically, there is a significant push back [1,[6][7][8][9] on assuming a balanced operation for transmission systems and subsequently modeling them in the positive-sequence domain. In particular: • Single phase transmission system models are not suitable to evaluate the impacts of distribution connected unbalanced loads and single-phase DERs on transmission systems [10][11][12][13].
• Integrated T&D models can have multiple voltage levels, diverse topologies (radial, lightly meshed, heavily meshed) and varying impedance characteristics. It is well known that Jacobian-based power flow approaches converge well for certain topologies and impedance characteristics only. For example, the Newton-Raphson method exhibits poor convergence characteristics for radial/lightly meshed systems with high R/X ratios [26,27]. • Integrated T&D systems operate unbalanced, and thus, the resulting matrices do not possess the symmetry observed for balanced systems. This renders the matrix decomposition difficult, or impractical, as the size of the system grows.

•
The large size and unbalanced nature of Integrated T&D systems implies that any matrix-based approach must be able to manage millions of rows and columns of a multi-phase system, and without significant advances in matrix-based computations, rendering the analysis complicated and slow.
Previous work discussed utilization of Graph Trace Analysis (GTA) power flow for integrated T&D systems, where the entire system is modeled as unbalanced, multi-phase, with varying topologies [1,17,28]. A matrix-free, GTA-based power flow was employed and found to be computationally efficient, robust, and capable of solving multiple topologies in a single model. The work here builds on previous studies and makes the following contributions: • Voltage stability of integrated T&D systems is discussed in the context of variations in distributed PV generation in a system. Specifically, how the PV variability relates to changes in static voltage stability margins, leading to empirical relations and conclusions. • A method to identify weak regions in a system is demonstrated. The method involves the use of selected load buses in different segments of the system, and employs a voltage stability heat map that provides a visual picture of vulnerable lines and buses in the system. This paper is organized as follows. Section 2 presents a brief review of the voltage stability margin concept, formulation, and related terminologies. Section 3 presents two case studies and analyzes the results. Finally, conclusions are presented.
All the concepts are demonstrated using sanitized models of actual utility systems. All of the computations are implemented in the academic version of the Distributed Engineering Workstation (DEW) [29] powered by the GTA power flow solver.

Static Stability Margin
The formulation of the stability margin follows from [30], where criteria for convergence of power flow for radial distribution circuits is discussed. Extending the concept to all power carrying circuit elements with finite impedance in an integrated T&D model, stability margin can be formulated by considering the element subsection shown in Figure 1. For simplicity, only a single-phase section is considered here, but stability margin can be calculated for each phase of a multi-phase line. This is important since modern transmission systems can exhibit phase imbalance due to the presence of both unbalanced loads and unbalanced Distributed Energy Resources (DERs).
Let us suppose that the kth circuit element shown in Figure 1 is connected between buses A and B. V A is the sending end voltage and V B refers to the load end voltage. The current flowing in the line is denoted by I k . The admittance of the element is Y k = G k + jB k . At the end of this section four types of loads are realized:

3.
Constant power equivalent load S Tk , which represents the sum of all the loads connected downstream of section k.

4.
Constant power equivalent load S Sk , which represents the sum of all the losses downstream of section k. This paper is organized as follows. Section 2 presents a brief review of the voltage stability margin concept, formulation, and related terminologies. Section 3 presents two case studies and analyzes the results. Finally, conclusions are presented.
All the concepts are demonstrated using sanitized models of actual utility systems. All of the computations are implemented in the academic version of the Distributed Engineering Workstation (DEW) [29] powered by the GTA power flow solver.

Static Stability Margin
The formulation of the stability margin follows from [30], where criteria for convergence of power flow for radial distribution circuits is discussed. Extending the concept to all power carrying circuit elements with finite impedance in an integrated T&D model, stability margin can be formulated by considering the element subsection shown in Figure 1. For simplicity, only a singlephase section is considered here, but stability margin can be calculated for each phase of a multiphase line. This is important since modern transmission systems can exhibit phase imbalance due to the presence of both unbalanced loads and unbalanced Distributed Energy Resources (DERs).
Let us suppose that the kth circuit element shown in Figure 1 is connected between buses A and B.
is the sending end voltage and refers to the load end voltage. The current flowing in the line is denoted by . The admittance of the element is = + . At the end of this section four types of loads are realized: 1.
-Constant power load.

2.
= + -Admittance representing constant Impedance load.  For section k, the power flow balance equation is given by: The power flowing into the receiving end and current through the element are given as: * = + | | 2 * (2) From (2) and (3), it follows that: Hence, and are given as: For section k, the power flow balance equation is given by: The power flowing into the receiving end and current through the element are given as: From (2) and (3), it follows that: Hence, P k and Q k are given as: where: Now solving (5) and (6) results in: where: Now let us define: where ∅ is given by: Thus: Only the positive sign is considered for the radical of (9) since that provides the desired physical solution. From (7)-(10), it follows that ∅ must be real and positive. For ∅ to be real, the following must hold: Since X and Z are positive: For ∅ to be positive, the following must hold: Equations (11) and (13) must be satisfied to have a valid power flow solution. It is found that when a system reaches a loading condition for which a steady state power flow solution does not exist, then (11) is violated for one or more circuit elements in the system. In this paper, the quantity Y 2 − 4XZ will be defined as the stability margin parameter, and given as: The parameter γ needs to be expressed in per-unit so that its value can have a sensible interpretation. In this regard, the value that (14) takes when V A = 1 p.u. and S k = 0 is defined as the base value γ 0 , for the per-unit expression. Note that in doing so γ 0 is calculated as: In our discussion γ is always considered as expressed in per-unit. As the loading on the system is increased, the stability margin of the circuit elements starts decreasing. When the stability margin approaches 0 for any circuit element, the steady state power flow solution ceases to exist.

Case Studies and Results
The first case study analyzes system stability in the context of PV generation variability using stability margin concepts. The second case study presents a method to evaluate the voltage stability of a system using voltage stability heat maps, by which buses and lines that are in trouble may be easily identified. For both case studies, the load is modeled as a constant current load where the power varies directly with voltage magnitude.

Case Study 1
Consider the T&D model shown in Figure 2. This system consists of 784,000 nodes and eight voltage levels, from transmission to distribution. For the analysis, four subcases are studied: • Stability margins at base load without PV.

•
Stability margins at five times increased loading (compared to base load) without PV.

•
Stability margins at five times increased loading (compared to base load) with PV. • Stability margins if loading is increased beyond five times.

Case Study 1
Consider the T&D model shown in Figure 2. This system consists of 784,000 nodes and eight voltage levels, from transmission to distribution. For the analysis, four subcases are studied:


Stability margins at base load without PV.  Stability margins at five times increased loading (compared to base load) without PV.  Stability margins at five times increased loading (compared to base load) with PV.  Stability margins if loading is increased beyond five times.    Adding 67,000 PVs to the distribution transformers throughout the system of Figure 2, of which approximately 6000 exist today, brings the system up to a level of energy independence, and with five times increased loading, the power flow analysis is performed again. Some PVs in this model have volt/var controllers, while most are modeled at unity power factor. The settings of all controllers during the study were kept the same as provided by the utility to keep consistency. The stability margins undergo a slight increase (due to increased voltage as PVs are added). The increase for the five lowest stability margins (which all lie in the area highlighted in Figure 2) are tabulated in Table 2. Note that the lines in Table 2 are the same as the ones in Table 1. For a selected bus where the transmission line with the least stability margin terminates, voltage stability nose curves with and without PV are plotted in Figure 3. Three different PV penetration levels, 6600, 21,200, and 67,000 PVs, respectively, are considered in Figure 3. That is, Figure 3 illustrates the impact of increased PV generation on voltages at the selected system bus. As distributed PVs are added, the voltages at the tip of the nose curve increase, resulting in an increase in the static stability margin and the shifting of the curve to the right and up. For a selected bus where the transmission line with the least stability margin terminates, voltage stability nose curves with and without PV are plotted in Figure 3. Three different PV penetration levels, 6600, 21,200, and 67,000 PVs, respectively, are considered in Figure 3. That is, Figure 3 illustrates the impact of increased PV generation on voltages at the selected system bus. As distributed PVs are added, the voltages at the tip of the nose curve increase, resulting in an increase in the static stability margin and the shifting of the curve to the right and up.    Figure 4 plots the maximum loading of the selected transmission line against the stability margin for different PV penetration levels. It shows that the relationship between an increase in penetration level of distributed PV and an increase in stability margin is non-linear. From figure 4 it can be seen that the change in stability margin of a transmission line is approximately linearly related to a change in loading of the transmission line within the region that ranges from 70-100% stability margin. Table 3. Linear regression equations and goodness of fit for all scenarios of Figure 4.

Scenario Linear Regression Equation for X ∈ ( , )
Goodness of Fit  Table 3 presents an evaluation of the linear regression used in the region of 70-100% stability margin, where X = Stability margin (%) and Y = Transmission line loading (KW). The goodness of fit as evaluated by using the standard R-squared ( 2 ) metric is found to be greater than 95% for all scenarios of Figure 4. Within this region, a sensitivity factor denoted by can be defined as: For different PV penetration levels in Figure 4, this sensitivity factor is computed to be:  evaluated by using the standard R-squared (R 2 ) metric is found to be greater than 95% for all scenarios of Figure 4. Within this region, a sensitivity factor denoted by σ can be defined as: Change in line loading Change in line stability margin Table 3. Linear regression equations and goodness of fit for all scenarios of Figure 4. For different PV penetration levels in Figure 4, this sensitivity factor is computed to be: This sensitivity, or slope, defines the line flow changes in kW per 1% change in stability margin. The above numbers show that these sensitivities at various PV penetration levels are similar. For every 1% drop in stability margin, the line loading increases by 289-355 kW for various PV penetration levels. It is also noted that the trend below 70% stability margin loses this linearity. The slope in this region (less than 70% stability margin) is almost the same for all PV penetration scenarios. This is found to be 595 kW (1.8% o f line capacity).

Stability Margins if Loading is Increased beyond Five Times
As loading is increased beyond five times the base load the power flow stops converging. Summarizing the different scenarios discussed for the system under study, it can be noted that:

•
As loading is increased on the system, stability margins decrease due to decrease in voltages.

•
As distributed PVs are added to the system, stability margins increase slightly due to increased voltages as a result of PV generation.
The change in the stability margin of a transmission line is approximately linearly related to the change in the loading of the transmission line for variations in the stability margin from 70-100%. The slope of the linear relation provides the change in line loading per percent change in stability margin.

Case Study 2
This case study uses the model shown in Figure 5 to demonstrate using stability margins in establishing a heat map that highlights vulnerable areas in the system, and corresponding weak lines and buses. This model consists of 1600 lines, 15 generators, and three voltage levels. All the PVs are modeled as a lumped negative load at the unity power factor.

Case Study 2
This case study uses the model shown in Figure 5 to demonstrate using stability margins in establishing a heat map that highlights vulnerable areas in the system, and corresponding weak lines and buses. This model consists of 1600 lines, 15 generators, and three voltage levels. All the PVs are modeled as a lumped negative load at the unity power factor.

Stability Analysis with Localized Loading Increases
To establish peak load carrying capability at different points of this system, load buses are selected at eight different locations, as shown in Figure 5. The locations of these buses are selected such that:

•
The area of the entire system is covered.

•
Potential hotspots, such as nodes with multiple transmission lines terminating, are analyzed.

•
The transmission bus has a significant level of PV generation in the distribution system.
The load at each of the selected buses is increased, one bus at a time, until power flow diverges, and the steady state voltage stability plot for the bus is constructed. In particular, the voltage, load, and stability at the tip of the nose curve are recorded for each of the eight buses. The power factor used for increasing the load at each bus is 0.9.
The steady state voltage stability nose curves from the localized loading analysis are plotted in Figure 6 for load locations 3 and 8. As the load is increased the GTA power flow is able to solve up to and beyond the tip of nose curve (the maximum load carrying point). The values recorded for the tips of the nose curves will be used as reference conditions in constructing a voltage stability heat map, to be considered in what follows. and stability at the tip of the nose curve are recorded for each of the eight buses. The power factor used for increasing the load at each bus is 0.9.
The steady state voltage stability nose curves from the localized loading analysis are plotted in Figure 6 for load locations 3 and 8. As the load is increased the GTA power flow is able to solve up to and beyond the tip of nose curve (the maximum load carrying point). The values recorded for the tips of the nose curves will be used as reference conditions in constructing a voltage stability heat map, to be considered in what follows.

Stability Analysis with Uniform Loading Increases across Entire System
Next, the load of the whole system is uniformly increased till the power flow diverges. As the load is increased, it is observed that beyond the scale factor of 1.5 the power flow diverges. The 14 lines with the lowest stability margins are noted, and these lines will be used in constructing a voltage stability heat map in the next section.

Heat Map Generation
The 14 lines with the lowest stability margins, determined in the previous section, are used here to generate what is referred to as a steady-state, voltage stability, heat map. The heat map is illustrated in Figure 7. In the heat map, the load buses, represented in the columns, are sorted in order from buses that have lines with the highest stability margins to buses that have lines with the lowest stability margins. In sorting the buses, the minimum stability margins of lines connected to a bus are noted and used for the sorting, where outliers are disregarded. For an example of an outlier, consider load bus-2 where the stability margin used for sorting is chosen to be 70.615 instead of 12.498, where the latter is treated as an outlier.
In Figure 7 the lines, represented by the rows, are sorted in order from highest to the lowest stability margins across the buses. Note that these stability margins are recorded at the tip of nose curves for the localized loading cases. The heat map shows that the most critical portion of the system

Stability Analysis with Uniform Loading Increases across Entire System
Next, the load of the whole system is uniformly increased till the power flow diverges. As the load is increased, it is observed that beyond the scale factor of 1.5 the power flow diverges. The 14 lines with the lowest stability margins are noted, and these lines will be used in constructing a voltage stability heat map in the next section.

Heat Map Generation
The 14 lines with the lowest stability margins, determined in the previous section, are used here to generate what is referred to as a steady-state, voltage stability, heat map. The heat map is illustrated in Figure 7. In the heat map, the load buses, represented in the columns, are sorted in order from buses that have lines with the highest stability margins to buses that have lines with the lowest stability margins. In sorting the buses, the minimum stability margins of lines connected to a bus are noted and used for the sorting, where outliers are disregarded. For an example of an outlier, consider load bus-2 where the stability margin used for sorting is chosen to be 70.615 instead of 12.498, where the latter is treated as an outlier.
In Figure 7 the lines, represented by the rows, are sorted in order from highest to the lowest stability margins across the buses. Note that these stability margins are recorded at the tip of nose curves for the localized loading cases. The heat map shows that the most critical portion of the system lies in the vicinity of load buses 6 and 7. This corresponds to the lower right corner of the model in Figure 5. The highest stability margins correspond to load bus-2 and load bus-8.
In Figure 7, the Max Net load row (fourth row from the top) of the heat map shows the maximum loading capability of each load bus (i.e., loading level after which the power flow starts diverging), where load bus-4 is chosen as the reference since it is capable of supporting the largest load among all the buses considered. Thus, bus loading is represented as a percentage of the maximum net load for load bus-4. For example, if load bus-4 supports a maximum load of 100 MW, then load bus-5, rated at 60%, would support a maximum load of 60 MW.
The highlighted cells with bold font in the heat map indicate that the line lies in the path connecting the two highlighted load buses (the impedance path concept), and none of the other load buses in the matrix are part of that particular path. For example, line-3 lies between load bus-5 and load bus-8, line-8 lies between load bus-2 and load bus-4, and so on.
The stability margins at the existing loading, or operating point, of the system are listed in the second column of the heat map, whereas the loads for the operating point are listed in the second row, labeled Net Load. By comparing these values against the maximum loads for each bus, the additional bus load that can lead to system instability can be computed. For example, the stability margin of line-1 is 68.4 MW and that of line-2 is 68.9 MW at the current operating condition. The loads at this current operating point at all load buses are given in the second row, i.e., 54 MW at bus-1, 14 MW at bus-2, and so on. If the maximum net load at load bus-2 reaches 138 MW, then voltage instability occurs.
Another insight provided by the heat map is the loss of renewable generation that can lead to instability. For example, bus-2 has 150 MW of renewable generation (shown in row 5), and if it loses 124 MW of this generation, then this is equivalent to its net load increasing to 138 MW (note from row 2 the current load is 14 MW). Noting that the maximum net load that bus-2 can support is 138 MW, this loss of PV generation would lead to instability. The first row gives the ratio of maximum net load to PV generation at each bus. For any bus, when this ratio becomes less than 1, it means that the loss of renewable generation can trigger instability as mentioned previously for bus-2. Some situations that may lead to a loss of significant PV include: • Loss of a large PV farm connected to the bulk system. Under certain situations, a rapid cloud cover across a large area where bulk PV farms are located could create a similar severe impact.

•
System frequency events causing PV inverters to trip. This is true for both the bulk system and distributed PV. This was reported by North American Electric Reliability Corporation (NERC) for both the Blue Cut fire (2016) and Canyon fire (2017) events. In both cases, a fault on the transmission line resulted in an extreme system frequency [31,32], causing PV inverters to enter the momentary cessation mode. Moreover, IEEE 1547-2018 also directs inverters to enter momentary cessation if certain system frequency limits are violated.
Similarly, for a line, the amount of additional loading that leads to instability can be computed by comparing the stability margin at the operating point against the stability margin at the maximum loading and using the sensitivity σ (MW change in loading per 1% change in stability margin). Note that, under Sensitivity in the right-hand side of the heat map, Line-1 will go unstable if its load is increased by 51.3 MW, which is calculated by multiplying the Sensitivity factor by the Stability Margin. Note line 14 in the heat map of Figure 7, which is the weakest line in the system. Line 14 has a sensitivity (i.e., σ) of 0.01 MW/percent change in stability margin. Multiplying this sensitivity by the current stability margin of Line 14, 47.8, shows that if the flow through line 14 increases by 0.47 MW, the system will go unstable.
In the last column, the MW change causing instability is expressed as the percentage of line capacity. For some lines, such as line-8, this change is almost half of the line capacity. On the other hand, for line-14, an increase of 0.4% of line capacity would cause instability. Now let us consider for the base case how far each load bus location is from instability. This is done by noting the load at each load location in the base case, and the case where loading is increased for that location until the system becomes unstable. Here, this is termed as the MW distance from

Calculating MW Distance from Instability at Each Load Bus Location
Now let us consider for the base case how far each load bus location is from instability. This is done by noting the load at each load location in the base case, and the case where loading is increased for that location until the system becomes unstable. Here, this is termed as the MW distance from instability. The computed values are shown in Table 4. Note from Table 4 that Load Bus-2 is the weakest bus in the system (i.e., it supports the smallest increase in load before instability occurs), and this correlates with the 15% loading for Load Bus-2 noted in Figure 7. Also note from Table 4 that Load Bus-4 is the strongest bus in the system. To evaluate the impact of a generator failure on the voltage stability heat map, a generator in the vicinity of Load Bus-3 is taken out of service. The rated generation of this generator is 50 MW. Following the same steps as in the previous section, a new voltage stability heat map is obtained as shown in Figure 8.
From Figure 8, the generator failure has a negative impact on the voltage stability margins and the percent loading for five of the buses. The largest bus percent loading reduction was 5% at Load Bus-3, which is the bus closest to the generator that was taken out of service.
Due to the generator failure, the MW distance from instability is also reduced for several load buses. As shown in Figure 9, this reduction is at a maximum for load bus-3.

Conclusions
The use of a matrix-free, highly efficient GTA-based power flow is demonstrated for analyzing voltage stability of integrated T&D systems in the context of PV generation variability. Using sanitized models of real-world utility systems, it is shown that as the loading is increased on the system, stability margins decrease due to a decrease in voltages. However, as distributed PVs start generating, stability margins increase. This increase is different for different PV penetration levels. The relationship between PV penetration change and stability margin change is found to be nonlinear. For the examples considered here, the change in the stability margin of a transmission line is approximately linearly related to the change in the loading of the transmission line for variations in Figure 9. Impact of generator failure on MW distances.

Conclusions
The use of a matrix-free, highly efficient GTA-based power flow is demonstrated for analyzing voltage stability of integrated T&D systems in the context of PV generation variability. Using sanitized models of real-world utility systems, it is shown that as the loading is increased on the system, stability margins decrease due to a decrease in voltages. However, as distributed PVs start generating, stability margins increase. This increase is different for different PV penetration levels. The relationship between PV penetration change and stability margin change is found to be non-linear. For the examples considered here, the change in the stability margin of a transmission line is approximately linearly related to the change in the loading of the transmission line for variations in the stability margin from 70-100%. The slope of the linear relation provides the change in line loading per percent change in stability margin.
A voltage-stability heat map was introduced to demonstrate the identification of weak lines and buses in the system. By choosing lines with the least stability margins for the uniform increase in loading case, and by selecting high connectivity buses associated with the low stability margin lines, a two-dimensional heat map was constructed. This provided useful insights into system voltage stability, such as the MW distance of a line to an instability or the loss of distributed generation below a bus that leads to an instability. The impact of a generation failure on the voltage stability heat map was illustrated with an example. System operators can use this methodology for planning and real-time operation of their system.