Copula-based modeling of stochastic wind power in Europe and implications for the Swiss power grid

Applied Energy 96 (2012) 33–44

Contents lists available at SciVerse ScienceDirect

Applied Energy journal homepage: www.elsevier.com/locate/apenergy

Copula-based modeling of stochastic wind power in Europe and implications for the Swiss power grid Simeon Hagspiel a,b,⇑, Antonis Papaemannouil b, Matthias Schmid b, Göran Andersson a a b

ETH Zürich, Power Systems Laboratory, Physikstrasse 3, 8092 Zürich, Switzerland swissgrid ag, Dammstrasse 3, 5070 Frick, Switzerland

a r t i c l e

i n f o

Article history: Received 22 June 2011 Received in revised form 11 October 2011 Accepted 22 October 2011 Available online 23 November 2011 Keywords: Stochastic wind power simulations Monte carlo method Copula theory Swiss power grid Probabilistic N  1 analysis Probabilistic transmission lines overload

a b s t r a c t Large scale integration of wind energy poses new challenges to the European power system due to its stochastic nature and often remote location. In this paper a multivariate uncertainty analysis problem is formulated for the integration of stochastic wind energy in the European grid. By applying Copula theory a synthetic set of data is generated from scarce wind speed reanalysis data in order to achieve the increased sample size for the subsequent Monte Carlo simulation. In the presented case study, European wind power samples are generated from the modeled stochastic process. Under the precondition of a modeled perfect market environment, wind power impacts dispatch decisions and therefore leads to alterations in power balances. Stochastic power balances are implemented in a detailed model of the European electricity network, based on the generated samples. Finally, a Monte Carlo method is used to determine power flows and contingencies in the system. An indicator is elaborated in order to analyze risk of overloading and to prioritize necessary grid reinforcements. Implications for the Swiss power grid are investigated in detail, revealing that the current system is significantly put at risk in certain areas by the further integration of wind power in Europe. It is the first time that the results of a probabilistic model for wind energy are further deployed within a power system analysis of the interconnected European grid. The method presented in this paper allows to account for stochastic wind energy in a load flow analysis and to evaluate deterministic indicators, such as the N  1 criterion, on a probabilistic basis. Thus, it constitutes an important extension of deterministic models. Ó 2011 Elsevier Ltd. All rights reserved.

1. Introduction 1.1. Motivation The composition of European power generation facilities is currently undergoing significant structural changes. Generation technologies deploying renewable energy sources now constitute the largest share of newly incorporated power plants in the European network [1]. Within the renewable energy sector, wind power is by far the most dynamically growing market, followed by solar power [2–4]. This trend is among other reasons enabled by the wind power’s ability to help meeting emission targets, mature level of technology, favorable distribution of wind resources throughout Europe and low production costs [5–8]. Analyzing the impact of increasing wind power penetration is thus

⇑ Corresponding author at: ETH Zürich, Power Systems Laboratory, Physikstrasse 3, 8092 Zürich, Switzerland. E-mail addresses: [email protected] (S. Hagspiel), Antonios.Papaemma [email protected] (A. Papaemannouil), [email protected] (M. Schmid), [email protected] (G. Andersson). 0306-2619/$ - see front matter Ó 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.apenergy.2011.10.039

of particular importance. Against this background, this paper specifically focuses on the integration of wind power into the European power system while disregarding the further development of other generation technologies. Wind power is by nature non-dispatchable because the primary energy source cannot be controlled [9]. Instead, it follows a stochastic distribution pattern, dependent on wind resource and wind turbine technology.1 Moreover, wind turbines are usually not built, where electricity is needed, but where the best conditions for this type of energy exist. These characteristics of wind power make its proper integration into the existing power systems difficult, especially regarding power balance, voltage stability and transmission network limitations. Thus, wind power constitutes a new challenge, especially for the highly meshed European interconnected grid. Several projects have addressed the physical as well as the commercial integration of wind energy in a pan-European framework [6,7,10–12]. However, they performed deterministic case studies for the integration of wind energy instead of applying a probabilis1 Note that reducing power output is possible. However, wind turbines are normally operated such that produced energy is maximized.

34

S. Hagspiel et al. / Applied Energy 96 (2012) 33–44

tic model which would much more realistically cope with the wind power’s characteristics [13,14]. On the other hand, a number of studies developed methodologies to probabilistically model wind power penetrating power systems [15–20]. Most of these publications include case studies, but for comparably small networks (such as test systems or national grids). A study focused on the Swiss power system as part of the European interconnected network considering stochastic wind power and probabilistic power flows around Europe has so far not been published. This paper intends to fill this gap by stochastically modeling the wind power Europe-wide and by investigating the effects on the Swiss power grid. 1.2. Structure of the paper The paper is structured as follows: In Section 2, different approaches for probabilistic power system analyses are discussed. The model to produce multivariate wind power distributions for various sites in Europe is provided in detail in Section 3, including a quality assessment of certain model components. In the case study presented in Section 4, wind power is first sampled from the modeled stochastic process, and then incorporated in a detailed model of the European power grid. A Monte Carlo method is used to determine changing power flows and contingencies in the system. Section 5 provides results for two different scenarios, namely wind power as by 2010 and 2030. Finally, relevant conclusions are duly drawn in Section 6. 2. Applied methodology The presented methodology considers a probabilistic system adequacy evaluation for several steady state grid conditions. Different methodologies can be applied for a probabilistic power flow analysis which are based on two general solution principles: either analytical or numerical [14,21]. The first step of both is to define the probability distribution functions (PDFs) of the stochastic input variables. When applying an analytical solution method, arithmetic operations, i.e. convolution techniques, have to be used to determine PDFs of system states and power flows. In this context two major problems arise, which is the non-linearity of the equations needed to calculate power flows on the one hand, and the dependence structure of the input variables on the other. In order to be able to overcome these difficulties, several assumptions are usually made [21]:    

Linearization of power flow equations. Total independent or linear-correlated power variables. Normal and discrete distributions. Constant network configuration and parameters.

The above list clearly indicates the drawbacks of the analytical solution method. Linearization results in low accuracy of the results when input variables are far from their corresponding means. As will be shown in the following, highly simplified correlation is not a feasible assumption when dealing with wind power. Furthermore wind speeds are not normally distributed. Several studies resolve the difficulty regarding dependence structures by making use of statistical moments and Cornish–Fisher expansion [22,23]. However, they model wind speed prediction errors instead of wind power infeeds in order to obtain normally distributed variables, and use linearized power flow equations. Other than that mentioned, a point estimate method can be used as presented in [24,19], solving the probabilistic load flow problem based on the first four central moments of each random variable and their spatial correlations.

In contrast to the analytical solution, the numerical approach adopts a Monte Carlo method to calculate the system states [15]. A deterministic power flow model is fed a large number of times with inputs of different combinations of nodal power values. The exact non-linear equations of the power flow problem can be calculated thus avoiding all the drawbacks listed for the analytical method and resulting in high accuracy. Furthermore, proved and tested software tools can be used to solve the power flows. However, these advantages usually come at the cost of large computation time. This paper presents a method to generate correlated European wind speed samples of arbitrary size by applying copula theory. The so generated data can subsequently be utilized in a Monte Carlo simulation. Hence, this approach constitutes a suitable and feasible option to analyze the effects of wind power penetrating the electricity network. The steps to be taken for the presented methodology are the following: 1. Acquire s wind speed samples, one for each geographic site of interest. Each sample consists of M data points, where M is usually small. 2. Analyze the data with respect to correlation and marginal distribution, i.e. calculate the correlation matrix RðR 2 Rss Þ and the density functions FX(x) = P(X 6 x). 3. Generate a synthetic set of data consisting of N data points using Copula theory. 4. Apply power curve models and nominal installed power to convert wind to electric power. Distributions for electric power infeed are thus obtained. 5. Use a power market model to determine new power balances for each data point k(k = 1, . . . , N). 6. Incorporate the new power balances in a model of the grid, i.e. generate a sample of N scenarios representing different grid situations. 7. Run a Monte Carlo simulation to determine power flows and/ or N  1 contingencies for each scenario k. 8. Evaluate the resulting distributions for the properties of interest (in this paper N  1 contingency). 9. Deduce the one-dimensional indicator v, which provides a basis for interpretation and prioritization. 3. Modeling wind power 3.1. Wind speed data Wind turbines convert wind energy into electrical energy, which is then fed into the power grid. As wind energy is kinetic energy of moving air masses, wind speed is the decisive factor for the generation of electric energy. The publicly available reanalysis database of the National Oceanic and Atmospheric Administration (NOAA) was deployed which contains reanalysis wind speed data at 10 m above ground level with a temporal resolution of 6 h [25]. The locations of the T62 Gaussian grid covering Europe are indicated by the yellow crosses in Fig. 1. It should be noticed that reanalysis data is produced at very coarse spatial and temporal resolutions, thus underestimating spatial and temporal variability and hardly covering extreme short-period peaks. This most surely leads to overestimating spatial correlation. However, the great advantage of this database is the good spatial coverage and consistency of contents [26]. It was thus decided that the data is suitable for the intended study. Data was extracted for 32 years (1979–2010) which is the maximal number of years available in this database. In a further step, wind speeds for two specific days per year were selected, one in January (winter-day) and one in July (summer-day). Hence, these

S. Hagspiel et al. / Applied Energy 96 (2012) 33–44

35

Fig. 1. Wind speed correlation for a reference point in northern Germany, winter and summer-day.

sets (consisting of 128 points each) represent typical wind speed conditions over Europe in winter and summer time. Note that the chosen periods overlap the instances for which the grid situations are analyzed in the case study presented in Section 4.2. 3.2. Correlation analysis As wind is a continuously moving source of energy, wind speeds for nearby locations are expected to be dependent of each other. By calculating Pearson’s linear correlation coefficient, which is a measure of the linear dependence of two stochastic variables X and Y, this assumption can be verified. It is calculated as in the following equation, where cov is the covariance of X and Y, and rX and rY are standard deviations

qX;Y ¼ corrðX; YÞ ¼

cov ðX; YÞ

rX rY

ð1Þ

From the linear dependence between multiple stochastic variables a correlation matrix R can be constructed containing all pairwise correlation coefficients (hence R 2 Rss ). R is characterized through symmetry and positive-definiteness, and all diagonal elements are by definition equal to 1. Correlation maps can be visual-

ized by extracting one specific line (or column, respectively) which corresponds to the location of interest. An example for a site in Northern Germany is shown in Fig. 1, for both the winter and the summer-day. For both cases, it can be seen that day-specific wind conditions show stronger correlation in east–west direction, which is a consequence of the major propagation direction of weather fronts in Europe. However, correlation of wind speeds is more pronounced on the winter-day. Further it can be seen that a large share of the total area is characterized by very low correlation with respect to the chosen reference point. Note that these findings hold true for other examined reference points, too.

3.3. Synthetic wind data generation For the final purpose of modeling stochastic wind power to be used in a Monte Carlo simulation it is necessary to increase the number of data points in order to obtain meaningful results. A synthetic set of data is thus to be generated representing the behavior of the original data with samples of arbitrary size. Wind speeds at multiple sites are represented by multivariate stochastic distributions. To generate such distributions with a variable number of samples, copula theory was identified to be a suit-

36

S. Hagspiel et al. / Applied Energy 96 (2012) 33–44

able method [9,15]. This method applies marginal distributions at each site together with the dependence structure of the distributions to generate a new set of data based on this information. We follow the approach described in [15] utilizing the wind speed data previously introduced. The first step is to generate uniformly distributed, interdependent distributions of arbitrary size using copula theory. Secondly, the inverse transform sampling method is applied to obtain the same marginal distributions as in the original data. 3.3.1. Random number generation and copula theory Generating uniformly distributions is easy to realize using pseudo-random number generation algorithms embedded in many software packages. However, the implementation of arbitrary multivariate distributions is a more difficult task. Copula theory is used to construct multivariate distributions such that different dependencies can be represented. The approach to formulating a multivariate distribution using a copula is based on the idea that a simple transformation can be made of each marginal variable in such a way that each transformed marginal variable has a uniform distribution (see next section). The dependence structure can be expressed as a multivariate distribution on the obtained uniforms. A copula essentially is a multivariate distribution on marginally uniform random variables. Mathematically, Sklar’s theorem is the foundation theorem for most applications of the copula. Sklar himself described copula as method to ‘‘join together one-dimensional distribution functions to form multivariate distribution functions’’ [27]. Definition 1 (Sklar’s theorem). : The random variables X and Y with cumulative distribution function FX and FY, are joint by copula C if their joint distribution can be written as:

F XY ðx; yÞ ¼ CðF X ðxÞ; F Y ðyÞÞ

ð2Þ

If FX and FY are continuous, then C is unique; otherwise the copula C is unique on the range of values of the marginal distributions. There are different families of copulas, for example Gaussian, Elliptical or Archimedian copulas. Several types of copula were applied to the data and tested graphically as suggested in [28]. The results indicated that a Gaussian copula model appropriately fits the data. The corresponding QQ-plots are shown at the end of this section (Fig. 3). In the context of modeling wind power distributions with copulas, it was shown in previous publications that the differences depending on the choice of copula family are rather small [9]. Given the results from graphical testing and encouraged by the fact that this type is comparably easy to implement for high-dimensional datasets, a Gaussian copula was applied in this paper. The Gaussian copula is constructed via Sklar’s theorem. C, as it was defined in Eq. (2), becomes for the Gaussian copula in the bivariate case:

C s ðu; v Þ ¼ Us ðU1 ðuÞ; U1 ðv ÞÞ

ð3Þ

with u, v 2 [0, 1] uniformly distributed variables. U is the univariate, and Us the bivariate standard normal distribution function:

  1 1 UðxÞ ¼ pffiffiffiffiffiffiffi exp  x2 2 2p  2  1 x  2sxy þ y2 Us ðx; yÞ ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffi exp  2 2ð1  s Þ 2p 1  s2

ð4Þ ð5Þ

For Gaussian copulas the rank correlation coefficient s can be directly calculated from the correlation coefficient as it was given in Eq. (1):

sðx; yÞ ¼ 2 sin

p 6

qðx; yÞ



ð6Þ

3.3.2. Marginal distributions For a multivariate distribution, i.e. a set of probability distributions that are linked with each other, the marginal distribution describes the probability distribution of a subset of random variables, while discarding the other random variables not contained in this subset. For example, in the dataset described in Section 3.1, a marginal distribution can be found for each site by analyzing the distribution of wind speed at this site only (without considering the other sites). The marginal distributions in turn are needed for the sampling algorithm. By applying the inverse transform sampling method the multivariate but so far uniform distributions can be transformed to their corresponding marginal distribution. Starting point for the inverse transform sampling method is the inverse probability integral transform which states that for a random variable (r.v.) X with an invertible cumulative distribution function (CDF) FX(x) = P(X 6 x), the r.v. FX(X) follows a uniform distribution on the interval [0, 1]. The proof of this statement is provided in the following equation. For: r 2 [0, 1]:

  h i 1 PðF X ðXÞ 6 rÞ ¼ P X 6 F 1 X ðrÞ ¼ F X F X ðrÞ ¼ r

ð7Þ

Thus, when taking a uniform distribution U in [0, 1], then it follows from FX(X) = U that F 1 X ðUÞ follows the distribution of X. It is important to notice that by applying the inverse CDF the dependence structure is not influenced. The sampling procedure is summarized as follows: 1. Find the stochastic dependence structure of the original data by calculating the correlation matrix R. 2. Generate random numbers Ui from the standard uniform distribution for each site i. This is done using a copula model in order to be able to include the stochastic dependence structure. 3. Apply the inverse cumulative distribution function F 1 X i ðU i Þ to each Ui. Three exemplary samples are illustrated in Fig. 2, where wind speeds of three specific locations are put in relation to each other. The 3000 points of the drawn sample are indicated by blue dots and the original data are represented by red crosses (128 points). A and B are close and wind speeds are strongly correlated, whereas C is distant from A and uncorrelated. Marginal distributions as well as correlations inherent in the original data are reflected in the sampled data. As explained above, we followed the approach suggested in [28] to graphically test different copula models, in our case the Gaussian. Fig. 3 shows the QQ-plots of the sampled and original wind speed data for the very same three locations depicted in Fig. 2. The points of each QQ-plot lie on the line, except for few outliers occurring for high wind speeds. The underlying distributions seem to be very similar. Overall it can be concluded that the chosen model adequately fits the original data. 3.4. Wind to electric power conversion Wind turbines convert the energy carried by the wind resource to electrical energy. The efficiency of this process is determined by the mechanical and electrical characteristics of the wind turbine and can be quantified by power curves. These curves indicate how much electric power will be measured at the output terminal when a certain wind speed occurs at hub height. In this study three classes of generation sites – Lowland, Highland and Offshore – are used for wind turbines respectively wind farms as suggested in [29]. By analyzing the corresponding geography, one class is assigned to each production site.

37

S. Hagspiel et al. / Applied Energy 96 (2012) 33–44

25

Theoretical Quantiles [m/s]

14

10 8 6

15 10 5 0

Synthetic data Reanalysis data 0

5

10

15

20

25

wind speed at location A [m/s]

25 Synthetic data Reanalysis data 20

15

20

25

5

0

5

10

15

Observed Wind Speeds [m/s]

15

10

25

5

0

5

10

15

20

25

wind speed at location A [m/s] Fig. 2. Scatter plots of three specific locations for both sample and original data.

20 15 10 5 0

Wind speeds are scaled to typical hub heights of modern wind turbines assuming a power law:

¼ v h0 ðh1 =h0 Þa

0

5

10

15

20

25

Observed Wind Speeds [m/s] Fig. 3. QQ-Plots of the sampled and original data for three specific locations. Accuracy is lost for high wind speeds.

ð8Þ

where h0 is the reference height, h1 the height of interest and a the shear exponent. Following [29], hub heights of 70 m were assumed for lowland and offshore, and 60m for upland turbines. a is 0.2 for lowland and 0.1 for upland and offshore installations. In general, every type of wind turbine is characterized by an own power curve. Power curves, as shown in Fig. 4, were defined as in [30] for each class and represent average state-of-the-art characteristics of a wind turbine in this sector, based on information regarding array efficiency, cut out wind speed, availability, electrical losses, etc. Furthermore, three additional power curves were constructed (dotted lines) that reflect the estimated future improvements of wind turbines until 2030. The range of nominal power output is expected to broaden for all classes thus increasing total energy gain, and hub heights increase to 100 m. The distributions of electric energy input, i.e. the probability that a certain amount of electricity is produced at a specific production site, is deduced by discrete convolution functions taking wind speed distribution and power curve as argument.

1

output power [ratio of nominal power]

1

10

10

0

0

vh

5

15

Theoretical Quantiles [m/s]

0

0

Observed Wind Speeds [m/s]

2

wind speed at location C [m/s]

20

4

Theoretical Quantiles [m/s]

wind speed at location B [m/s]

12

0.8

0.6

0.4

2010 offshore 2010 lowland 2010 upland 2030 offshore 2030 lowland 2030 upland

0.2

0 0

5

10

15

20

wind speed [m/s] Fig. 4. Power curve models.

25

30

35

38

S. Hagspiel et al. / Applied Energy 96 (2012) 33–44

For two discrete functions f ; g : D ! C, defined on D # Z, the discrete convolution is defined as:

ðf  gÞ ¼

X k2D

f ðkÞgðn  kÞ ¼

X

f ðn  kÞgðkÞ

ð9Þ

k2D

As shown in the lower graph of Fig. 5, the resulting distribution shows a large concentration of occurrences at zero and nominal output power. This is typical for the electricity production pattern of a wind turbine. However, the distribution clearly depends on wind resource as well as on turbine technology. It should be noticed that no forced outage rates (FOR) are included in the model, i.e. occurrences of zero output power in cases, where the wind turbine is taken out of operation (e.g. for maintenance work). However, it was shown in previous publications that those do not significantly impact the results of a probabilistic power flow analysis [31]. 3.5. Nominal wind power capacities in Europe In order to determine wind power capacities, the TradeWind consortium elaborated a detailed survey together with European National Wind Associations. Europe was subdivided into more than 100 zones for which the associations determined present and estimated future levels of nominal wind power. For the estimation of future trends, three scenarios were defined up to 2030, representing low, medium and high penetration levels, where medium is the most likely outcome and low and high are the lowest and highest credible outcomes [32]. Compared to 10.5% in

2010, these three scenarios constitute a 24.2%, 32.5% and 42.3% share of wind power with respect to total European capacity in 2030, given that all other capacities do not change. In Fig. 6 this data is reproduced in an illustrative way for the year 2010 and 2030 medium scenarios. In the modeling process, one location of the T62 Gaussian grid as it was presented in Fig. 1 is assigned to each zone. Consequently, wind power distributions in absolute terms can be calculated by simple multiplication. Similarly, one representative node of the European interconnected grid is assigned to each zone thus determining the point of infeed for wind power in this particular area. In reality wind parks will certainly be more spatially dispersed and additional nodes will serve as infeed for wind power. However, the aggregation into more than 100 zones maintains a sufficient level of detail, especially regarding the Swiss power grid which is far away from large wind farms. For people interested in similar studies it should nevertheless be noticed that a higher resolution might be necessary as soon as the focus shifts to other areas, where more wind power is built.

3.6. Validation of model components A procedure was developed to validate the following model components: reanalysis wind speed data, model for wind to electric power conversion and data on nominal wind power capacities. The 2008 wind speed time series was extracted from NOAA’s reanalysis database and entirely applied to the wind to power conversion model previously described. By multiplication with the nominal wind power capacities in each zone the time series of output power are obtained. A simple integration yields annual energy produced from wind. Done per country, the output can be compared to statistical data from Eurostat for the same time frame [1]. This comparison is shown in Fig. 7. Total deviation remains well below 10%, and deviation on national level is significant for a few countries only. Especially for Great Britain and Ireland the model overestimates wind energy. The impact of this error is certainly significant when aiming at regional studies for these countries. However, in the European power grid, Great Britain and Ireland are connected to the continent via DC link which has a limited transfer capacity of 2000 MW. Hence, for load flows in central Europe and Switzerland in particular, the impact of overestimated wind power in Ireland and Great Britain is limited. Overall it can be concluded that the tested model components are of adequate quality.

4. Case study for wind power in the European grid 4.1. Dispatch model

Fig. 5. Wind to electric power conversion via convolution.

After having generated the wind power samples using the method described in the previous sections a dispatch model determines, where the additional power produced from wind shall be compensated. This econometric model is built on merit order curves constructed from generation capacities and marginal costs in each country. At each moment, supply meets the inelastic demand thus determining the generation mix. In this context – due to very low marginal costs – wind power is assumed to always provide electricity to the grid. Further, the market aims at minimizing costs and therefore allows for trading activities across borders. In order to set up the dispatch algorithm several simplifying assumptions were made which are provided in the following list. A more sophisticated optimization model such as nodal pricing or other could also be utilized (for instance, [33,34] use nodal pricing for

39

S. Hagspiel et al. / Applied Energy 96 (2012) 33–44

Installed nominal wind power: 1000MW in 2010 (yellow) 2000MW in 2030 (blue) EU cumulative: 84GW in 2010 262GW in 2030

Fig. 6. Nominal installed wind power in European countries, 2010 and 2030 scenario.

Wind energy in 2008 [TWh]

50 Eurostat Statistics Results from Model

40 30 20 10

Au s Be tria lg Bu ium lg a C ria ro C a ze ch C tia R ypr ep us u D blic en m Es ark to Fi nia nl an F d G ra G er nce re m at an Br y ita G in re H ece un ga ry Ita La ly L Lu ith tvia xe ua m nia bo ur g N et M he al rla ta n N ds or w Po ay l Po and rtu g Ire al R lan om d an Se ia Sl rbi ov a Sl akia ov en i Sp a Sw ain Sw e itz de er n la nd

0

Fig. 7. Comparison of modeling results and statistics.

wind power market penetration). However, modeling the dispatch in great detail is not the focus of this study.  Perfect market environment: In a short term perspective, power plant owners bid their electric power at marginal costs. After having determined the total demand for electricity, cheapest bids are accepted and electricity is dispatched until demand is covered.  Marginal costs: Marginal costs of different types of power plants were evaluated in [35], where they take the values presented in Table 1.  Generation capacities: Generation capacities for different plant types in European countries are given ENTSO-E’s System Adequacy Forecast (SAF) [36]. The SAF not only publishes current generation capacities, but also presents two scenarios for the development until 2025 in 5 years intervals. This opens up the possibility to study different market

Table 1 Marginal costs of different electricity production technologies. Production technology

Marginal costs [Euro/MW h]

Wind Hydro Nuclear Lignite Hard coal Gas Renewables (w/o wind & hydro) Mixed oil/gas Oil Fossil non-attributable Not clearly identifiable

2 3 11 18 19 27 33 52 55 58 62

environments as they are estimated by European TSOs, both for a conservative as well as for a best estimate scenario. For

40

S. Hagspiel et al. / Applied Energy 96 (2012) 33–44

this paper generation capacities were assumed as published for 2010, adjusted by availability factors.  Net Transfer Capacity: Only limited capacities are available to transport electric energy from country A to B. These limitations are formalized by the Net Transfer Capacity (NTC) concept which is the maximum exchange program between two areas [37]. NTC values are taken as published by ENTSOE for Winter 2009–2010 and Summer 2010 [38,39]. It should be noted that only the commercial exchange is affected by applying NTC. In contrast, physical power flows are determined by grid configuration and the physics of the network (i.e. power flows along the least-resistance path).

4.2. Power Grid model The changes in power production and consumption determined within the dispatch model are integrated in a physical model of the European power grid. For the present study a network model is used which contains around 6000 nodes and 9000 branches including information on network topology, properties of lines and transformers, load and generators. Therefore, it covers the entire European interconnected power grid in great detail [40]. Snapshots were taken from operations data and mapped to the network model. Three instances in time were chosen representing different load situations that are crucial for the Swiss power grid:  Winter off-peak (3rd Wednesday in January 2011, 3:30am): Maximum import situation (approx. 3000 MW import).  Summer peak (3rd Wednesday in July 2010, 10:30am): maximum export situation (approx. 5000 MW export).  Winter peak (3rd Wednesday in January 2011, 10:30am): country balance of approx. zero MW. These three snapshots cover the whole range of the Swiss power exchange balance and represent the base case scenarios for the following probabilistic network analysis. Note that wind speed data was modeled for the typical weather conditions observed over Europe on the very same days. Further note that snapshot scenarios taken from operations data contain information on the operational status (on/off) of each element in the system. However, for the system adequacy evaluation it was ensured that all elements are put into operation, such that the full network potential can be tapped.

5. Results 5.1. Computational performance For the present paper Instant Security Probing of Electrical Networks/Interactive Power Flow Analyser (ISPEN/IPFA) is used to perform power flow and N  1 contingency analyses. To solve a given problem, ISPEN applies an AC Newton–Raphson algorithm combined with the matrix sparsity method [43]. The stochastic input is generated using Matlab. Matlab is also used to control the consecutive feeding of ISPEN with input such that a large number of scenarios can be evaluated automatically (a feature which is not implemented in ISPEN). After having optimized the communication channel between Matlab and ISPEN a performance analysis was done for the process. The results are presented in Fig. 8, where the evaluation time is illustrated in dependence of the sample size. As expected, computation time for calculating power flows (N case) as well as N  1 contingencies is linearly increasing at rates of approximately 2.14 s/scenario and 6.31 s/scenario, respectively. These results were obtained on a 2.26 GHz Pentium notebook. Initialization time is small and roughly equal for both types of analysis. Noticeably, the process can be parallelized. If k redundant processors are at hand, calculation times could be reduced by the same factor k (without considering initialization). It should be noted that changing the subset of elements defined for the N  1 analysis significantly influences computational effort since the procedure is iterative. At the same time, defining a subset large enough to cover all system elements potentially influencing the flows in the area of interest is crucial. Critical N  1 contingencies could otherwise be missed. 5.2. Implications for the Swiss power grid Two scenarios shall be discussed in this section, namely wind power as by 2010 and by 2030 (most likely outcome). The model described in the previous sections is used to generate wind power production samples of 3000 data points each, and to calculate effects on N  1 contingencies arising in the Swiss power grid. The wind 2010 scenario results are provided in Section 5.2.1, whereas Section 5.2.2 deals with the wind 2030 scenario. These results are then used in a further step to assess criticality of elements in the Swiss power system, as discussed in Section 5.2.3. Some remarks should be made regarding the 2030 scenario: Besides the nominal European wind power capacities, all other 2000

4.3. Power flow and contingency analysis

1600

Evaluation Time [sec]

Based on the equations of a conventional AC power flow analysis, a dedicated software tool is used to evaluate stochastic input scenarios with respect to power flows and N  1 contingencies. For the latter, a set of elements in the system is consecutively omitted, and the impact of this event on all other elements in the system determined (i.e. power flows for the changed network topology are recalculated). Computational effort is generally high for such an iterative approach, which gave reason for limiting the set of omitted elements to the Swiss power grid and some nearby elements in the neighboring countries [41]. N  1 analysis is of particular importance to TSOs since it is ENTSO-E’s key rule to operating the electricity grid in a secure way [42]. The results section therefore lays a focus on N  1 loading. Notably the model would also allow to analyze loading levels in undisturbed grid configurations (N loading), an option which would entail reduced computational effort.

Evaluation Time N Evaluation Time N−1 Least Square Regression N Least Square Regression N−1

1800

1400 1200 1000 800 600 400 200 0

0

50

100

150

200

250

300

Sample size [#] Fig. 8. Computational performance of the calculation.

350

41

S. Hagspiel et al. / Applied Energy 96 (2012) 33–44

5.2.3. Criticality assessment With the ultimate purpose of using loading distributions to assess and prioritize criticality of system elements and/or network expansion plans, the probability (i.e. risk) of overloading can be determined for each line and transformer. To this end, a criticality indicator v is defined as the product of probability that a line is overloaded in the N  1 case (pN1 > 100%), weighted by the maximum loading occurring on this element (maxN1) and summed over the three snapshots (SP, WP and WOP). By weighting with the maximum loading it can be assured that the above mentioned low probability – high risk events are well taken into account. v is calculated for each network element x:

vðxÞ ¼

X SP;WP;WOP

pN1>100% ðxÞ  maxN1 ðxÞ

ð10Þ

Element Loading [%]

100 80 60

001 002 003 004 005 006 007 008 009 010 011 012 013 014 015 016 017 018 019 020

System Element

115 110 105 100 95 90 85 80 75 021 022 014 023 024 025 026 027 028 029 016 002 030 017 018 004 005 031 032 003

70

System Element

140 120 100 80 60 40 20 033 017 018 022 012 013 034 029 035 036 023 028 021 031 037 038 010 015 039 002

5.2.2. Wind 2030 scenario Compared to the 2010 scenario a very similar set of elements is critical in the 2030 case, as shown in Fig. 10. However, loading levels are generally higher, with additional elements exceeding the 100% mark and values up to more than 180%. At the same time distribution ranges widen, indicating increased volatility. Kurtosis is greater than 3 for 76.2% of the distributions. Shapes are notably very distinct for the different elements. Specifically it is not the elements with the greatest mean that enclose the most extreme values. Hence, reducing the distributions to their means (often done when shapes are similar) would be a dangerous pitfall.

120

40

Element Loading [%]

5.2.1. Wind 2010 scenario The results of the probabilistic load flow analysis are stochastic distributions of N  1 loadings on each element. These are visualized by means of boxplots. On each box, the central mark is the median, the edges of the box are the 25th and the 75th percentiles, and the whiskers extend to the most extreme data points not considered outliers. Points are defined as outliers if they are larger than q2 + 1.5(q2  q1) or smaller than q1  1.5(q2  q1), where q1 and q2 are the 25th and 75th percentiles. These are plotted individually. Outliers are not considered as errors but kept for the analysis. Analyzing the shape of the distributions, it was found that for 65.3% kurtosis is greater than for normally distributed data (i.e. >3), indicating heavy-tailed distributions. Special attention should thus be given to low probability – high risk events. The results for the wind 2010 scenario are shown in Fig. 9 for the three cases Summer Peak (SP), Winter Peak (WP) and Winter Off-Peak (WOP). Elements are ranked by their median values. Note that for each case only the 20 worst case elements are presented (the Swiss power grid consists of more than 250 elements). System elements are numbered consecutively as they appear, and equally as they reappear.

140

Element Loading [%]

system properties are held constant. Specifically, we do not consider factors such as alterations in the electricity generation mix (except from wind power) or modifications in the grid infrastructure which would be out of the scope of this study. In the future, these and other factors will certainly be subject to change and therefore influence the results presented in this study. This is why scenario 2030 should not be seen as a prognosis. It rather serves as an extreme scenario highlighting the effects of increasing wind power penetration in the European power system. Noticeably there might be structural changes in stochastic wind speed characteristics, e.g. due to global warming. Such issues and their implications should be subject of future research. Based on analysis of the reanalysis data, introduced in Section 3.1, we assume wind speed characteristics to be equal to those observed during the last 30 years.

System Element Fig. 9. Distributions of N  1 contingencies in the Swiss power grid. Wind scenario 2010.

Fig. 11 depicts indicator v for the two wind power scenarios 2010 and 2030 including all elements in the Swiss electric network – unless v does not appear to be zero. The colors of the stacked bars allow to identify the snapshots (SP, WP, WOP) for which the overloading occurs. Note that the indicators v for the 2010 wind scenario, as shown in Fig. 11a, are directly deduced from the distributions presented in Fig. 9, as can be retraced by the consistent element numbering. In the case of the Swiss power grid it becomes clear that Summer Peak is the most critical system state under the influence of stochastic wind power, whereas criticality of winter peak and off-peak is lower (though not negligible). From the very number of listed elements as well as from the range of v it can be concluded that higher wind power penetration in Europe significantly impacts the Swiss power grid and its secure operation. In Fig. 12 critical elements are highlighted for the wind 2030 scenario in a map of the Swiss power grid to make the results more

42

S. Hagspiel et al. / Applied Energy 96 (2012) 33–44

120 100 80 60

019

020

018

017

016

014

015

013

012

011

010

009

007

008

006

005

004

003

002

001

40

System Element

System Element

Element Loading [%]

140

Element Loading [%]

115 110 105 100 95 90 85 80

022 002 004 003 001 006 005 011 007 010 008 009 021 012 013 033 015 034 035 036 019 029 042 043 044 045 046 028 018 017 023 031

75

Summer Peak Winter Off−Peak Winter Peak

0

0.5

1

1.5

2

Criticality Indicator χ [−] 003

032

031

005

004

018

017

030

002

016

029

028

027

026

025

024

023

014

022

021

70

System Element

120 100 80 60 40

002

039

015

010

038

037

031

021

028

023

036

035

029

034

013

012

022

018

017

033

20

System Element Fig. 10. Distributions of N  1 contingencies in the Swiss power grid. Wind scenario 2030.

comprehensible. Highlighting colors were chosen to emphasize different ranges of v. This graphical representation reveals that most critical elements are located at or close to Swiss borders. Several lines and transformers on the German border are heavily loaded mainly in winter when large imports from Germany penetrate the Swiss network. Furthermore, the two corridors to Austria have a high probability of being affected. Lines to Italy are the most critical, especially during summertime when export as well as transit flows significantly put the corresponding lines at risk. To France it is the relatively weak connection in the region of Geneva which responds to additional wind power with N  1 loadings largely exceeding the 100% threshold. 6. Conclusions This paper presented a methodology to stochastically model the integration of wind power in future power systems, with focus on

System Element

Element Loading [%]

140

002 005 022 013 012 001 008 003 004 007 009 010 021 033 018 017 011 034 028 029 036 035 023 032 047 048 031 048 049 043 050 014 051 026 052 053 054 019 038 039 055 056 057 058 059 037 060 061 025

Summer Peak Winter Off−Peak Winter Peak

0

0.5

1

1.5

2

Criticality Indicator χ [−] Fig. 11. Priority ranking based on the criticality indicator v for the wind 2010 and 2030 scenarios.

the European interconnected power grid. The methodology can be used in transmission planning, in integrated system adequacy forecasting, where generation and transmission systems interact with each other, as well as for operational purposes. A Gaussian copula was used to generate multivariate wind speed distributions reproducing both marginal distributions as well as stochastic dependence. The model was tested graphically

S. Hagspiel et al. / Applied Energy 96 (2012) 33–44

43

∋ ∋ ∋ ∋ Fig. 12. Critical elements in the Swiss power grid. Wind 2030 scenario.

by QQ-plots indicating that the model appropriately fits the reanalysis data. In a next step, distributions of wind speed were converted to distributions of wind power. The latter modeling component coupled with the 2008 time series for reanalysis wind speeds and European nominal wind power capacities yielded the electricity produced from wind turbines in 2008. National and European amounts were compared to Eurostat statistics, revealing that wind energy in the system was accurately modeled with deviations from measured values in the range of less than 10%. Through the coupling of an econometric and a physical model of the European interconnected system, stochastic power flows in the grid were calculated, based on perfect market behavior and triggered by stochastic European wind power previously sampled. The results showed that despite very low wind power capacity in Switzerland, the country is impacted by additional power flows resulting from a large number of wind turbines in other parts of Europe. The probability that NTC limits at the Swiss borders are commercially used up to their limits is high. Moreover, this probability is expected to increase when more wind power is integrated in the future. Thus, as physical flows were found to often deviate from scheduled ones, the secure operation of the grid is severely affected. The probabilistic approach of the model not only allowed to identify critical elements, but also to quantitatively assess the level (i.e. risk) of criticality on these elements. To this end, the onedimensional indicator v was defined based on N  1 contingency analysis results, measuring the criticality of loading conditions. v thus establishes a mean to prioritize expansion plans on a probabilistic basis. The results revealed that by integrating additional wind energy in the European power grid – as it is expected until 2030 – the number of critical elements as well as the level of N  1 loading in the Swiss power grid significantly increases. System elements located at or close to Swiss borders, especially towards Austria and Italy, are affected the most. Special attention should thus be given to these areas when it comes to operating procedures or grid expansion planning. The major drawbacks of the methodology are twofold: Firstly, computational effort is significant since a large number of power

flow calculations has to be performed subsequently in a Monte Carlo simulation. However, by applying copula theory the number of scenarios can be controlled, as shown in Section 3.3. Furthermore, the process could be parallelized on multiple processors. Using modern computers and appropriate software tools computation times can thus be brought to an acceptable level (see Section 5.1). Secondly, probabilistic modeling is not as straightforward as in the deterministic case. In this context, the paper presents a clear and structured approach to probabilistically model wind power, to implement the sampled data in a probabilistic system adequacy evaluation, and to interpret the results.

Acknowledgments The authors thank ETH Zürich and equally Swissgrid for having provided the institutional framework allowing to perform this study. Swissgrid’s network development group particularly supported the project.

References [1] Energy statistics. Eurostat; 2010a. . [2] Wind in power, 2009 European statistics. EWEA; 2010. [3] Wind power barometer. EurObserv’ER; 2010b. [4] IEA Wind Energy – Annual Report 2008. Internation Energy Agency (IEA); 2009. [5] Energiewirtschaftliche Planung für die Netzintegration von Windenergie in Deutschland an Land und Offshore bis zum Jahr 2020 (DENA-Netzstudie I). Tech. rep.; Deutsche Energie-Agentur GmbH; 2005. [6] Integration erneuerbarer Energien in die deutsche Stromversorgung im Zeitraum 2015–2020 mit Ausblick 2025 (DENA-Netzstudie II). Tech. rep.; Deutsche Energie-Agentur GmbH; 2010. [7] Hulle FV, Tande JO, Uhlen K, Warland L, Korps M, Meibom P, et al. Integrating wind – developing europe’s power market for the large-scale integration of wind power. Tech. rep.; TradeWind; 2009. [8] Leuthold F, Jeske T, Weigt H, von Hirschhausen C. When the wind blows over Europe – a simulation analysis and the impact of grid extensions. Electricity Markets Working Papers, 2009. [9] Leuthold F, Jeske T, Weigt H, von Hirschhausen C. When the wind blows over Europe – a simulation analysis and the impact of grid extensions. Electricity Markets Working Papers, 2009.

44

S. Hagspiel et al. / Applied Energy 96 (2012) 33–44

[10] Leuthold F, Jeske T, Weigt H, von Hirschhausen C. When the wind blows over Europe – a simulation analysis and the impact of grid extensions. Electricity Markets Working Papers, 2009. [11] Integrating wind power in the European power systems – prerequisites for successful and organic growth. Tech. rep.; UCTE; 2004a. [12] Wind Power in the UCTE Interconnected System. Tech. rep.; UCTE; 2004b. [13] CIGRE Working Group 601 of Study Committee C4. Review of the current status of tools and techniques for risk-based and probabilistic planning in power systems. Tech. rep.; CIGRE; 2010. [14] Buygi MO, Shanechi HM, Balzer G, Shahidehpour M. Transmission planning approaches in restructured power systems. In: IEEE Bologna power tech conference; 2003. [15] Papaefthymiou G, Kurowicka D. Using copulas for modeling stochastic dependence in power system uncertainty analysis. IEEE Trans Power Syst 2009;24(1):40–9. [16] Haghi HV, Bina MT, Golkar M, Moghaddas-Tafreshi S. Using copulas for analysis of large datasets in renewable distributed generation: PV and wind power integration in Iran. Renewable Energy 2010;35(9):1991–2000. [17] Klöckl B, Papaefthymiou G, Pinson P. Probabilistic tools for planning and operating power systems with distributed energy storage. Elektrotech Informationstechnik 2008;125(12):460–5. [18] Klöckl B, Papaefthymiou G. Multivariate time series models for studies on stochastic generators in power systems. Electric Power Syst Res 2010;80(3):265–76. [19] Morales J, Baringo L, Conejo A, Minguez R. Probabilistic power flow with correlated wind sources. IET Generat, Transm Distrib 2010;4(4):641–51. [20] Pinson P, Papaefthymiou G, Klöckl B, Nielsen H, Madsen H. From probabilistic forecasts to statistical scenarios of short-term wind power production. Wind Energy 2009;12(1):51–62. [21] Chen P, Chen Z, Bak-Jensen B. Probabilistic load flow: A review. In: The third international conference on electric utility deregulation and restructuring and power technologies, IEEE; 2008. p. 1586–91. [22] Usaola J. Probabilistic load flow in systems with wind generation. IET Generat, Transm Distrib 2009;3(12):1031–41. [23] Usaola J. Probabilistic load flow with correlated wind power injections. Electric Power Syst Res 2010;80(5):528–36. [24] Morales J, Perez-Ruiz J. Point estimate schemes to solve the probabilistic power flow. IEEE Trans Power Syst 2007;22(4):1594–601. [25] NCEP-DOE Reanalysis 2: Gaussian Grid. NOAA Earth System Research Laboratory; 2011. .

[26] Saha S, Moorthi S, Pan HL, Wu X, Wang J, Nadiga S, et al. The ncep climate forecast system reanalysis. Bull Am Meteorol Soc 2010;91(8):1015–57. [27] Sklar A. Fonctions de repartition a n dimensions et leurs marges. Publications de l’Institut de Statistique de L’Universite de Paris; 1959. [28] Genest C, Favre A. Everything you always wanted to know about copula but were afraid to ask. J Hydrol Eng 2007;12(4):347–68. [29] McLean J. Tradewind project, WP2.4 – characteristic wind speed time series. Tech. rep.; TradeWind; 2008a. [30] McLean J. Tradewind project, WP2.6 – equivalent wind power curves. Tech. rep.; TradeWind; 2008b. [31] Billinton R, Gao Y, Karki R. Application of a joint deterministic–probabilistic criterion to wind integrated bulk power system planning. IEEE Trans Power Syst 2010;25(3):1384–92. [32] van der Toorn G. Tradewind project, WP2.1 – wind power capacity data collection. Tech. rep.; TradeWind; 2007. [33] Morales J, Conejo A, Perez-Ruiz J. Simulating the impact of wind production on locational marginal prices. IEEE Trans Power Syst 2011;26(2):820–8. [34] Leuthold F, Weigt H, von Hirschhausen C. Efficient pricing for european electricity networks the theory of nodal pricing applied to feeding-in wind in Germany. Utilities Policy 2008;16(4):284–91. [35] Korpas M, Warland L, Tande JOG, Uhlen K, Purchala K, Wagemans S. Tradewind project, WP3.2 – grid modelling and power system data. Tech. rep.; TradeWind; 2007. [36] System adequacy forecast 2010–2025. Tech. rep.; ENTSO-E; 2010. [37] Technical guidelines for net transfer capacity determination. Tech. rep.; UCTE; 2004. [38] Indicative values for Net Transfer Capacities (NTC) in Europe. Winter 2009– 2010, working day, peak hours. ENTSO-E; 2009. [39] Indicative values for Net Transfer Capacities (NTC) in Europe. Summer 2010, working day, peak hours. ENTSO-E; 2010. [40] Kurzidem M. Analysis of flow-based market coupling in oligopolistic power markets. Ph.D. Thesis; ETH Zurich; 2010. [41] Hagspiel S. Impact of stochastic wind energy on switzerland’s power grid. Master’s Thesis; ETH Zurich; 2011. [42] Policy 5: Emergency operations. ENTSO-E; 2010. [43] ISPEN Power Flow Software homepage. 2010. .