Grid-scale Fluctuations and Forecast Error in Wind Power
GGrid-scale Fluctuations and Forecast Error in Wind Power
G. Bel, C. P. Connaughton, M. Toots, and M. M. Bandi Department of Solar Energy and Environmental Physics, Blaustein Institutes for Desert Research,Ben-Gurion University of the Negev, Sede Boqer Campus 84990, Israel Centre for Complexity Science, University of Warwick, Coventry, CV4 7AL, UK Collective Interactions Unit, OIST Graduate University,1919-1 Tancha, Onna-son, Okinawa, Japan 904-0495 (Dated: September 5, 2018)The fluctuations in wind power entering an electrical grid (Irish grid) were analyzed and foundto exhibit correlated fluctuations with a self-similar structure, a signature of large-scale correlationsin atmospheric turbulence. The statistical structure of temporal correlations for fluctuations ingenerated and forecast time series was used to quantify two types of forecast error: a timescale error( e τ ) that quantifies the deviations between the high frequency components of the forecast and thegenerated time series, and a scaling error ( e ζ ) that quantifies the degree to which the models fail topredict temporal correlations in the fluctuations of the generated power. With no a priori knowledgeof the forecast models, we suggest a simple memory kernel that reduces both the timescale error( e τ ) and the scaling error ( e ζ ). INTRODUCTION
Renewable power generation, unlike conventionalpower, exhibits variability owing to natural fluctuationsin the energy source [1]. Wind power, in particular,shares spectral features of the turbulent wind from whichit derives energy [2–4]. This variability in power outputadds a cost to renewable power [5, 6] that is absent in con-ventional power sources. Whereas distributed wind farmsare expected to smooth the fluctuations [7], power enter-ing the electrical grid still exhibits large amplitude fluctu-ations [1]. Large ramps in power fluctuations present thepossibility of grid destabilization [8] and blackout, a con-stant source of concern for system operators [7, 9]. Thisrisk increases the cost of operating reserves [10] neededon standby to return a grid back to operation in the eventof failure. Naturally, forecast models constitute essentialtools in estimating the magnitude of fluctuations before-hand and in planning for the optimal operating reservesrequired on call. Yet, no standards for forecast accuracycurrently exist [11].Extant works on wind power forecast error, rangingfrom the turbine to the grid scale, focus on modeling theforecast error distribution [12–17]. Since a probabilitydistribution is time-independent, it contains no informa-tion on temporal error variations. Several studies haveconsidered the dependence of the mean and the varianceof the error on the duration for which the power is pre-dicted (ranging from minutes to hours) [18, 19]. Otherworks have considered the different distributions of er-rors for mean power over different durations [14, 15, 20].However, none of these studies account for the fluctuationcorrelations of atmospheric turbulence [21] transferred tothe generated power in the analysis of forecast error norfor the temporal correlations of the errors.Whereas power fluctuations at the scale of an individ-ual turbine [3] and a wind farm [2, 4] have been shown to exhibit self-similar scaling, such fluctuations from in-dividual wind farms are expected to smooth out beforethey enter the electrical grid. Using data from the Irishgrid operator EIRGRID [22], we show that wind powerentering the grid exhibits correlated fluctuations with aself-similar structure. Such scaling points to large-scalecorrelations in atmospheric turbulence influencing the ag-gregate wind power entering the grid.In this article, we exploit these correlations at the gridscale and draw upon the Statistical Theory of Hydrody-namic Turbulence to quantify two types of forecast er-ror. The first is a timescale error ( e τ ) that quantifies thetimescales over which the forecast models fail to predicthigh-frequency power fluctuations. This timescale errorsets a bound on the numerical resolution of forecast mod-els and would already be known to system operators whoown and run the forecast models. However, details of themodels are not available to potential customers in energyspot markets [23] who could use this information to fac-tor in the risk associated with a non-supply of promisedwind power by producers. The second type of error wequantify is a scaling error ( e ζ ) that establishes a differ-ence in the self-similar scaling of fluctuations as observedfor actual generated power vis `a vis the power that wasforecast to be generated. This error could be potentiallyuseful to model developers, and if such an error resultsfrom large-scale correlations in atmospheric turbulence,incorporating them into models is not subject to limi-tations arising from numerical resolution. Having estab-lished the errors, we then employ a simple memory kernelupon the forecast time series and show that the errors canbe easily reduced with a minimal computational cost.Two raw time series are provided by EIRGRID: thewind power generated nationwide across all Ireland en-tering the grid p g ( t ), and the power forecast by EIR-GRID’s models p f ( t ) for the same period. The time seriessampled at 15 minute intervals span a five-year period a r X i v : . [ phy s i c s . d a t a - a n ] M a r FIG. 1: (color online) a) Raw time series (for 10 days) ofthe generated power p g ( t ) (black open circles), forecast power p f ( t ) (red open squares), and the instantaneous forecast error p d ( t ) (blue open triangles) in megawatts (MW). Every thirddata point is plotted for easy visibility. b) The probabilitydensity function of the raw instantaneous forecast error Π( p d )(solid black circles) has exponentially decaying, fat tails rela-tive to a Gaussian distribution (solid black line) of the samemean and standard deviation as Π( p d ). (2009 − p g ( t ) and forecastpower p f ( t ), and their instantaneous difference p d ( t ) ≡ p f ( t ) − p g ( t ), which we define as the instantaneous fore-cast error, are shown in Fig. 1a for a 10-day period,permitting a few immediate qualitative observations.Firstly, p g ( t ) exhibits correlated fluctuations. Secondly, p f ( t ), while closely following p g ( t ), misses the high fre-quency (relative to the sampling rate of time series) com-ponents. The instantaneous forecast error p d ( t ) exhibitscorrelated fluctuations with a coefficient of variation–standard deviation/mean = 116 . / . ∼
5, implyinglarge magnitude fluctuations in p d ( t ) (i.e., a broad dis-tribution). DISTRIBUTION OF FORECAST ERROR
As a point of comparison with prior works [12–17], wenote that the probability density function (PDF) of theraw instantaneous forecast error Π( p d ) (fig. 1b) exhibitsfat exponential tails that decay slower than a Gaussianfunction of the same mean and standard deviation asΠ( p d ). Indeed, a PDF with exponential tails may beexpected for reasons detailed in the following. Being ascalar product of an instantaneous force ( (cid:126)f ( t )) and ve-locity ( (cid:126)v ( t )), the statistics of temporal variation in powerare determined by the product of two random variables p ( t ) ≡ (cid:126)v ( t ) · (cid:126)f ( t ) [26]. The statistics of the product( Z = XY ) of two normally distributed random variables( X and Y ) was first studied by C. Craig [27] (henceforthreferred to as Craig’s-XY distribution). An asymptoticanalysis reveals that Craig’s-XY distribution is logarith-mically singular about zero with exponentially decayingtails, and its asymmetry (skewness) depends on the in-stantaneous cross-correlation between X and Y [28, 29].Whereas this asymptotic analysis is possible only when X and Y are Gaussian, the structure of Craig’s XY-distribution itself is more generally observed, even when X and/or Y are non-Gaussian [28].The basic structure of Craig’s-XY distribution is ex-pected for the PDF of power p ( t ) (e.g., compare Fig. 2c in[3] with Fig. 2 in [29]) and its estimation error δp ( t ). Sup-pose for a single wind turbine, the errors in estimating orforecasting velocity and force are δ(cid:126)v ( t ) and δ (cid:126)f ( t ), respec-tively, then p ( t ) + δp ( t ) ≡ ( (cid:126)v ( t ) + δ(cid:126)v ( t )) · ( (cid:126)f ( t ) + δ (cid:126)f ( t )).Expanding the RHS permits decomposition into power p ( t ) = (cid:126)v ( t ) · (cid:126)f ( t ) and its error δp ( t ) = (cid:126)v ( t ) · δ (cid:126)f ( t ) + δ(cid:126)v ( t ) · (cid:126)f ( t ) + δ(cid:126)v ( t ) · δ (cid:126)f ( t ) . (1)Consequently, δp ( t ) is a random variable whose statis-tics are determined by the sum of three terms (Eq. 1),each being a product of two random variables. One ex-pects δp ( t ) to exhibit features of Craig’s-XY distribu-tion irrespective of whether or not (cid:126)v , δ(cid:126)v , (cid:126)f , and δ (cid:126)f areGaussian. In fact, given that the velocity distribution ofatmospheric turbulence is known to follow the Weibulldistribution [30, 31], a numerical approach may becomenecessary.The power forecast error statistics can be easily scaledfrom the turbine to the grid scale. If M wind farms feedpower to the grid, and the i th farm has N i turbines, thecumulative error in estimating wind power then equalsthe instantaneous forecast error p d ( t ) for the grid, shownin Fig. 1a, and is given by: p d ( t ) = Σ Mi =1 Σ N i j =1 δp j ( t ) . (2)Summing the power error statistics over all turbines(across all farms) causes an averaging of fluctuations,starting with the most probable ones that occur aroundzero, thus smoothing the logarithmic singularity. Allthese generic features are readily observed for Π( p d ) inFig. 1b. Beyond contributing to the extant literature[12–17], the structure of Π( p d ) provides no useful infor-mation for our analysis and will not be discussed fur-ther. Understanding temporal variability (fluctuations)and uncertainty (error) requires analysis of the temporalevolution of the distributions, their moments and multi-point temporal correlation functions. We therefore pro-ceed through a statistical analysis of the temporal cor-relations in the fluctuating time series for generated andforecast power. DATA ANALYSIS
The time series was analyzed in two stages, with trendsin the series being identified in the first stage, followedby an analysis of the fluctuations around the trends inthe second stage. Trend removal permits a focus on sys-tematic differences between p g ( t ) and p f ( t ), ignoring dif-ferences due to new wind farms and seasonal variabilityof the wind power. Trend identification was performedsuch that the cross-correlation between the generated andforecast power trends was maximal. We used the fastFourier transform (FFT) for each of the series and definedthe trends by inverting the transform using only the fre-quencies with maximal amplitudes. The number of max-imal amplitudes was set by the requirement of the high-est cross-correlation between p g ( t ) and p f ( t ). Keepingthe zero frequency (to preserve the signal mean) and fivemore frequencies resulted in a peak cross-correlation of0 . P G ( t ), forecast power by P F ( t ) and theirinstantaneous difference by P D ( t ) ≡ P F ( t ) − P G ( t ).The characteristic fluctuation timescales for the de-trended time series were first computed from their re-spective autocorrelation functions defined as: C X ( τ ) = ( P X ( t ) − P X )( P X ( t + τ ) − P X )( P X ( t ) − P X ) (3)where P X is a time-average subtracted from the signal(de-trending does not render a zero signal mean since the FIG. 2: (color online) a) The five-year trend for p g ( t ) (blacksolid line) and p f ( t ) (red dashed line) is subtracted from theraw time series in subsequent analysis. b) Log-linear scale:autocorrelation functions C G ( τ ) (open black circles), C F ( τ )(open red squares) and C D ( τ ) (open blue triangles) for P G ( t ), P F ( t ) and P D ( t ), respectively, exhibit exponential decorrela-tion with respective characteristic timescales obtained fromfit to data of τ G = 80 .
94 points (20.24 hours), τ F = 81 points(20.24 hours) and τ D = 25 .
86 points ( ∼ zero frequency component was preserved). The subscript X should be replaced with G for generated power, F forforecast power, and D for instantaneous forecast error,respectively. The three autocorrelation functions (fig. 2b)exhibit exponential decay for short times with a datafit following the functional form C X ( τ ) ∼ A X e − ( τ/τ X ) ,where A X (cid:39) .
0, owing to C X ( τ ) being normalized,and τ X represents the characteristic decorrelation timefor each time series, yielding τ G = 80 .
94 data points( ∼ .
24 hours) for generated power, τ F = 81 points(also ∼ .
24 hours) for forecast power, and τ D = 25 . ∼ . τ X ; deviations were ap-parent only for long-term behavior spanning a week (orlonger timescales) when the decorrelation had alreadyoccurred. The correlation time of high frequency fluctu-ations ( (cid:46)
20 hours) is much shorter than the slow varyingtrend (over months to years). Hence the de-trending pro-tocol (in particular, the number of maximal amplitudes)does not influence the analysis to follow–a fact verifiedand reported upon later.Autocorrelation functions for the generated ( C G ( τ ))and forecast( C F ( τ )) power exhibit nearly identical scal-ing and the same characteristic decay timescales ( τ G = τ F = 20 .
24 hours), suggesting the accurate capture ofcorrelations in generated power by the forecast models.Yet, the autocorrelation function C D ( τ ) for instanta-neous forecast error P D ( t ) informs us that some corre-lations are not captured. In particular, we qualitativelyknow that P F ( t ) misses the high frequency componentsof P G ( t ), and they end up in P D ( t ), thereby contribut-ing to its two-point correlator. This correlation deficitsuggests that higher order moments of the two-point cor-relator are necessary to capture the statistical structureof the missing fluctuations. TEMPORAL STRUCTURE FUNCTIONS
Statistical analysis of higher order correlations is awell-developed, mature tool within the Statistical The-ory of Hydrodynamic Turbulence in which higher or-der two-point correlators are studied through
StructureFunctions . Kolmogorov’s theory of 1941 (K41) [32] laysthe foundation for structure functions through the cel-ebrated “4/5th law”: S ( r ) ≡ (cid:104) (∆ v || ( r )) (cid:105) ≡ (cid:104) ( v || ( R + r ) − v || ( R )) (cid:105) = − εr , where the third moment of lon-gitudinal velocity differences ( (cid:104) (∆ v || ( r )) (cid:105) ) between twopoints spatially separated by a longitudinal distance r ,is proportional to the product of the average turbulentdissipation rate ( ε ) and the longitudinal spacing r [33].The n th order structure function encodes all cross-terms up to order n of the two-point correlator for a givenstationary signal. The physical relevance of structurefunctions may be appreciated by considering a stationary,fluctuating signal x ( t ) with zero mean. The differencebetween the two values of this signal taken time τ apart(∆ x ( τ ) ≡ x ( t + τ ) − x ( t )) is collected at various windows(of duration τ ) along the time series. ∆ x ( τ ) is thereforea random variable with statistics of its own, and the n thorder structure function defined as S n ( τ ) = (cid:104) (∆ x ( τ )) n (cid:105) is the n th moment for its PDF Π(∆ x ( τ )). The moment ���� FIG. 3: (color online) Structure functions of order n = 1 − S Gn ( τ ) and (b) forecast power S Fn ( τ )plotted versus τ in log-log scale exhibit self-similar scaling S Xn ( τ ) ∝ τ ζ Xn ( X is G for generated and F for forecast power).The scaling is robust for (a) the generated power over 1.4decades (40 time steps). (b) In contrast, for forecast power,the first- and second-order structure functions exhibit scalingup to τ = 40 time steps, but for n >
2, no scaling is observedfor τ ≤
10 time steps. Self-similar scaling is restored over alimited range of timescales 10 < τ < S n ( τ ) varies with the time difference τ between signals,and its scaling (if any) reveals temporal variations in thestatistical structure of fluctuations in the signal to the n th order.Tails of the PDF Π(∆ x ( τ )) exert themselves with theincreasing order n of the structure function, thus necessi-tating more data to resolve higher order structure func-tions. A weak test for resolving the n th order structurefunction involves splitting the time series into smallerwindows and testing for identical scaling on the trun-cated series. However, this test only assures stationarityof the statistics. A strong test for the ability to resolvethe n th order structure function requires that first, themoment’s integrand (∆ x ) n Π(∆ x ) → | ∆ x | → ∞ [34](required due to finiteness of data), and second, the PDFΠ(∆ x ) should decay faster than 1 / | ∆ x | n +1 for | ∆ x | → ∞ or else the integral (cid:82) (∆ x ) n Π(∆ x ) d x would diverge forlarge | ∆ x | [35] (test for existence of a PDF’s n th mo-ment). Whereas the two conditions are not indepen-dent, the second condition is theoretical and does notdepend upon the available statistics. When conductingdata analysis, even when the second condition is satis-fied, insufficient data can lead to noise and prevent theintegrand (∆ x ) n Π(∆ x ) from satisfactorily converging tozero. The first condition is therefore dependent on finite-ness of data. Based on both weak and strong tests, weconclude that the EIRGRID data can resolve structurefunctions up to order n = 12; however, only results up to n = 10 are presented.Since even-order structure functions take only positivevalues, they converge faster than ones with odd order. Toovercome this distinction between odd and even orders,we compute the n th order structure function of the abso-lute value of differences: S Xn ( τ ) ≡ (cid:104)| P X ( t + τ ) − P X ( t ) | n (cid:105) where subtraction of mean P X ( t + τ ) and P X ( t ) is as-sumed. While ensuring the same convergence rate foreven- and odd-order statistics, it also collates all datain the positive quadrant permitting easy visualization.Analysis of fractional-order structure functions allowsbetter testing for anomalous scaling [36]. Fractional-order structure functions are only defined for absolutevalues of signal differences [36]–another reason why wecalculate structure functions of absolute differences. Wealso calculated the structure functions of orders n =0 . − . . RESULTS
Figure 3 plots the structure functions of order n =1 −
10 (fractional-order structure functions are calcu-lated but not shown) for the absolute value of signaldifferences of the generated power | ∆( P G ( τ )) | (fig. 3a)and forecast power | ∆( P F ( τ )) | (Fig. 3b). Self-similaror power-law scaling is observed for the generated powerstructure functions over 1 . τ ≤ n = 1and 2. For n >
2, no scaling is observed for timescales τ ≤
10. The scaling is restored over a limited range oftimescales 10 < τ <
40 (0 . S Xn ( τ ) ∝ A Xn τ ζ Xn (4) where ζ Xn is the scaling exponent. For simple mono-fractal scaling, ζ Xn ∝ n . However, fluctuations with amulti-fractal character exhibit a nonlinear dependence ofthe scaling exponent ζ Xn with respect to n . Super- (sub-)linear variation of ζ Xn versus n implies temporal expan-sion (compression) of fluctuations [37]. Scaling exponentsfor all the structure functions were computed from thelog derivative, ζ Xn = d log( S Xn ( τ )) d log( τ ) , which provides a morereliable estimate of the exponent than a power-law fit[36, 38]. The pre-factor A Xn in Eq. 4 is subsequentlyobtained from fit to data. In Fig. 3, all the data (redsolid circles) were divided by A Xn such that all fits (solidblack lines) commence from both mantissa ( τ ) and ordi-nate ( S Xn ( τ )) at unity, for an easy comparison of ζ Xn withorder n . All the data in Fig. 3, Fig. 4a, and Fig. 5b there-fore follow the scaling relation: S Xn ( τ ) ∝ τ ζ Xn ( A Xn ≡ S Fn ( τ ) for n > τ ≤
10 con-firms the qualitative observation made in Fig. 1a thatforecast models do not capture high frequency fluctu-ations. More importantly, Fig. 3b ascribes a precisebound on the time ( τ = 10, 2.5 hours) out to which thehigh frequency fluctuations are missed. Finally, scalingpresence for S Fn ( τ ) , n = 1 , C G ( τ ) and C F ( τ ) and their identical characteristic decay times, τ G and τ F , observed in Fig. 2b. This is to be expectedon the grounds that the second-order structure function S ( τ ) ≡ (cid:104) (∆ x ( τ )) (cid:105) = (cid:104) x ( t + τ ) (cid:105) + (cid:104) x ( t ) (cid:105)− (cid:104) x ( t ) x ( t + τ ) (cid:105) shares a direct correspondence with the autocorrelationfunction where the cross-term is identical to the numer-ator of Eq. 3. The failure of S Fn ( τ ) for n > τ = 10 reveals one typeof forecast error in the models; we call this the timescaleerror e τ .Before proceeding to the second type of error aris-ing from scaling mismatch, we define the cross-structurefunction X F Gn ( τ ) ≡ (cid:104)| P F ( t + τ ) − P G ( t ) | n (cid:105) . X F Gn ( τ ) rep-resents n th order moments for the PDF of the relativemagnitude of fluctuations between P G ( t ) and P F ( t + τ ),and their cross-terms correspond to higher-order two-point cross-correlators between the generated and fore-cast power. This function is plotted in Fig. 4a. Again,we notice that scaling is absent at early times ( τ ≤ < τ < X F Gn ( τ ) exhibits no scaling for n = 1 and 2, un-like the forecast structure functions (Fig. 3b). Although S Fn ( τ ) exhibits scaling for order n = 1 and 2, its expo-nent ζ Fn (cid:54) = ζ Gn ; this scaling deficit is reflected in X F Gn ( τ )for n = 1 and 2. �� �� FIG. 4: (color online) a) Log-log scale: cross-structure func-tions X FGn ( τ ) versus τ (red solid circles) exhibit no scaling atearly times τ ≤
10, with scaling restored for 10 < τ < ζ Xn versus the order of struc-ture function n for generated G (red solid circles), forecast F (blue solid squares) and modified forecast M (black solid tri-angles) structure functions, and cross-structure functions F G (green solid inverted triangles), and their respective second-order polynomial fits: solid red line for ζ Gn , small dashed blueline for ζ Fn , medium dashed black line for ζ Mn and long dashedgreen line for ζ FGn . DISCUSSION
Having established the various structure functions, wenow consider the behavior of their scaling exponents ζ Xn ( X ≡ G for generated, F for forecast and F G for thecross-structure function). Figure 4b plots ζ Xn versusthe order n together with their polynomial fits to thequadratic order. ζ Gn = 10 − + 0 . n − . n scalesalmost linearly (mono-fractal) with a small, but mea- surable, quadratic deviation towards multi-fractal behav-ior. The exponent ζ Fn = 0 .
007 + 0 . n − . n exhibitsa slightly more pronounced quadratic deviation (multifractal behavior) relative to ζ Gn . On the other hand, ζ F Gn = 10 − + 0 . n − . n scales almost linearly with n , implying mono-fractal scaling.We now consider the measurement error for the afore-mentioned scalings. Firstly, given that all de-trendingprotocols suffer from an ad hoc choice of a de-trendingtimescale, we tested the scalings for dependence on thede-trending procedure by varying the number of max-imal amplitudes. Ignoring the condition for maximalcross-correlation between p g ( t ) and p f ( t ), the number ofmaximal amplitudes contributing to the trends was var-ied. The scalings were invariant up to the inclusion of 15maximal amplitudes into the trend, beyond which, co-efficients for the polynomial fits started varying in thesecond decimal place. Having ascertained the robustnessof our choice for the five maximal amplitudes at which thecross-correlation peaks, we focused on a second source ofscaling measurement error, namely statistical variability.Since the scalings are analyzed up to τ = 100 data points,the de-trended time series were split into eight indepen-dent windows (each with 21912 data points), and thestructure functions were re-computed for each window.The variation in the log derivative ( ζ Xn = d log( S Xn ( τ )) d log( τ ) )for the eight independent measurements was taken as thepossible scatter in the scaling estimation, thereby pro-viding a confidence interval for the polynomial fits. Thescatter was found to be ζ Xn ± .
01 in both the measuredvalue of ζ Xn and the corresponding polynomial fits (foreach of the polynomial coefficients) for each of the eightindependent datasets, revealing that the polynomial fitswere meaningful only to the linear order for ζ Gn and ζ F Gn .The quadratic-order polynomial coefficient for ζ Fn , de-spite being larger than the scatter of ± .
01, is not usefulowing to the fact that the corresponding quadratic termsfor ζ Gn and ζ F Gn are smaller than the scatter magnitude.Despite qualitatively observing a quadratic deviationfor ζ Xn in Fig. 4b, our inability to ascribe significance to itarises from the fact that the multi-fractal component (de-viation from linear scaling) of the scalings is miniscule.This is significant in light of several studies that havedemonstrated multi-fractal scaling for wind power fluc-tuations at the turbine [3, 4] and farm scales [39]. Turbu-lence theory traces the source of multi-fractal behavior tointermittent fluctuations that can arise from two sourcesin the atmospheric context. The first, known as inter-nal intermittency, occurs at the small scales of turbulentflow. These intermittent fluctuations would be naturallyreflected in the power generated at the turbine and farmscales. However, when adding together power generatedby geographically distant wind farms, internal intermit-tency should smooth out [40] since it is a small-scale ef-fect and cannot extend across geographically distributedwind farms. Furthermore, the sampling interval (15 min-utes) for EIRGRID data is not expected to resolve anyeffects that may arise from internal intermittency, whichoccur at much shorter timescales (high frequencies).The second source of intermittency, known as externalintermittency, occurs at the edge of any free-stream [41]and arises in the atmospheric context due to coupling be-tween the atmospheric boundary layer turbulence and aco-moving weather system [21]. External intermittency,which can be experienced in the form of wind gusts, is ofgreater relevance in the present analysis as it can bothcorrelate distributed farms through the weather systemand occur at timescales longer than the 15-minute sam-pling interval for EIRGRID data. The nearly fractal scal-ing of ζ Gn informs us that both internal and external in-termittency are being smoothed to the point of renderinggrid-level power fluctuations almost mono-fractal.The self-similar scaling of S Gn ( τ ) over several hoursdoes strongly point to the influence of large-scale turbu-lent structures on power fluctuations at the grid level.The 20-hour characteristic decorrelation time ( τ G ) forgenerated power in Fig. 2c, if taken as the large eddyturnover time of atmospheric turbulence, also lends cre-dence to such an argument. Finally, independent proofin support of this argument also comes from Katzensteinet al. [40] who show that an individual wind farm ex-hibits f − / ( f being the frequency) scaling for the windpower spectrum (equivalent to τ / scaling of the second-order structure function in the time domain). However,as wind power from various farms is summed, the spec-trum steepens (please see Fig. 3 in [40]). Such spectralsteepening can be clearly attributed to the smoothingof high frequency (short timescale) fluctuations corre-sponding to small eddies. But the low frequency (longtimescale) fluctuations corresponding to large-scale ed-dies lose no power spectral density, clearly indicating theinfluence of large-scale turbulent structures.We finally consider the forecast error due to the scalingmismatch. We define the scaling error as e ζ ≡ ζ Fn − ζ Gn .Under this definition, if the time series for forecast andgenerated power were identical, then S Gn ( τ ) ≡ S Fn ( τ ),implying ζ Gn ≡ ζ Fn , and therefore e ζ = 0. Another typicalcase arises if forecast models fail completely, resulting ina flat time series with no fluctuations, ζ Fn = 0 resultingin an error e ζ = − ζ Gn . Using the polynomial fits for ζ Xn (see Fig. 4b) to linear order, we obtain e ζ = (7 × − + 0 . n ) − (10 − + 0 . n ) = − .
003 + 0 . n . Thiscan be cross-validated against the difference ζ Gn − ζ F Gn =(10 − +0 . n ) − (10 − +0 . n ) = 0 . n . Since ζ Xn → n →
0, the 0th order term falling within the scatter maybe taken to be zero. Both estimates of error are identicalin linear order ( e ζ = 0 . n ).The analysis thus far demonstrates the importance oftemporal correlations in wind power and their role in es-timating forecast errors. It is reasonable to ask whetherthis knowledge could help improve the forecast time se- ���� FIG. 5: (color online) a) Log-linear scale: γ opt versus order ofstructure function shows no improvement for n < n ≥ γ opt at n = 4. b) Log-log scale: Structure functions S Mn ( τ )versus τ (red solid circles) for the modified forecast time seriesshow considerable improvement over their counterparts S Fn ( τ )in Fig. 3b. ries, despite having no knowledge of the models em-ployed. In particular, to capture the short-term corre-lations missed by the forecast, we introduce a modifiedforecast that is based on the original forecast, convolutedwith an exponentially decaying memory kernel derivedfrom the generated power time series. The modified fore-cast power is given by P M ( t ) = (cid:82) t P F ( τ ) e − γ ( t − τ ) dτ .The memory duration (1 /γ ) was chosen so as to min-imize the relative difference between the structure func-tions of the generated and forecast power. As expected(as shown earlier, the low-order structure functions of thegenerated and forecast power are very similar), we foundthat the optimal γ varies with the order of the structurefunction. For n <
4, the memory-modified forecast showsno improvement in the agreement between S Gn and S Fn .For n ≥
4, the modified forecast exhibits better agree-ment with the structure functions of the generated poweras shown in Fig. 5b. The optimal γ ( γ opt ) was found tobe γ ≈ .
06 and γ ≈ .
37, as shown in Fig. 5a, plot-ted in log-linear scale to show the variation in γ opt for n ≥
4. The simple scheme, suggested here, not only triesto rectify the timescale error e τ , but also attempts tostatistically align the temporal correlations by improv-ing the scaling error e ζ .As is apparent from Fig. 5b, the structure functions( S Mn ( τ ) ≡ (cid:104)| ∆ P M ( τ ) | n (cid:105) ) for modified forecast time seriesare substantially improved over their unmodified coun-terpart (fig. 3b). Firstly, scalings are restored at highfrequencies ( τ ≤ ζ Mn =0 .
01 + 0 . n − . n . To linear order, the scaling er-ror e ζ = ζ Mn − ζ Gn = 0 . n − . n = 0 . n , a considerableimprovement over the original forecast time series. Be-ing computationally inexpensive, and given that spinningand non-spinning reserves must act within 10 minutesof failure, with replacement reserves acting within 20-60minutes, there are tangible benefits to incorporating sucha memory kernel into models to monitor instabilities inreal-time. Furthermore, it might be possible to improvethe forecast models using different parameterizations ofthe regional climate models or weather models. SUMMARY
In summary, wind power exhibits significant tempo-ral correlations even at the grid level, where fluctuationsare expected to average out [7] as power is fed fromgeographically distributed wind farms. Previous stud-ies show that the temporal correlations of the wind areessential to studying wind-generated large-scale oceancurrents [42]; a similar appreciation of large-scale cor-relations in atmospheric turbulence within the contextof wind power is called for. Fluctuations, albeit pos-ing a problem to system operators, possess a statisticalstructure through temporal correlations, which could beexploited to quantitatively analyze the error in forecastmodels. The technique proposed here is only limited bythe sampling rate of the time series. Beyond potentiallyserving as a standard for quantifying wind-power forecastaccuracy, it could have applications for any renewable en-ergy source with temporally correlated fluctuations pos-sessing a statistical structure.MT and MMB were supported by the OIST Gradu-ate University with subsidy funding from the CabinetOffice, Government of Japan. CPC was hosted by OISTGraduate University while performing this work. GB wassupported through the European Union Seventh Frame- work Programme (FP7/2007-2013) under grant number[293825]. The authors gratefully acknowledge EIRGRIDfor permission to use their data and N. Ouellette for sci-entific discussions. [1] D. J. C. MacKay,
Sustainable Energy - Without the HotAir (UIT Cambridge Ltd., Cambridge, UK, 2009).[2] J. Apt, J. Power Sources. , 369 (2007).[3] P. Milan, M. W´’achter, and J. Peinke, Phys. Rev. Lett. , 138701 (2013).[4] R. Calif, F. G. Schmitt, and Y. Huang, Geophys. Res.Abstracts , 15443 (2014).[5] C. Lueken, G. E. Cohen, and J. Apt, Environ. Sci. Tech-nol. , 9761 (2012).[6] W. Katzenstein and J. Apt, Energy Policy , 233(2012).[7] R. Wiser, Z. Yang, M. Hand, O. Hohmeyer, D. Infield,P. H. Jensen, V. Nikolaev, M. O’Malley, G. Sinden, andA. Zervos, Wind Energy, In IPCC Special Report on Re-newable Energy Sources and Climate Change Mitigation (Cambridge University Press, Cambridge, UK and NewYork, USA, 2011).[8] J. O. G. Tande, Appl. Energy , 395 (2000).[9] M. H. Albadi and E. F. El-Saadany, Electr. Pow. Syst.Res. , 627 (2010).[10] A. Fabbri, T. Gomez San Roman, R. Abbad, andV. H. M. Quezada, IEEE Trans. Power Syst. , 1440(2005).[11] A. Costa, A. Crespo, J. Navarro, G. Lizcano, H. Madsen,and E. Feitosa, Renew. Sustain. Energy Rev. , 1725(2008).[12] R. Doherty and M. O’Malley, IEEE Trans. Power Syst. , 587 (2005).[13] H. Bludszuweit, J. A. Dominguez-Navarro, and A. Llom-bart, IEEE Trans. Power Syst. , 983 (2008).[14] B.-M. Hodge and M. Milligan, NREL Report No.:NREL/CP-5500-50614 (2011).[15] B. M. Hodge, E. G. Ela, and M. Milligan, Wind Engi-neering , 509 (2012).[16] B. M. Hodge, D. Lew, M. Milligan, H. Holttinen, E. Sil-lanp¨a¨a, S .and G´omez-L´azaro, R. Scharff, L. S´’oder, X. G.Lars´en, G. Giebel, D. Flynn, et al., Wind power fore-casting error distributions: An international compari-son (National Renewable Energy Laboratory, Tech. Rep,2012).[17] J. Wu, B. Zhang, Z. Li, Y. Chen, and X. Miao, Elec.Pow. and Energy Syst. , 100 (2014).[18] H. Madsen, P. Pinson, G. Kariniotakis, H. A. Nielsen,and T. S. Nielsen, Wind Engineering , 475 (2009).[19] M. Lange, J. Sol. Energy Eng. , 177 (2005).[20] Note that [14] suggests the Cauchy distribution for theerrors. However, this distribution is not suitable becauseall its moments are undefined [43].[21] G. G. Katul and C. R. Chu, Phys. Fluids , 2480 (1994).[22] .[23] C. Weber, Energy Policy , 3155 (2010).[24] .[25] This does not apply for Ireland since its grid is isolated from mainland Europe.[26] One can define power as p ( t ) ≡ A(cid:126)v ( t ) (where A is asuitable constant) and apply the same arguments withoutloss of generality.[27] C. Craig, Ann. Math. Stat. , 1 (1936).[28] M. M. Bandi and C. P. Connaughton, Phys. Rev. E. ,036318 (2008).[29] M. M. Bandi, S. G. Chumakov, and C. P. Connaughton,Phys. Rev. E. , 016309 (2009).[30] V. J. Seguro and T. W. Lambert, Wind Engineering andIndustrial Aerodynamics , 75 (2000).[31] A. H. Monahan, Journal of Climate , 497 (2006).[32] A. N. Kolmogorov, Dokl. Akad. Nauk SSSR , 16(1941).[33] U. Frisch, Turbulence: The legacy of A. N. Kolmogorov (Cambridge University Press, UK, 1995).[34] M. M. Bandi, W. I. Goldburg, J. R. Cressman, andA. Pumir, Phys. Rev. E. , 026308 (2006).[35] G. Samorodnitsky and M. S. Taqqu, Stable Non-Gaussian Random Processes (Chapman and Hill, New York, 1994).[36] S. Y. Chen, B. Dhruva, S. Kurien, K. R. Sreenivasan,and M. A. Taylor, J. Fluid. Mech. , 183 (2005).[37] B. Mandelbrot and R. L. Hudson,