Resilience of electricity grids against transmission line overloads under wind power injection at different nodes

A steadily increasing fraction of renewable energy sources for electricity production requires a better understanding of how stochastic power generation affects the stability of electricity grids. Here, we assess the resilience of an IEEE test grid against single transmission line overloads under wind power injection based on the dc power flow equations and a quasi-static grid response to wind fluctuations. Thereby we focus on the mutual influence of wind power generation at different nodes. We find that overload probabilities vary strongly between different pairs of nodes and become highly affected by spatial correlations of wind fluctuations. An unexpected behaviour is uncovered: for a large number of node pairs, increasing wind power injection at one node can increase the power threshold at the other node with respect to line overloads in the grid. We find that this seemingly paradoxical behaviour is related to the topological distance of the overloaded line from the shortest path connecting the wind nodes. In the considered test grid, it occurs for all node pairs, where the overloaded line belongs to the shortest path.

supplemented by power from conventional generators to keep control over the balancing of generated and consumed power. Therefore, in view of plans to substitute a considerable fraction of the total generated power by renewable energy 18 , it becomes important to study favourable embedding strategies of renewable energy sources into existing power grids 19,20 .
In this paper, we study how wind power feeding at different nodes of a power grid affects its stability against overloads of transmission lines. Such overloading can lead to overheating and line overload. Specifically we address the following question: If a given amount of conventionally generated power shall be replaced by wind power, where are the most favourable locations for wind farms if single line overloads shall be avoided? The methodology of our approach follows partly previous work 19 and is illustrated in Fig. 1. Wind power is injected at the generator nodes of an IEEE test grid by replacing nodes of conventional generators. Thereby, heterogeneities in the power production and consumption, as well as in the transmission line properties are taken into account, and we embed the wind power in a topological environment typical for a generator node. For the distribution of the fluctuating wind power we choose a Weibull distribution. This is motivated by empirical results for wind velocities and the so-called "power curve" 21 , which describes the average relation between wind velocity and power. A quasi-static response of the grid state with respect to the wind fluctuations is considered and the balancing of generated and consumed power is ensured by scaling the power input from the remaining (non-substituted) conventional generators. For our analysis, the IEEE RTS-96 test grid 22 is taken as an example.
After describing the methods more specifically, we first address the question how strongly the resilience against line overloads varies with the location of wind power injection, if exactly one conventional generator in the grid is replaced by a wind farm. It turns out that the grid resilience can be quite sensitive to the location of the injection node. For a wind farm with average power production of 200 MW, we find the overload probabilities to vary over more than two orders of magnitude. Nodes with highest overload probabilities have at most one generator node as neighbour and they have a comparatively low total capacity of their emanating transmission lines. We then study how simultaneous wind power input at different nodes affects the grid resilience against line overloads. Therefore, pairs of the conventional generators are replaced by wind farms with fluctuating power generation. We find that a higher wind power input at one of the injection nodes can increase the threshold power for line overload at the other node. This surprising behaviour is correlated with the distance of the overloaded line from the shortest path connecting the pair of wind nodes. Finally, we show that spatial correlations between wind power fluctuations at the two injection nodes need to be taken into account in order to identify the best pairs with lowest overload probabilities. We conclude the paper with a summary of the key results, their impact for applications and an outlook for further investigations.

Methods
Feasibility regions. The power flow in the test grid is treated based on a linearised version of the ac power flow equations, which is an often applied procedure in the electric engineering literature 23 . In this approximation, three simplifications are made: (i) the resistances of the transmission lines [jk] between nodes j and k are neglected in comparison to their reactances, implying that their properties can be fully described by for all pairs of nodes l and m that are not connected, i j at the nodes = … j N 1, , are considered to be constant and here set to one, and (iii) the voltage phase The IEEE RTS-96 test grid consists of 30 generator nodes (red circles, dark grey), 41 load nodes (yellow circles, bright grey) and 108 transmission lines. Wind power is injected by replacing conventional generators as indicated (blue/grey circles). If the injected wind power becomes too strong, line overload occurs, as indicated by the red (dark) overloaded line. The radii of the circles is drawn proportional to the total power generated/consumed at the respective nodes (as listed in the IEEE data set, i.e. before replacement of conventional generators), and the thickness of the transmission lines is proportional to their maximum capacity. The two paths of connecting lines coloured in yellow (bright grey) mark shortest paths between the two wind feeding nodes. differences between nodes connected by a transmission line are assumed to be much smaller than one, in reference units, the (real) powers P j ejected from nodes j are then connected to the voltage phase angles by the dc power flow equations In the IEEE test grid data, powers of consumers (demands) are listed for a typical situation, and yield the net power = + P P P j j j (g) ( d) at each node j. Due to the lossless transmission in the dc power flow, there is a balancing of total generated and consumed power ∑ = P 0 j j , which follows also from the symmetry of the b jk and Eq. (1). By choosing one voltage phase angle as the reference, e. g., θ N = 0, Eq. (1) are solved to express the remaining N − 1 phase angles θ j in terms of the N − 1 independent P j (from power balancing, , the transmission line [jk] is considered to be overloaded. In studying the stability of the grid under additional injection of wind energy, we follow previous work 19 and assume that the fluctuations of wind power occur on time scales short compared to that of load fluctuations and long compared to time scales needed for power adjustment of conventional generators. Accordingly, the consumed powers P j (d) are taken to be constant. For including wind energy, we replace n of the P j (g) by powers > g 0 l from wind farms, = … l n 1, , , and rescale the remaining ones with a common factor to ensure ∑ = P 0 j j . This rescaling can be viewed as a simple means to account for droop control or regulation response 19 . As for the values of the P jk [ ] max in Eq. (2), we use the short-time emergency ratings of the IEEE Reliability Test System 1996 (IEEE RTS-96). This test grid was developed for comparative and benchmark studies 22 and provides 15 tables of data. In Table 1 we list the entries in the IEEE RTS-96 data set, from which the relevant information for this work was extracted.
Given a replacement, the phase angles from the solutions of Eq. (1) become linear functions of the wind powers, , the non-substituted P j (g) and the set of susceptances of the transmission lines. Inserting these solutions into the conditions (2) yields, for each transmission line, feasibility regions in an n-dimensional Cartesian space  + n spanned by the wind powers. For the line [jk], this region is bounded by the two limiting planes α . Considering the set of all such limiting planes, there is a subset, which, together with the coordinate planes, confines a convex n-polytope around the origin of the wind power space, which is not intersected by any of the limiting planes. This polytope constitutes the feasibility region with respect to a line overload anywhere in the grid. If the n wind powers … g g ( , , ) n 1 lie inside this region, no overload occurs. Otherwise, at least one transmission line overloads.
In this work, we consider either one (n = 1) or two (n = 2) wind farms. Then the polytope becomes a line (n = 1) or polygon (n = 1). The line is given by a threshold value g i (c) for each possible wind node i, and the polygons are denoted as P ij for the possible pairs (i, j) of injection nodes.

Statistics of wind power fluctuations and line overload probabilities.
To capture the statistics of wind power fluctuations is a difficult task that requires a good model for the transformation of wind speed into wind power (as performed by wind mills) and a description of the wind velocity statistics in the turbulent flow of the atmosphere, which shows long-ranged temporal and spatial correlations. There is continuing progress in the modelling of these issues 13,16,[24][25][26] but this progress has not yet matured to a state of established standard models.
Here we base our description on empirical findings for the distribution of wind velocities and on the known average relation between wind speed and power characterised by the power curve 21,27 . As for the spatial correlations between wind powers at different nodes, we consider the two extremes of completely uncorrelated and completely correlated fluctuations. This allows us to gain insight into the importance of such correlations for the probability of line overloads.

Source location
Power demand P j (d)    To be specific, we take the Weibull distribution , as reported in the literature 28,29 . In contrast to the shape parameter, the scale parameter λ v depends significantly on the location of the wind farm, because locations with comparable  k 2 v and stronger mean wind speed must have larger λ v . This variation, however, is not relevant in our modelling approach, because we consider a situation, where a given mean amount of wind power ḡ (tot) shall be injected into the grid. Different λ v then amount to different wind farm sizes.
According to the power-curve, the power g of a wind farm is on average proportional to the cube of the wind speed v over the most relevant range of velocities 27 , where wind turbines operate. Taking ∝ g v 3 , the PDF ρ g ( ) of wind powers becomes a Weibull distribution as well with shape . If one single wind farm is included into the grid, the scale parameter λ λ = i will be fixed for injection node i by the condition x a y 0 1 This equation is solved numerically for λ i . Knowing λ i , we obtain the line overload probability for wind power feeding at node i.
In the case of two wind farms, we need to specify the joint probability density ρ g g ( , ) 2 1 2 for the powers g 1 and g 2 at the two wind farms. For a given pair of wind nodes (i, j), this gives the line overload probability where I g g ( , ) ij 1 2 is the indicator function of the feasibility region ij  , i. e., ; , ; , For completely correlated wind velocity fluctuations, i. e. equal wind speeds at the two wind farm locations, we can write λ = g g 1 1 0 and λ = g g 2 2 0 , where g 0 is a reference power, as, e. g., given by the transformation of speed into power by one wind mill (of the same type in the two farms). Different scale parameters λ λ ≠ 1 2 take into account that the farms at the two locations can have a different size. Hence, we have λ λ = g g 2 1 1 2 as a constraint, which implies that the conditional probability of power g 2 for given power The joint probability density then becomes In both the uncorrelated and correlated case, the two scale parameters λ λ ( , ) The minimisation is performed by using the MATLAB Interior Point algorithm 31 and corresponds to an optimisation of the wind farm sizes, if wind power with total mean amount ḡ ij (tot) shall be injected at the nodes i and j. Hence, when comparing line overload probabilities for different pairs of wind feeding nodes, this optimisation is always implicitly assumed. We finally note that the ratio λ λ / j i in the correlated case is given by the slope of the longest line between the origin and one of the other corner points of  ij .

Results
Grid resilience under wind power injection at a single node. As described in the Methods, a critical power g i (c) defines the feasibility region for wind power injection at a single node i: ] i (c) no overload occurs in the grid, while for > g g i (c) at least one transmission line overloads. Figure 2 shows the corresponding Scientific RepoRts | 7: 11562 | DOI:10.1038/s41598-017-11465-w overload probabilities Π i calculated from Eq. (7), where the λ i were determined from Eq. (5) for a wind farm size with an average power production for all pairs. This corresponds to a power production of a mid-sized wind farm and also to the typical power production of the individual conventional generators in the test grid. The overload probabilities Π i vary by more than two orders of magnitude dependent on the selected injection node, showing that the protection against line overload can be an important factor for optimising wind power integration into existing grids. We found two structural features of the nodes to be particularly relevant for the magnitude of the associated overload probabilities. First, a fluctuation of high wind power should be more easily compensated by the grid, if it can be distributed among many strong lines in the immediate neighbourhood of the injection node. Analysing this correlation, we found the logarithms of the overload probabilities to linearly correlate with the capacity-weighted node degree with a Pearson correlation coefficient of −0.922. Secondly, conventional generator nodes in the neighbourhood of the wind node can help to stabilise the grid against line overload, because their regulated response implies a decrease of their power supply in the environment of the wind node upon increase of injected wind power. The most vulnerable nodes with high Π i in Fig. 2 have at most one conventional generator as neighbour and a small capacity-weighted node degree. The capacity-weighted node degrees for these nodes are even very similar in values, leading to seemingly repeated motifs of the corresponding Π ij in Fig. 2. Grid resilience under wind power injection at two nodes. For wind power injection g g ( , ) 1 2 at two nodes i and j, the feasibility regions are polygons ij  : If ∈ g g ( , ) ij . Let us first consider the areas F ij of the polygons, which give a rough measure for overload probabilities in the case of uncorrelated powers injected at the two wind farms [cf. Eqs. (8) and (9)]. The examples in Fig. 3(a) show that these areas can be quite different for given topological distance L ij . The histogram of all F ij in Fig. 3(c) displays a broad unimodal distribution with the mode at about (800 MW) 2 . The large spread of areas implies strong variations of overload probabilities for uncorrelated wind power feeding (see below) and hence reflects the relevance of proper node selection already seen in the case of single-node injection.
Interesting information on correlation properties of the two wind nodes with respect to line overload is provided by the shapes of the polygons. If wind power injection at one node would not affect the critical power injection for line overload at the other node, the polygons become simple rectangles. For comparison, two of such virtual rectangles, ( c) are indicated by the dashed lines for the polygon drawn in large scale in Fig. 3(a). Here, the value g i (0) is the critical value for power injection at node i, if no power is generated at node j, and g j (0) has the analogous meaning. The value g i (0) differs from the critical value g i (c) for single node injection, because the latter refers to a situation, when there is a nonzero conventionally generated power at node j. With increasing topological distances L ij between the wind nodes, correlations decrease and the polygons tend to exhibit a more rectangular shape, see Fig. 3(a).
More detailed information on the correlations can be inferred by analysing the edges of the polygons. The two edges along the coordinate axes = g 0 1 and = g 0 2 are due to the constraints ≥ g 0 2 and ≥ g 0 1 , respectively, while the other edges are associated with overloads of transmission lines. If the wind powers g 1 and g 2 at the injection nodes i and j are driven out of the feasibility region due to high wind velocities at one or both injection nodes, one For slope angles near zero (nearly horizontal edge) and π ± ( /2) (nearly vertical edge), a change of power at one injection node has almost no effect on the threshold power for line overload at the other node (as it is strictly the case for the edges in the rectangles). A negative slope means that an increase of wind power at one of the injection nodes decreases the threshold power at the other injection node for the overload of line kl [ ]. We refer to this expectable behaviour as negative node coupling. Interestingly, we also see in Fig. 3(a) that an increase of wind power at one injection node can rise the threshold power for line overload at the other injection node, i.e. it makes the grid more stable. This seemingly paradoxical behaviour occurs for positive slope angles and we then say that the wind nodes exhibit a positive coupling with respect to the overload of line kl [ ]. In a way, this resembles another peculiar behaviour known as Braess paradox 32 , where the addition of a transmission line, or an increase of its capacity, makes the grid less stable or weakens the flow.
Histograms of slope angles m atan( ) occur only for sufficiently small L ij . This is shown in Fig. 4(a), where all angles m atan( )

ij kl
[ ] are plotted against the distances L ij . With increasing L ij , the angle distribution separates into two peaks around zero and π ± ( /2) (in a repeated scheme, the figure may be viewed as periodically continued along the slope angle axis). Overall, the couplings become small with increasing distance L ij . The reason for the possibility to obtain positive couplings is that flows are directional and accordingly power injection from the wind nodes can compensate each other along a transmission line. If the overloaded line is on the shortest path connecting the two farms, one would therefore expect a compensating effect to occur with higher probability. More generally, we expect the positive couplings to be the more likely the closer the overload line lies on the shortest path. To quantify this feature, we introduce the following measure for the distance of the link kl [ ] to the shortest path connecting nodes i and j: This measure corresponds to the excess length when subtracting the distance L ij from the length of the path connecting nodes i and j via link kl [ ], see the inset of Fig. 4 . In Fig. 4  Finally, we analyse the impact of spatial correlations between the wind powers at the injection nodes on the line overload probabilities Π ij . We choose the same average total wind power production of = g 200 MW ij (tot) as considered in Fig. 2 for the single node injection. This average total power is generated now by two wind farms, where the individual farm sizes are optimised to yield the lowest possible overload probability Π ij (see Methods). For comparing the uncorrelated with the correlated case, we present the overload probabilities in an array in Fig. 5(a), where for each possible pair i j ( , ) of injection nodes, Π ij is plotted in a colour scale (see the scale bar in the figure). In the upper triangle of the array above its diagonal, the Π ij are given for the uncorrelated case, and in the lower triangle for the correlated case.
As a first result from Fig. 5(a) we find that the smallest Π ij are by about two orders of magnitude smaller than in Fig. 2 for single-node injection, which demonstrates the advantage of decentralised wind power generation for the grid stability (for same average total wind power feeding). The injection nodes 4, 10, 16, 20, and 30, yielding the lowest line overload probability under single node injection, cf. Figure 2, give also comparatively low overload probabilities when paired with another injection node, both in the uncorrelated and correlated case. Corresponding rows and columns of the Π matrix in Fig. 5(a) remain in the red (dark grey) regime of overload probabilities in the range − − − 10 10 5 3 (with one exception for node number = i 20 in the correlated case, which shows higher Π ij when paired with nodes = j 21, 22, or 23). That pairs of nodes tend to give low Π ij if one of its nodes yields a low overload probability as a single wind power source, is a consequence of our optimisation of the wind farm sizes: If a "good node" i (Π i low) is paired with a "bad node" j (Π j high), Π ij can be avoided to become very high by increasing the farm size at node i.
These common features for the correlated and uncorrelated case, however, do not imply that the best pairs of injection nodes are the same for uncorrelated and correlated wind powers. The five pairs yielding lowest Π ij are listed in Table 2. While one of the nodes with numbers 4, 10, and 30 appears in all pairings for both the uncorrelated and correlated case, its "pairing node" is always different, i.e. none of the five best pairs in the uncorrelated case agrees with one of the five best pairs in the correlated case. This demonstrates the relevance of wind power correlations in the search for optimal wind feeding nodes. Pairs of the same rank in Table 2 have about 2-3 times lower overload probabilities for correlated wind powers. For completeness, we show in Fig. 5(b) the ranking of all node pairs with respect to grid resilience against line overload.

Discussion and Outlook
In this work we studied the line overload probabilities under fluctuating wind power injection at different nodes of an IEEE test grid in a model setup, where one or two conventional generators are replaced by wind feeding nodes with given total average power production. A quasi-static response of the power flow to wind fluctuations was assumed and we calculated this flow from the dc power flow equations under the constraint of total balance as in Fig. 2 (now generated by two wind farms). In each of the two arrays, the upper/lower triangle shows the results for uncorrelated/correlated wind powers at the two farm locations.  Table 2. Node pairs for wind power injection with the lowest line overload probabilities. The five most optimal pairs of nodes for wind power feeding yielding the lowest probabilities for line overload in the gird. Different optimal pairs are obtained for uncorrelated and correlated wind powers at the two wind farm locations.
between generated and consumed power. This was achieved by a corresponding rescaling of the controllable generators. To describe the distribution of wind speeds, we used the Weibull distribution with a typical shape parameter = k 2 reported in the literature. By resorting to the known average relation between wind speed and power, this was transformed into a Weibull distribution of wind power with shape parameter = k 2/3. For two injection nodes we considered the two limiting cases of uncorrelated and completely correlated wind powers at the wind farm locations to get insight into the role of large-scale spatial correlations of wind fluctuations. Scale parameters of the Weibull distributions were optimised by adjusting the wind farm sizes to yield lowest overload probabilities in the case of two-node injection (at given total wind power injection).
The main findings of our study can be summarised as follows: (i) The overload probabilities vary strongly with the location of injection nodes as a consequence of the heterogeneous grid properties, showing their importance for optimal integration of wind energy. A reduction of overload probabilities by about two orders of magnitude is seen if the same average amount of wind power is injected via two nodes rather than a single node. This gives an estimate of the benefit of decentralising wind energy integration for avoiding line overloads. (ii) An analysis of the structure of feasibility regions in the two-node injection case allows one to obtain valuable insight into couplings between injection nodes with respect to transmission line overloads, which are independent of the detailed wind statistics. In particular, many positive couplings exist, where an increase of wind power at one injection node increases the threshold power for line overload at the other node. It was found that node pairs with positive couplings need to be topologically close to each other and that overloaded lines for positively coupled injection nodes tend to lay along the shortest path between these nodes. The first feature can be understood from the decrease of coupling strength with injection node distance, and the second feature from the likelihood that flows emanating from the injection nodes mutually compensate each other along transmission lines close to the shortest path. (iii) Large-scale spatial correlations between wind fluctuations seem to be a relevant factor for optimal integration of wind power. We found them to change the best pairs of nodes yielding lowest overload probabilities, and for these pairs to reduce the risk of line overload by a factor of two to three.
Generally, it is an important goal to provide proper tools and measures for estimating risks of transmission line overloads under increasing integration of renewable energy sources into power grids, and to develop risk-minimising strategies for the design of new grid structures that are better suitable for transporting electric energy between a large number of small power sources. The methodology used in this study, where we followed previous work reported in ref. 19 is just one step towards these goals. An important issue to be clarified in future studies is how far the quasi-static approach can be considered to be valid. This question is intimately connected with the magnitude of various time scales, in particular those specifying the power response to a change in wind speed and the correlations of wind speed fluctuations. It needs to be checked also how far the primary stability control on the scale of seconds 11 can be effectively accounted for by a simple power rescaling of controllable generators. We have started to investigate these problems by conducting dynamical studies based on swing equations for the IEEE RTS-96 test grid, i.e. the same test grid as used in this study. When succeeding to identify a time scale of the validity of the quasi-static approach it should become possible to connect the calculated values of the line overload probabilities to time intervals, where one can expect a line overload to occur. In the present study we have refrained to pursue any attempts in this direction and regarded the calculated overload probabilities solely as a relative measure for comparison of different wind power injections. Moreover, we have not considered the specific function that specifies the power curve of one wind turbine or one wind farm. It is known that this average function has a saturation level, above the so-called rated speed 21 , which is specific for each wind turbine and needs to be taken into account for calculating the costs associated with maintenance. For estimating overload probabilities with the proposed methodology, the saturation effect can be considered in the transformation of the velocity distribution (4) to the power distribution. However, a significant effect on the overload probabilities is only to be expected if the derivation from the Weibull distribution are pronounced in the feasibility regions.
Apart from these necessary checks, we believe that our concept of positive and negative node couplings with respect to flow through transmission lines can be a valuable general tool to qualify the mutual influencing of power injections. The concept can be generalised to more than two power injections ( > n 2). In this case, the feasibility region for power injection at n given nodes is an n-dimensional polytope with − n ( 1)-dimensional hyperplanes confining it, and the set of coefficients β rs , = … kj lj have opposite signs, and negative coupling otherwise. One may ask then, for example, how strongly the character of node coupling changes with the addition of injection nodes and how far the combined effect of many nodes on the flow can be traced back to pairwise couplings.
The results of our work can help to reveal the influence brought by renewable energy sources on the electricity grid. For applications, it is also important to estimate costs and returns, and to develop optimisation strategies with respect to the expected gain. Such optimisations will depend also on specific objectives of an investor and hence should involve expertise from economists working in this field. To illustrate how the overload probabilities can enter corresponding optimisation problems, let us consider a simple approach, where we introduce two different kinds of costs, fixed costs, and maintenance costs that vary in time. The fixed costs λ c ( ) . The total costs when implementing renewable sources at nodes i and j can then be defined as  where τ is the time horizon for the investment. The term κ Π ij ij is the average rate of overloads, with κ ij a proportionality factor. This proportionality factor depends on the injection nodes due to the costs depending on the